%這是****貨運(yùn)量****的神經(jīng)網(wǎng)絡(luò)公式計算程序

clear
clf reset;
figure(gcf);
%setfsize(500,200);
echo on
clc
%初始數(shù)據(jù)及進(jìn)行歸一化處理——————————————————————————————————
                         %初始數(shù)據(jù)是兩個，分別是客運(yùn)量和貨運(yùn)量
clc
primary=[12815       15543       19326       22864       26150       28468 ...
         30882       39375       45759       49589       52560       48726 ...
         51083       56495       62767       83606       92090      101370 ...
        107317      108654      111847      112872      116997      126007      128980;
          2690        2998        3012        3042        3616        3728 ...
          3988        9397       17680       19426       24128       24354 ...
         22879       24162       28957       36439       40593       45052 ...
         47400       45224       45338       45815       47151       55705       63532]'; 
         %矩陣轉(zhuǎn)置，第一列是客運(yùn)量，第二列是貨運(yùn)量。
   %歸一化
    maxium=max(primary);%求出兩列的最大值，客運(yùn)量最大值：maxium(1),貨運(yùn)量最大值:maxium(2).
    minium=min(primary);%求出兩列的最小值，客運(yùn)量最小值: minium(1),貨運(yùn)量最小值:minium(2).
    p1=0.7*(primary(:,1)-minium(1))./(maxium(1)-minium(1))+0.2;%代表歸一化后的客運(yùn)量列向量
    p2=0.7*(primary(:,2)-minium(2))./(maxium(2)-minium(2))+0.2;%代表歸一化后的貨運(yùn)量列向量
    %歸一化，將來計算仿真值的原值，需用上面的公式！！！！！
    
second=[p1';
        p2'];%將向量又轉(zhuǎn)置，第一行p1'代表客運(yùn)量的歸一化后的向量，第二行p2'代表貨運(yùn)量的歸一化后的向量。
    
 %初始數(shù)據(jù)及進(jìn)行歸一化處理結(jié)束！———————————————————————————————
 
 %開始劃分?jǐn)?shù)據(jù)，分為訓(xùn)練集和檢驗集——————————————————————————————
                     %總樣本量N=25,取序列長度n=5,滾動間隔為n0=2。
                                                               
                     % 貨運(yùn)量———————————本程序只劃分貨運(yùn)量數(shù)據(jù)集———————————————
                     for i=1:10
                        y1(i)=p2(2*i); 
                        y2(i)=p2(2*i+1);
                        y3(i)=p2(2*i+2);
                        y4(i)=p2(2*i+3);
                        y5(i)=p2(2*i+4);%y1等是時間序列行向量。
                        t(i)=p2(2*i+5);
                        %上為訓(xùn)練集，下為測試集——————————————
                        cs1(i)=p2(2*i-1); 
                        cs2(i)=p2(2*i);
                        cs3(i)=p2(2*i+1);
                        cs4(i)=p2(2*i+2);
                        cs5(i)=p2(2*i+3);%cs1等是測試集時間序列行向量。
                        tcs(i)=p2(2*i+4);
                        
                     end
                     y=[y1;
                        y2;
                        y3;
                        y4;
                        y5];%y是貨運(yùn)量的矩陣，y1,y2,y3,y4,y5為5次輸入，8個樣本。    
                     p=y;
                     t=t;%目標(biāo)輸出
                     
                     cs=[cs1;
                         cs2;
                         cs3;
                         cs4;
                         cs5]; %cs等是測試集。 
                      
                      
                      %預(yù)測2003年和2004年的貨運(yùn)量——————————————
                      p2003=[p2(21);
                             p2(22);
                             p2(23);
                             p2(24);
                             p2(25)];
                        
                      p2004=[p2(22);
                             p2(23);
                             p2(24);
                             p2(25)];

                     
                     
                     

  %數(shù)據(jù)劃分結(jié)束，分為訓(xùn)練集和檢驗集，訓(xùn)練集是p,t.—————————————————————
  %5個輸入，樣本為10個。
  %輸出為1個向量，樣本為10個
%pause
clc
plot(t,'+');
title('Training Vectors');
xlabel('Input Vector P');
ylabel('Target Vector T');
%確定輸入個數(shù)、隱含層神經(jīng)元、輸出層神經(jīng)元數(shù)目————————————————————————
clc
r=5;          %輸入個數(shù)
s1=8;         %第一層神經(jīng)元數(shù)
s2=5;         %第二層神經(jīng)元數(shù)
s3=1;         %第三層神經(jīng)元數(shù)
%確定各個層的權(quán)重和閥值———————————————————————————————————
[w1,b1,w2,b2,w3,b3]=initff(p,s1,'tansig',s2,'tansig',t,'purelin');
echo off

%k=pickic;
k=2;
if k==2
   w1=rands(s1,r);
   b1=rands(s1,1);
   w2=rands(s2,s1);
   b2=rands(s2,1);
   w3=rands(s3,s2);
   b3=rands(s3,1);
end
echo on
clc
%確定顯示頻率、最大步數(shù)、學(xué)習(xí)率———————————————————————————————
df=10;
me=8000;
eg=0.001;
lr=0.4;                  %學(xué)習(xí)速率通常初始值在0.3～0.6之間。
tp=[df me eg lr];
%開始訓(xùn)練網(wǎng)絡(luò)————————————————————————————————————————
[w1,b1,w2,b2,w3,b3,ep,tr]=trainlm(w1,b1,'tansig',w2,b2,'tansig',w3,b3,'purelin',p,t,tp);
%pause

%實時繪制誤差曲線——————————————————————————————————————
clc
plottr(tr,eg);          
%pause

%根據(jù)已經(jīng)訓(xùn)練好的網(wǎng)絡(luò)來仿真輸出———————————————————————————————
echo off
a=simuff(p,w1,b1,'tansig',w2,b2,'tansig',w3,b3,'purelin');
%比較目標(biāo)輸出和仿真輸出圖形————————————————
acs=simuff(cs,w1,b1,'tansig',w2,b2,'tansig',w3,b3,'purelin');
%仿真出測試集的輸出
   figure;
   plot(a);
   hold on
   plot(t,'r');
   title('貨運(yùn)量');
   %比較目標(biāo)輸出和仿真輸出間相對誤差—————————————
   %將歸一化后的目標(biāo)向量和仿真向量還原成原值，進(jìn)行比較————————————————————
   %使用式p2=0.7*(primary(:,2)-minium(2))./(maxium(2)-minium(2))+0.2;反計算出原始數(shù)據(jù)！——   
   a=(a-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);
   t=(t-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);
   rate=((a-t)./t).*100;
   acs=(acs-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);
   tcs=(tcs-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);
   csrate=((acs-tcs)./tcs).*100;
   
   
   %預(yù)測2003年和2004年的貨運(yùn)量——————————————
        a2003=simuff(p2003,w1,b1,'tansig',w2,b2,'tansig',w3,b3,'purelin');
        
                           
                      p2004=[p2(22);
                             p2(23);
                             p2(24);
                             p2(25);
                             a2003];
        a2004=simuff(p2004,w1,b1,'tansig',w2,b2,'tansig',w3,b3,'purelin');
         %換算回來
        a2003=(a2003-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);
        a2004=(a2004-0.2).*((maxium(2)-minium(2))/0.7)+minium(2);


   
   
   
   
   
   
   
   
   
   
   
   %將訓(xùn)練和測試集的結(jié)果按循序排列，并與目標(biāo)比較，若測試集效果好，則網(wǎng)絡(luò)理想——————————   
   for i=1:10
      predict(2*i-1)=acs(i);      
      predict(2*i)=a(i);  
      target(2*i-1)=tcs(i);      
      target(2*i)=t(i); 
      relateerr(2*i-1)=csrate(i);      
      relateerr(2*i)=rate(i);
   end
   
   state=['                             這是貨運(yùn)量的神經(jīng)網(wǎng)絡(luò)模型']
   explain=['               預(yù)測值                          目標(biāo)值               相對誤差']
   compare=[predict' target' relateerr']
   
   figure;
   plot(target,'r-*');
   hold on
   
   l=length(predict);
   predict(l+1)=a2003;
   predict(l+2)=a2004;
   plot(predict,'b-+');
   a2003
   a2004