灰色预测

1.GM（1,1）

变量一阶，微分方程一阶

1.数列预测

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)

for i=2:n

z(i)=0.5*(x1(i)+x1(i-1));

end

B=[-z(2:n)',ones(n-1,1)];

Y=x0(2:n)';

u=B\Y

x=dsolve('Dx+a*x=b','x(0)=x0');

x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});

yuce1=double(subs(x,'t',[0:n-1]));

digits(6),y=vpa(x) %为提高预测精度，先计算预测值，再显示微分方程的解

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 %计算级比偏差值

说明：delta,rho<0.2，预测可接受

2.灾害预测

挑出异常点组成新的数列，大于上限或者小于下限

对时间求cumsum，来预测下一个发生灾害的时间点

clc,clear

a=[390.6,412,320,559.2,380.8,542.4,553,310,561,300,632,540,406.2,313.8,576,587.6,318.5]';

t0=find(a<=320);

t1=cumsum(t0);n=length(t1);

B=[-0.5*(t1(1:end-1)+t1(2:end)),ones(n-1,1)];Y=t0(2:end);

r=B\Y

y=dsolve('Dy+a*y=b','y(0)=y0');

y=subs(y,{'a','b','y0'},{r(1),r(2),t1(1)});

yuce1=subs(y,'t',[0:n+1])

digits(6),y=vpa(y) %为提高预测精度，先计算预测值，再显示微分方程的解

yuce=diff(yuce1);

yuce=[t0(1),yuce]