clc
clear
a=[ 26.00 27.50 26.63 27.00 27.50 27.25 27.25 28.50 27.88 28.23 28.45 28.23 27.90 28.00 28.00 27.60 27.90 27.75 28.30 28.65 28.65 28.50 29.00 28.80 26.60 27.75 26.80 26.75 27.25 27.00]';
x1=a(1:9);  y1=a(12);
x2=a(4:12); y2=a(15);
x3=a(7:15); y3=a(18);
x4=a(10:18); y4=a(21);
x5=a(13:21); y5=a(24);
x6=a(16:24); y6=a(27);
x7=a(19:27); y7=a(30);
%x8=a(22:30); y8=a(33);
%x9=a(25:33); y9=a(36);
%x10=a(28:36);y10=a(39);
%x11=a(31:39);y11=a(42);
%x=[x1 x2 x3 x4 x5 x6 x7 ];
y=[y1 y2 y3 y4 y5 y6 y7];
%%%%%%進(jìn)行輸入數(shù)據(jù)的歸一化
[x1n,meanx1,stdx1]=prestd(x1);
[x2n,meanx2,stdx2]=prestd(x2);
[x3n,meanx3,stdx3]=prestd(x3);
[x4n,meanx4,stdx4]=prestd(x4);
[x5n,meanx5,stdx5]=prestd(x5);
[x6n,meanx6,stdx6]=prestd(x6);
[x7n,meanx7,stdx7]=prestd(x7);
%[x8n,meanx8,stdx8]=prestd(x8);
%[x9n,meanx9,stdx9]=prestd(x9);
%[x10n,meanx10,stdx10]=prestd(x10);
%[x11n,meanx11,stdx11]=prestd(x11);
xn=[x1n x2n x3n x4n x5n x6n x7n];
I=size(xn,1);P=size(xn,2);
%%%%%%%%%%%%使目標(biāo)向量歸一化
yn=y./29;
%%%%%%%%%%%%給出學(xué)習(xí)率、動(dòng)量因子等參數(shù)
eta=0.02;%學(xué)習(xí)率取的適當(dāng)可以加快誤差的減小
aerfa=0.965;%動(dòng)量因子
max_perf_inc=1.04;
lr_dec=0.7;
lr_inc=1.05;c=0.5;
%假設(shè)小波函數(shù)節(jié)點(diǎn)數(shù)為：H個(gè)；樣本數(shù)為P;
%輸出節(jié)點(diǎn)數(shù)為:J個(gè)；輸入節(jié)點(diǎn)數(shù)為：I個(gè)；
H=2*I+1;
J=1;
%%%%%%%%%%%%%%%%%%%%初始化權(quán)系數(shù),閾值，小波參數(shù)；
whi=rand(I,H);
wjh=rand(H,J);
b1=rand(H,1);
b=rand(H,1);
a=rand(H,1);
b2=rand(J);
%%%%%%%%%%%保存的誤差；
for p=1:P
     Err_NetOut=[];
     flag(p)=1;count(p)=0;
     while flag(p)>0  %進(jìn)入小波神經(jīng)網(wǎng)絡(luò)的運(yùn)行過程
           flag(p)=0;
           count(p)=count(p)+1;
           xhp1(p)=0;
           %%%%%%%%%%%求隱層的輸入值
           for h=1:H
               for i=1:I
                   xhp1(p)=xhp1(p)+whi(i,h)*xn(i,p);%求權(quán)值和輸入值的乘積和
               end
               ixhp(h,1)=xhp1(p)+b1(h,1);%求隱層的輸入值
               xhp1(p)=0;
           end
           %%%%%%%%%%%求隱層的輸出值
           for h=1:H
               oxhp(h,1)=fai((ixhp(h,1)-b(h,1))/a(h,1));%將小波函數(shù)(mallet)作為隱層的傳遞函數(shù)求隱層的輸出值
           end
           %%%%%%%%%%%求輸出層的輸入值
           ixjp1=0;
           for h=1:H
               ixjp1=ixjp1+wjh(h,1)*oxhp(h,1);
           end
           ixjp=ixjp1+b2(1);%求輸出層的輸入值
           %%%%%%%%%%求輸出層的輸出值
           ixjp2=c.*ixjp;
           oxjp(p)=fnn(ixjp2);   %求輸出層的輸出（運(yùn)用的是sigmoid函數(shù)(fnn)加了一個(gè)伸縮因子c）
           %%%%%%%%%%%求誤差并保存
           wuchayy=1/2*sumsqr(oxjp(p)-yn(p));%求輸出值與目標(biāo)值差值的平方和
           Err_NetOut=[Err_NetOut wuchayy];%保存每次的誤差;
           %調(diào)節(jié)學(xué)習(xí)率eta==================================================
           if count(p)==1
              eta=eta;
           else 
                   if Err_NetOut(count(p))/Err_NetOut(count(p)-1)>max_perf_inc
                         eta=eta*lr_dec;
                   else eta=eta*lr_inc;
                   end
          end    
          %%%%%%%%%%%%%求網(wǎng)絡(luò)參數(shù)的梯度值
          detaj=-c*(oxjp(p)-yn(p))*oxjp(p)*(1-oxjp(p));%在誤差函數(shù)中對(duì)輸出層的輸入值求偏導(dǎo)并取負(fù)值
          for h=1:H
              detawjh(h,1)=eta*detaj*oxhp(h,1);%誤差函數(shù)對(duì)隱層和輸出層之間的權(quán)值的偏導(dǎo)
          end
          detab2=eta*detaj;%誤差函數(shù)對(duì)輸出層的閾值的偏導(dǎo)
          for h=1:H
               detah(h,1)=detaj*wjh(h,1)*diffai((ixhp(h,1)-b(h,1))/a(h,1))/a(h,1);%difmallet()函數(shù)是mallet（）函數(shù)的導(dǎo)數(shù)，該行是求誤差函數(shù)對(duì)隱層輸入值的偏導(dǎo)
          end
          for h=1:H
              for i=1:I
                  detawhi(i,h)=eta*detah(h,1)*xn(i,p);%求誤差函數(shù)對(duì)輸入層和隱層之間的權(quán)值的偏導(dǎo)
              end
          end
          detab1=eta*detah;
          detab=-eta*detah;%求誤差函數(shù)對(duì)小波函數(shù)中的平移因子的偏導(dǎo)
          for h=1:H
               detaa(h,1)=-eta*detah(h,1)*((ixhp(h,1)-b(h,1))/a(h,1));%求誤差函數(shù)對(duì)隱層的小波函數(shù)中的收縮因子的偏導(dǎo)
          end
          detac=ixjp*detaj/c
          %%%%%%%%%%%%%%%%%%引入動(dòng)量因子aerfa，修正各個(gè)系數(shù)==========================================
          wjh=wjh+(1+aerfa)*detawjh;
          whi=whi+(1+aerfa)*detawhi;
          a=a+(1+aerfa)*detaa;
          b=b+(1+aerfa)*detab;
          b1=b1+(1+aerfa)*detab1;
          b2=b2+(1+aerfa)*detab2;
          c=c+(1+aerfa)*detac;
          %%%%%%%%%%%%%%%%%%%將輸出恢復(fù)到原始范圍=======================================
          oxjpn(p)=oxjp(p)*29;
          %%%%%%%%%%%%%%%%%%網(wǎng)絡(luò)循環(huán)計(jì)算的判斷條件
          if Err_NetOut(end)>1e-8
               flag(p)=flag(p)+1;
          else
                figure;
                plot(Err_NetOut);
                title('誤差曲線x(p)');
                disp('目標(biāo)達(dá)到');
                disp(oxjpn(p)); 
                disp(count);
          end
          if count>2000
                figure;
                plot(Err_NetOut);
                title('誤差曲線');
                disp('目標(biāo)未達(dá)到');
                disp(oxjpn(p)); 
          end
         %break;  %跳出避免死循環(huán)
         disp(oxjpn(p));
     end
 end