%問題：用BP網絡預測某省用水總量
%現有某省1978～2002年的客運量數據，即時間序列有Ｎ=25個歷史觀察值,
%采用長度ｎ=5為一期,n0=2為滾動間隔。即用前5年的數據預測下一年的數據,
%每隔2年開始一期長度為5年的觀察數據。要求預測出2003年、2004年的客運量。
%分析：20期數據樣本分成兩組,每組10個，分別做為訓練集（1978年開始）和
%測試集（1979年開始）。設計BP網絡有5個輸入，1個輸出。采用雙隱含層，
%隱含層神經元數目采用試差法確定,經試驗,第1、2隱含層神經元數目分別為8、5，
%即BP網絡設計為5-8-5-1。
%將表1中的輸入、輸出數據比例化到區間[0.2～0.9]內,并將初始權值隨機化,取η=0 6,
%用1978年開始的10個訓練集訓練網絡。經1979年開始的10個測試集測試。 
%1978～2002年的客運量數據
%年份	客運量	
%1978 	12815	
%1979	15543	
%1980	19326	
%1981	22864	
%1982	26150	
%1983	28468	
%1984	30882	
%1985	39375	
%1986	45759	
%1987	49589	
%1988	52560	
%1989	48726	
%1990	51083	
%1991	56495	
%1992	62767	
%1993	83606	
%1994	92090	
%1995	101370	
%1996	107317	
%1997	108654	
%1998	111847	
%1999	112872	
%2000	116997	
%2001	126007	
%2002	128980	

close all
clc;
clear;
%按年份排列的原始客運量數據
p0=[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];

 %繪制客運量隨年份的變換曲線
year=1978:2002;
plot(year,p0,'b+')
hold on
plot(year,p0,'b-.')
pause

 %將原始數據歸一化到[0.2  0.9]
 %規一化公式：規一化后數據＝z1+(z2-z1)*(x-xmin)/(xmax-xmin)
 %[z1 z2]為規格化后的數據范圍，如本例z1=0.2，z2=0.9；x為待規格化數據；
 %xmin、xmax分別為數據序列的最小、最大值。
  p1=0.2+(0.9-0.2).*(p0-min(p0))./(max(p0)-min(p0));

 %繪制規一化后客運量隨年份的變換曲線
 %plot(year,p1)
%pause

%按每5年為一期組織樣本，作為訓練集,p為輸入，t為目標。從第1978年客運量開始
j=0;
for i=1:10;
       p(:,i)=p1(i+j:i+j+4);  %取第i個訓練樣本，給p
       t(:,i)=p1(i+j+5);
       j=j+1;
end
p
t
pause

%按每5年為一期組織樣本，作為測試集,testp為輸入，testt為目標。從1979年客運量開始
j=0;
for i=1:10;
       testp(:,i)=p1(i+j+1:i+j+5);
       testt(:,i)=p1(i+j+6);
       j=j+1;
end

%建立BP網絡，兩層隱含層，隱層神經元數為8、 5，輸出為1個單元，訓練函數為trainlm
net = newff(minmax(p),[8 5 1],{'tansig' 'tansig' 'purelin'},'trainlm');
%可以設計為一個隱含層的BP網絡，隱層神經元數為14；因為數據規一化處理過，
%所以輸出層傳遞函數也可以用'tansig'
%net = newff(minmax(p),[14  1],{'tansig'  'tansig'},'trainlm');
%可以用'trainlm'，'traingd',' traingdm',' traingda',' traingdx'分別訓練網絡，看不同訓練函數的效果.

%設置訓練參數
net.trainParam.lr=0.6;    %學習率
net.trainParam.epochs = 10000; %最大訓練次數
net.trainParam.goal = 0.0000001;  %目標誤差

%訓練網絡
net = train(net,p,t);

%仿真
y = sim(net,testp);         %將1979起的10組測試數據輸入網絡進行測試
E=testt-y;                       %計算測試集網絡輸出和目標的誤差
mse=MSE(E)                 %計算均方誤差

%預測2003年客運情況
p03=p1(21:25)'              %取預測數據，1998－2002共5年，轉置成列矢量
y03=sim(net,p03)         %用sim仿真，計算2003年客運量（規一化值） 

y2003=min(p0)+(y03-0.2)*(max(p0)-min(p0))/0.7;   %求出客運量實際值

p2=[p0 y2003];                          %將2003年數據加入原始數據p0中
year=[year 2003];                      %將2003年加入年份中
figure
plot(year,p2,'g-.')
hold on
plot(year,p2,'g+')                     %繪制包括2003年預測數據的客運量曲線
hold on
plot(2003,y2003,'ro')                %用紅色0繪制2003年預測的數據點。
pause
disp('2003年客運量預測值為：'),y2003
disp('程序運行結束');