function main()
clc
clear all
close all
warning off
samnum=20;
testsamnum=20;
forcastsamnum=2;
hiddenunitnum=8;
indim=3;
outdim=2;
A=[20.6 22.4 25.4 27.1 29.4 30.1 30.9 34.1 36.4 38.1 39.1 40.0 41.9 44.6 47.3 52.9 55.7 56.8 59.1 60.6];
B=[0.6 0.75 0.85 0.90 1.05 1.35 1.45 1.60 1.70 1.85 2.15 2.20 2.25 2.35 2.50 2.60 2.70 2.85 2.95 3.10];
C=[0.09 0.11 0.11 0.14 0.20 0.23 0.23 0.32 0.32 0.34 0.36 0.36 0.38 0.49 0.56 0.59 0.59 0.67 0.69 0.79];
D=[5126 6217 7730 9145 10460 11387 12353 15750 18304 19836 21024....
    19490 20433 22598 25107 33442 36836 40548 42927 43462];
E=[1237 1379 1385 1399 1663 1714 1834 4322 8132 8936 11099 11203....
    10524 11115 13320 16762 18673 20724 20803 21804];
p=[A;B;C];
t=[D;E];
[samin,minp,maxp,tn,mint,maxt]=premnmx(p,t);
rand('state',sum(100*clock))
noisevar=0.01;
noise=noisevar*randn(2,samnum)
samout=tn+noise;
testsamin=samin;
testsamout=samout;
maxepochs=5000;
lr=0.035;
e0=0.65*10^(-3);
w1=0.5*rand(hiddenunitnum,indim)-0.1;
b1=0.5*rand(hiddenunitnum,1)-0.1;
w2=0.5*rand(outdim,hiddenunitnum)-0.1
b2=0.5*rand(outdim,1)-0.1;
errhistory=[];
for i=1:maxepochs
    hiddenout=logsig(w1*samin+repmat(b1,1,samnum));
    networkout=w2*hiddenout+repmat(b2,1,samnum);
    error=samout-networkout;
    sse=sumsqr(error)
    errhistory=[errhistory sse];
    if sse<e0
        break
    end
    delta2=error;
    delta1=w2'*delta2.*hiddenout.*(1-hiddenout);
    dw2=delta2*hiddenout';
    db2=delta2*ones(samnum,1);
    dw1=delta1*samin';
    db1=delta1*ones(samnum,1);
    w2=w2+lr*dw2;
    b2=b2+lr*db2;
    w1=w1+lr*dw1;
    b1=b1+lr*db1;
end
hiddenout=logsig(w1*samin+repmat(b1,1,testsamnum));
networkout=w2*hiddenout+repmat(b2,1,testsamnum);
a=postmnmx(networkout,mint,maxt);
x=1990:2009;
newk=a(1,:);
newh=a(2,:);
figure;
subplot(2,1,1);
plot(x,newk,'r-o',x,D,'b--+');
legend('网络输出客运量','实际客运量');
xlabel('年份');
ylabel('客运量/万人');
title('神经网络客运量学习和测试对比图');
subplot(2,1,2);
plot(x,newh,'r-o',x,E,'b--+');
legend('网络输出货运量','实际货运量');
xlabel('年份');
ylabel('货运量/万吨');
title('神经网络货运量学习和测试对比图');
pnew=[73.4 75.6;4.0 4.1;0.9 1.1];
pnewn=tramnmx(pnew,minp,maxp);
hiddenout=logsig(w1*pnewn+repmat(b1,1,forcastsamnum))
anewn=w2*hiddenout+repmat(b2,1,forcastsamnum)
anew=postmnmx(anewn,mint,maxt)