%嚴格按照BP網絡計算公式來設計的一個matlab程序,對BP網絡進行了優化設計
%優化1：設計了yyy，即在o(k)計算公式時，當網絡進入平坦區時(<0.0001)學習率加大，  出來后學習率又還原 P68 3.30
%優化2：v(i,j)=v(i,j)+deltv(i,j)+a*dv(i,j); 參照了P67 3.29

clear all
clc
inputNums=3;              %輸入層節點
outputNums=3;             %輸出層節點
hideNums=10;              %隱層節點數
maxcount=20000;           %最大迭代次數
samplenum=3;              %一個計數器，無意義
precision=0.001;          %預設精度
yyy=1.3;                  %yyy是幫助網絡加速走出平坦區 

alpha=0.01;                %學習率設定值
a=0.5;                    %BP優化算法的一個設定值，對上組訓練的調整值按比例修改
error=zeros(1,maxcount+1);  %error數組初始化；目的是預分配內存空間
errorp=zeros(1,samplenum);%同上

v=rand(inputNums,hideNums);      %3*10;v初始化為一個3*10的隨機歸一矩陣; v表輸入層到隱層的權值
deltv=zeros(inputNums,hideNums); %3*10;內存空間預分配
dv=zeros(inputNums,hideNums);    %3*10;  


w=rand(hideNums,outputNums);     %10*3;同V
deltw=zeros(hideNums,outputNums);%10*3
dw=zeros(hideNums,outputNums);   %10*3


samplelist=[0.1323,0.323,-0.132;0.321,0.2434,0.456;-0.6546,-0.3242,0.3255];  %3*3;指定輸入值3*3(實為3個向量)
expectlist=[0.5435,0.422,-0.642;0.1,0.562,0.5675;-0.6464,-0.756,0.11];  %3*3;期望輸出值3*3(實為3個向量)，有導師的監督學習

count=1;
while (count<=maxcount)          %結束條件1迭代2000次
    c=1;
    while (c<=samplenum)
      for k=1:outputNums 
        d(k)=expectlist(c,k);    %獲得期望輸出的向量,d(1:3)表示一個期望向量內 的值
      end
    
      for i=1:inputNums
        x(i)=samplelist(c,i);    %獲得輸入的向量（數據）,x(1:3)表一個訓練向量
      end

%Forward();
     for j=1:hideNums
         net=0.0;
           for i=1:inputNums
               net=net+x(i)*v(i,j);%輸入層到隱層的加權和∑X(i)V(i)               3.11
           end
         y(j)=1/(1+exp(-net));     %輸出層處理f(x)=1/(1+exp(-x))單極性sigmiod函數 3.14
     end
     for k=1:outputNums
         net=0.0;
         for j=1:hideNums
             net=net+y(j)*w(j,k);
         end
         if count>=2&&error(count)-error(count+1)<=0.0001
            o(k)=1/(1+exp(-net)/yyy);
         else o(k)=1/(1+exp(-net));     %同上
         end
     end

%BpError(c)反饋/修改;
    errortmp=0.0;
    for k=1:outputNums
        errortmp=errortmp+(d(k)-o(k))^2;  %第一組訓練后的誤差計算
    end
    errorp(c)=0.5*errortmp;               %誤差E=∑(d(k)-o(k))^2 * 1/2      3.15
%end
      
      %Backward();
    for k=1:outputNums
        yitao(k)=(d(k)-o(k))*o(k)*(1-o(k)); %輸入層誤差偏導                  3.25a
    end
    for j=1:hideNums
        tem=0.0;
        for k=1:outputNums
            tem=tem+yitao(k)*w(j,k);         %為了求隱層偏導，而計算的∑
        end
         yitay(j)=tem*y(j)*(1-y(j));         %隱層偏導                      3.25b
     end

   %調整各層權值
   for j=1:hideNums
       for k=1:outputNums
           deltw(j,k)=alpha*yitao(k)*y(j);    %權值w的調整量deltw(已乘學習率)3.26a
           w(j,k)=w(j,k)+deltw(j,k)+a*dw(j,k);%權值調整，這里的dw=dletw(t-1)，實際是對BP算法的一個
           dw(j,k)=deltw(j,k);                %改進措施--增加動量項目的是提高訓練速度          
       end                                                                 %3.29
   end
   for i=1:inputNums
       for j=1:hideNums
           deltv(i,j)=alpha*yitay(j)*x(i);     %同上deltw
           v(i,j)=v(i,j)+deltv(i,j)+a*dv(i,j); 
           dv(i,j)=deltv(i,j);
       end
   end
    c=c+1;
    end%第二個while結束；表示一次BP訓練結束
    
    double tmp;
    tmp=0.0;
    for i=1:samplenum
      tmp=tmp+errorp(i)*errorp(i);%誤差求和
    end
    tmp=tmp/c;
    error(count)=sqrt(tmp);%誤差求均方根,即精度

    if (error(count)<precision)%另一個結束條件
      break;
     end
    count=count+1;%訓練次數加1
end%第一個while結束
error(maxcount+1)=error(maxcount);
p=1:count;
pp=p/50;
plot(pp,error(p),'-'); %顯示誤差