%用于函數逼近的BP算法程序
%作者魏海坤，摘自《神經網絡結構設計的方法與理論》
% function main()
SamNum = 100;               %....訓練樣本數
TestSamNum = 101;           %....測試樣本數
HiddenUnitNum = 10;         %....隱節點數
InDim = 1;                  %....樣本輸入維數
OutDim = 1;                 %....樣本輸出維數
%........................根據目標函數獲得樣本輸入輸出
rand('state',sum(100*clock))
NoiseVar = 0.1;
Noise = NoiseVar * randn(1,SamNum);
SamIn = 8 * rand(1,SamNum) - 4;
SamOutNoNoise = 1.1 * (1 - SamIn + 2*SamIn.^2) .* exp(-SamIn.^2/2);
SamOut = SamOutNoNoise + Noise;

TestSamIn = -4:0.08:4;
TestSamOut = 1.1 * (1 - TestSamIn + 2*TestSamIn.^2) .* exp(-TestSamIn.^2/2);

MaxEpochs = 30000;            %....最大訓練次數
lr = 0.003;                   %....學習率
E0 = 0.5;                     %....目標誤差

W1 = 0.2 * rand(HiddenUnitNum,InDim) - 0.1;     %....輸入層到隱層的初始權值
B1 = 0.2 * rand(HiddenUnitNum,1) - 0.1;         %....隱節點初始偏移
W2 = 0.2 * rand(OutDim,HiddenUnitNum) - 0.1;    %....隱層到輸出層的初始權值
B2 = 0.2 * rand(OutDim,1) -0.1;                 %....輸出層初始偏移
%.........................................................
W1Ex = [W1 B1];                                 %....輸入層到隱層的初始權值擴展
W2Ex = [W2 B2];                                 %....隱層到輸出層的初始權值擴展

SamInEx = [SamIn; ones(1,SamNum)];             %....樣本輸入擴展
ErrHistory = [];                                %....用于記錄每次權值調整后的訓練誤差
for i = 1:MaxEpochs
    
    %....正向傳播計算網絡輸出
    HiddenOut = logsig( W1Ex * SamInEx );
    HiddenOutEx = [HiddenOut; ones(1,SamNum)];
    NetworkOut = W2Ex * HiddenOutEx;
    
    %....停止學習判斷
    Error = SamOut - NetworkOut;
    SSE = sumsqr(Error);
    
    %....記錄每次權值調整后的訓練誤差
    ErrHistory = [ErrHistory SSE];
    
    if SSE < E0
        break
    end
    
    %....計算反向傳播誤差
    Delta2 = Error;
    Delta1 = W2' * Delta2 .* HiddenOut .* ( 1-HiddenOut );
    
    %....計算權值調節量
    dW2Ex = Delta2 * HiddenOutEx';
    dW1Ex = Delta1 * SamInEx';
    
    %....權值調節
    W1Ex = W1Ex + lr * dW1Ex;
    W2Ex = W2Ex + lr * dW2Ex;
    
    %....分離隱層到輸出層的初始權值，以便后面用
    W2 = W2Ex(:,1:HiddenUnitNum);
end

figure
hold on
grid
plot(SamIn,SamOut,'r+');
plot(TestSamIn,TestSamOut,'b--');
xlabel('Input x');
ylabel('Output y');

%....顯示計算結果
W1 = W1Ex(:,1:InDim)
B1 = W1Ex(:,InDim+1)
W2
B2 = W2Ex(:,1+HiddenUnitNum)

%....測試
TestHiddenOut = logsig( W1 * TestSamIn + repmat( B1,1,TestSamNum ));
TestNNOut = W2 * TestHiddenOut + repmat( B2,1,TestSamNum );
plot(TestSamIn,TestNNOut,'g-');

%....繪制學習誤差曲線
figure
hold on
grid
[xx,Num] = size(ErrHistory);
plot(1:Num,ErrHistory,'k-');