第十五章：预测方法

15.1 微分方程模型

描述实际对象的某些特性随时间（或空间）而演变的过程、分析它的变化规律、预测它的未来性态、研究它的控制手段。

dxy=@(t,x)[-0.0544*x(2)+54000*(t>=0 & t<1)+6000*(t>=2 & t<3)+13000*(t>=5 & t<6)
    -0.0106*x(1)];  
[t,xy]=ode45(dxy,[0:36],[0,21500])
subplot(211), plot(t,xy(:,1),'r*',t,xy(:,2),'gD')
xlabel('时间t'),  ylabel('人数'), legend('美军','日军')
subplot(212),  plot(xy(:,1),xy(:,2))  
xlabel('美军人数x'),  ylabel('日军人数y')  

15.2 灰色预测模型

使用的不是原始数据序列，而是生成的数据序列。只适用于中短期的预测，指数增长的预测。

15.2.1 GM(1, 1)预测模型

一阶微分方程，且只含一个变量

1. 数据的检验与处理
2. 建立模型
3. 检验预测值
   （1）残差检验
   （2）级比偏差值检验
4. 预测预报

clc,clear
x0=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]'; %注意这里为列向量
n=length(x0);
lamda=x0(1:n-1)./x0(2:n)  %计算级比
range=minmax(lamda')  %计算级比的范围
x1=cumsum(x0)  %累加运算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
u=B\Y  %拟合参数u(1)=a,u(2)=b
syms x(t)
x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1)); %求微分方程的符号解
xt=vpa(x,6) %以小数格式显示微分方程的解
yuce1=subs(x,t,[0:n-1]); %求已知数据的预测值
yuce1=double(yuce1); %符号数转换成数值类型，否则无法做差分运算
yuce=[x0(1),diff(yuce1)]  %差分运算，还原数据
epsilon=x0'-yuce    %计算残差
delta=abs(epsilon./x0')  %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'  %计算级比偏差值u(1)=a

15.2.2 GM(2 ,1)、DGM和Verhulst模型

非单调的摆动发展序列或有饱和的S形序列
（1）GM(2, 1)模型

clc,clear
x0=[41,49,61,78,96,104]; %原始序列
n=length(x0); 
x1=cumsum(x0)  %计算1次累加序列
次a_x0=diff(x0)' %计算1次累减序列
z=0.5*(x1(2:end)+x1(1:end-1))'; %计算均值生成序列
B=[-x0(2:end)',-z,ones(n-1,1)]; 
u=B\a_x0   %最小二乘法拟合参数
syms x(t)
x=dsolve(diff(x,2)+u(1)*diff(x)+u(2)*x==u(3),x(0)==x1(1),x(5)==x1(6)); 
xt=vpa(x,6) 
yuce=subs(x,t,0:n-1); %求已知数据点1次累加序列的预测值
yuce=double(yuce) 
x0_hat=[yuce(1),diff(yuce)]; %求已知数据点的预测值
x0_hat=round(x0_hat) %四舍五入取整数
epsilon=x0-x0_hat    %求残差
delta=abs(epsilon./x0)  %求相对误差

（2）DGM(2, 1)模型

clc,clear
x0=[2.874,3.278,3.39,3.679,3.77,3.8];
n=length(x0);
a_x0=diff(x0)';  %一阶向前差分
B=[-x0(2:end)',ones(n-1,1)]; 
u=B\a_x0  
syms x(t)
d2x=diff(x,2); dx=diff(x); 
x=dsolve(d2x+u(1)*dx==u(2),x(0)==x0(1),dx(0)==x0(1))
xt=vpa(x,6) 
yuce=subs(x,'t',0:n-1); 
yuce=double(yuce) 
x0_hat=[yuce(1),diff(yuce)] 
epsilon=x0-x0_hat 
delta=abs(epsilon./x0)  

（3）灰色Verhulst预测模型
描述具有饱和状态的过程，即S形过程

clc,clear
x0=[4.93   2.33   3.87   4.35   6.63   7.15   5.37   6.39   7.81   8.35];
x1=cumsum(x0);  
n=length(x0);
z=0.5*(x1(2:n)+x1(1:n-1));   
B=[-z',z'.^2];
Y=x0(2:end)';
u=B\Y     
syms x(t)
x=dsolve(diff(x)+u(1)*x==u(2)*x^2,x(0)==x0(1));  
xt=vpa(x,6) 
yuce=subs(x,'t',[0:n-1]); 
yuce=double(yuce) 
x0_hat=[yuce(1),diff(yuce)] 
epsilon=x0-x0_hat    
delta=abs(epsilon./x0)  
xlswrite('book4.xls',[x0',x0_hat',epsilon',delta'])

15.3 差分方程

15.3.1 商品销量预测

y0=[11 16 25 12 12 18 26 14 13 20 27 15 15 24 30 15 16 25 32 17]';
y=y0(9:20);
x=[y0(5:16),y0(1:12),ones(12,1)];
z=x\y
for t=21:25
    y0(t)=z(1)*y0(t-4)+z(2)*y0(t-8)+z(3);
end
yhat=y0(21:25)   

15.3.2 养老保险

clc, clear
M=600; N=420; p=200; q=2282;
eq=@(x) x^M-(1+q/p)*x^(M-N)+q/p;
x=fzero(eq,[1.0001,1.5]) 

15.4 马尔科夫预测

15.4.1 马尔科夫链的定义

某一系统在已知现在情况的条件下，系统未来时刻的情况只与现在有关，而与过去的历史无直接关系。

15.4.2 转移概率矩阵及柯尔莫哥洛夫定理

clc,clear
format rat  %数据格式是有理分数
fid=fopen('msdata.txt','r');
a=[];
while (~feof(fid))
    a=[a fgetl(fid)];  %把所有字符串连接成一个大字符串行向量
end
for i=0:1
    for j=0:1
        s=[int2str(i),int2str(j)]; 
        f(i+1,j+1)=length(findstr(s,a));  
    end
end
fs=sum(f,2);  %求f矩阵的行和
f=f./repmat(fs,1,size(f,2))  %求状态转移频率

15.4.3 转移概率的渐进性质——极限概率分布

format rat  
p=[0.8 0.1 0.1;0.5 0.1 0.4;0.5 0.3 0.2];
a=[p'-eye(3);ones(1,3)];  %构造方程组ax=b的系数矩阵
b=[zeros(3,1);1];  %构造方程组ax=b的常数项列
p_limit=a\b  
format %恢复到短小数的显示格式

15.5 时间序列

15.5.1 平稳性Daniel检验

对于显著水平$\alpha$，由时间序列$a_i$计算(t, $R_t$)的Spearman秩相关系数q
若|T|>$t_{\alpha /2}$(n-2)，则拒绝$H_0$，认为序列非平稳。且当q>0时，认为序列有上升趋势，q<0时，认为序列有下降趋势。又当|T|>$t_{\alpha /2}$(n-2)时，接受$H_0$，认为序列平稳。

15.5.2 税收收入AR预测模型

clc, clear
a=[15.2	 15.9	18.7	22.4	26.9	28.3	30.5
   33.8	 40.4	50.7	58	    66.7	81.2	83.4];
a=a'; a=a(:); a=a'; 
Rt=tiedrank(a) 
n=length(a); t=1:n; 
Qs=1-6/(n*(n^2-1))*sum((t-Rt).^2)   %计算q的值
T=Qs*sqrt(n-2)/sqrt(1-Qs^2)   %计算T统计量的值
t_0=tinv(0.975,n-2)     %计算上alpha/2分位数
b=diff(a)   
m=ar(b,2,'ls')  
bhat=predict(m,b')  
bhat(end+1)=forecast(m,b',1) 
ahat=[a(1), a+bhat']  
delta=abs((ahat(1:end-1)-a)./a)  
xlswrite('yu.xls',ahat), xlswrite('yu.xls',delta,'Sheet1','A3') 

15.6 插值与拟合

15.6.1 导弹运动轨迹问题

clc,clear
syms x y
syms z positive  %定义符号变量z为正
format long g  
a=load('daodan.txt'); 
d=a(:,[2:end]);  
n=size(a,1); sol=[]; 
for i=1:n
eq1=x^2+y^2+z^2-d(i,1)^2;  
eq2=x^2+(y-4500)^2+z^2-d(i,2)^2;
eq3=(x+2000)^2+(y-1500)^2+z^2-d(i,3)^2; 
[xx,yy,zz]=solve(eq1,eq2,eq3);  
sol=[sol;double([xx,yy,zz])];  
end
sol 
pp1=csape(a(:,1),sol(:,1)) 
xishu1=pp1.coefs(end,:)  
pp2=csape(a(:,1),sol(:,2))  
xishu2=pp2.coefs(end,:)  
pp3=csape(a(:,1),sol(:,3))  
xishu3=pp3.coefs(end,:)  

15.6.2 录像机计数器的用途

clc, clear
a=load('jishu.txt');  
t=a([1,3],:); t=t(:); 
n=a([2,4],:); n=n(:); 
xishu=[n.^2,n]; %构造系数矩阵
ab=xishu\t
n0=4450
that=ab(1)*n0^2+ab(2)*n0

15.7 神经元网络

15.7.1 BP神经网络

1. 初始化网络及学习参数
2. 提供训练模式，训练网络，直到满足学习要求
3. 前向传播过程
4. 反向传播过程

15.7.2 RBF神经网络

1. 确定基函数的中心C
2. 确定基函数的宽度d
3. 确定从隐含层到输出层的连接权值w

15.7.3 年径流预报实例

clc, clear
a=load('jingliu.txt'); 
a=a';
P=a([1:4],[1:end-1]); [PN,PS1]=mapminmax(P); 
T=a(5,[1:end-1]);  [TN,PS2]=mapminmax(T); 
net1=newrb(PN,TN)  %训练RBF网络
x=a([1:4],end); xn=mapminmax('apply',x,PS1); 
yn1=sim(net1,xn); y1=mapminmax('reverse',yn1,PS2) 
delta1=abs(a(5,20)-y1)/a(5,20)  %求预测的相对误差
net2=feedforwardnet(4);  %初始化BP网络，隐含层的神经元取为4个
net2 = train(net2,PN,TN); 
yn2= net2(xn); y2=mapminmax('reverse',yn2,PS2)