%BP網絡的第一階段學習期（訓練加權系數wki,wij）
%初始化
lr=0.05;err_goal=0.0002;%lr為修正權植的學習速率；err_goal為期望誤差最小植
max_epoch=1000; a=0.9;%max_epoch為訓練的最大次數
Oi=0;Ok=0;%置隱含層和輸出層各神經元輸出初值為零
%提供一組訓練集和目標值（1輸入，1輸出）
pi=3.1416;b=1;%b是要逼近的函數中的一參數，可選1，2，4，8
X=-2:0.1:2;T=1+sin((b/4)*pi*X);
%初始化Wki,Wij(M為輸入節點j的數量；q為隱含層節點i的數量；L為輸出節點k的數量）
[M,N]=size(X);q=6;[L,N]=size(T);%N為訓練集對數量
Wij=rand(q,M);Wki=rand(L,q);
Wij0=zeros(size(Wij));Wki0=zeros(size(Wki));
for epoch=1:max_epoch
    %計算隱含層各神經元輸出
    NETi=Wij*X;
    for j=1:N
        for i=1:q
            Oi(i,j)=1/(1+exp(-NETi(i,j)));
        end
    end
    %計算輸出層各神經元輸出
    NETk=Wki*Oi;
    for i=1:N
        for k=1:L
            Ok(k,i)=NETk(k,i);
        end
    end
    %計算誤差函數
    E=((T-Ok)'*(T-Ok))/2;
    if(E<err_goal) break;end
    %調整輸出層加權系數
    deltak=2*Ok.*(1-Ok).*(T-Ok);
    W=Wki;
    Wki=Wki+lr*deltak*Oi';
    Wki0=W;
    %調整隱含層加權系數
    deltai=Oi.*(1-Oi).*(deltak'*Wki)';
    W=Wij;
    Wij=Wij+lr*deltai*X'+a*(Wij-Wij0);
    Wij0=W;
end
epoch%顯示計算次數
%BP網絡的第二階段工作期（根據訓練好的Wki,Wij和給定的輸入計算輸出）
X1=X;%給定輸入
%計算隱含層各神經元輸出
NETi=Wij*X1;
for j=1:N
    for i=1:q
        Oi(i,j)=1/(1+exp(-NETi(i,j)));
    end
end
%計算輸出層各神經元輸出
NETk=Wki*Oi;
for i=1:N
    for k=1:L
        Ok(k,i)=NETk(k,i);
    end
end
SSE=mae(E);
Ok,q,SSE%顯示網絡輸出層的輸出,隱含層節點i的數量，平均絕對誤差
plot(X,T,'+',X,Ok)