灰色预测

                               **灰色预测算法简介**

1、灰色预测模型（数模必掌握）
灰色预测模型使用范围：
①数据样本点个数少，6-15个
②数据呈现指数或曲线的形式
③只适合做中短期预测，不适合长期预测。
灰色预测原理比较简单，详细的可以参考司守奎《数学建模算法与应用》。
需要注意的几点是：
（1）灰色预测的使用范围
（2）灰色预测中的“级比”如果级比不在范围要对数据进行处理。
（3）司老师书中的代码，并没有运行出后面的运行结果，如果想运行出预测的结果，看下面的说明。
（4）在使用灰色预测的时候要考虑残差等（见代码的最后三行）
（5）代码直接复制粘贴文本文档的文件就可以了。
（6）文本文档是给出了两种代码，不要复制错了，第一个是司老师书中的。第二个是实践，某城市 1986～1992 年道路交通噪声平均声级数据见表1。

                         表1 城市交通噪声数据/dB(A)

序号 年份 Leq 序号 年份 Leq
1 1986 71.1 5 1990 71.4
2 1987 72.4 6 1991 72.0
3 1988 72.4 7 1992 71.6
4 1989 72.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); %累加运算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
u=B\Y
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]);%为提高预测精度，先计算预测值，再显示微分方程的解。%%%这里，把[0:n-1]修改就可以了，如果预测后5年的，就改成n+4。
yuce1=double(yuce1);%符号数转换成数值类型，否则无法做差分运算
yuce=[x0(1),diff(yuce1)] %差分运算，还原数据
yuce2=subs(x,t,[0:n+4]);%为提高预测精度，先计算预测值，再显示微分方程的解。
Yuce2=double(yuce2);%符号数转换成数值类型，否则无法做差分运算
yucexin=[x0(1),diff(yuce2)] %差分运算，还原数
【直接把绿色总n-1修改，运算会提示错误，但是不影响，如果想消除错误提示，可以添加黄色区域的，计算出两个结果，一个结果是不预测的一个结果是预测的。】
epsilon=x0'-yuce %计算残差
delta=abs(epsilon./x0') %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda' %计算级比偏差值
灰色预测模型使用简单，基本上可以实现黑箱操作。该模型在实际数学建模竞赛中，只需要把矩阵 换成自己的矩阵即可。改变“xt=vpa(x,6)”中的“6”即可改变显示精度，例如：xt=vpa(x,7)（一般六位显示精度满足数学建模数据精度的需要）。
------下面是代码：
-----------------------代码一----------------------
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
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';%计算级比偏差值


----------------------代码二------------------------------
clc,clear;  
syms a b;  
c=[a b]';  
A=[89677,99215,109655,120333,135823,159878,182321,209407,246619,300670];  
B=cumsum(A);  %原始数据累加  
n=length(A);  
for i=1:(n-1)  
    C(i)=(B(i)+B(i+1))/2; %生成累加矩阵  
end  
%计算待定参数的值  
D=A;D(1)=[];  
D=D';  
E=[-C;ones(1,n-1)];  
c=inv(E*E')*E*D;  
c=c';  
a=c(1);b=c(2);  
%预测后续数据  
F=[];F(1)=A(1);  
for i=2:(n+10)  %只推测后10个数据，可以从此修改  
    F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;  
end  
G=[];G(1)=A(1);  
for i=2:(n+10)  %只推测后10个数据，可以从此修改  
    G(i)=F(i)-F(i-1);  %得到预测出来的数据  
end  
t1=1999:2008;  
t2=1999:2018;  %多10组数据  
G  
h=plot(t1,A,'o',t2,G,'-'); %原始数据与预测数据的比较