灰色预测MATLAB程序

合集下载

灰色系统G(1,1)预测步骤【模板带代码】

=3499.075e -0.1062t

3641.075
编写程序
u=alpha(2)/alpha(1) v=X0(1)-u v=3499.075 u=—3641.075
(5)进行参差检验
1）根据预测公式，计算
v=3499.075 u=—3641.075
Xˆ
1
k
1

X
0
1
for n=0:10
X2(n+1)=v*exp(-alpha(1)*n)+u
end
X2
2.0690
2）累减生成序列
Xˆ X3 =1.0e+003 * (0) 0.1420 0.4079 0.4536
0.5044
0.7713 0.8577 0.9537 1.0605
源程序：X3(1)=X2(1)
for m=1:10
kesi =
4.4388 339.0664 176.2445 203.6132
0 0.1998 1.2682 0.2130 0.8524 0.1330
0.0089
0.3767
0.2203
0.4155
{0%，19.98%，126.82%，0.89%，37.67% ，22.03% ，41.55% ，21.30%，85.24%，13.30%}
e=
179.4592 111.5134 74.1747 175.0204 159.6072 29.2461 215.2168 33.1910
3.2147 24.1540
源程序：S0=0.6745*X0std e=abs(daita0-daita0mean) 对所有的 e 都小于 S0 ，故小参差概率 P(k S0) 1 0.95

灰色预测MATLAB程序

灰色预测作用：求累加数列、求a b的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.1 72.4 72.4 72.1 71.4 72.0 71.6];format long; %设置计算精度if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x';endn=length(x); %取输入数据的样本量z=0;for i=1:n %计算累加值，并将值赋予矩阵bez=z+x(i,:);be(i,:)=z;endfor i=2:n %对原始数列平行移位y(i-1,:)=x(i,:);endfor i=1:n-1 %计算数据矩阵B的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:));endfor j=1:n-1 %计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1 %构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,: );%显示输出预测值的累加数列endvar(1,:)=ago(1,:) %显示输出预测值for i=1:n %如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列alpha %显示输出参数数列var %显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程clc,clearx0=[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)=bx=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)}) %代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]); %求已知数据的预测值y=vpa(x,6) %其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据。

灰色预测MATLAB程序

作用：求累加数列、求a b的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[ ];format long; %设置计算精度if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x';endn=length(x); %取输入数据的样本量z=0;for i=1:n %计算累加值，并将值赋予矩阵bez=z+x(i,:);be(i,:)=z;endfor i=2:n %对原始数列平行移位y(i-1,:)=x(i,:);endfor i=1:n-1 %计算数据矩阵B的第一列数据c(i,:)=*(be(i,:)+be(i+1,:));endfor j=1:n-1 %计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1 %构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,: );%显示输出预测值的累加数列endvar(1,:)=ago(1,:) %显示输出预测值for i=1:n %如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列alpha %显示输出参数数列var %显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程clc,clearx0=[ ]'; %注意这里为列向量n=length(x0);lamda=x0(1:n-1)./x0(2:n) %计算级比range=minmax(lamda') %计算级比的范围x1=cumsum(x0) %累加运算B=[*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];Y=x0(2:n);u=B\Y %拟合参数u(1)=a,u(2)=bx=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)}) %代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]); %求已知数据的预测值y=vpa(x,6) %其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据。

灰色预测模型matlab程序精确版

%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型function gmcal=gm1(x)if nargin==0x=[1019,1088,1324,1408,1601]endformat long gsizex=length(x);%求数组长度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，用于计算Byn1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k,yn1sizez1=length(z1);%size(yn1);z2 = z1';z3 = ones(1,sizez1)';YN = yn1'; %转置%YNB=[z2 z3];au0=inv(B'*B)*B'*YN;au = au0';%B,au0,auafor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,ua%输出预测的 a u 和 u/a的值constant1 = x(1)-ua;afor1 = -afor;x1t1 = 'x1(t+1)';estr = 'exp';tstr = 't';leftbra = '(';rightbra = ')';%constant1,afor1,x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra)%输出时间响应方程%******************************************************%二次拟合k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(-(k2-1)*afor);endend%ze1sizeze1=length(ze1);z4 = ones(1,sizeze1)';G=[ze1' z4];X1 = x1';au20=inv(G'*G)*G'*X1;au2 = au20';%z4,X1,G,au20Aval = au2(1);Bval = au2(2);%Aval,Bval%输出预测的 A,B的值strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra)%输出时间响应方程nfinal = sizex-1 + 1;（其中+1可改为+5等其他数字，即可预测更多的数字）%决定预测的步骤数5 这个步骤可以通过函数传入%nfinal = sizexd2 - 1 + 1;%预测的步骤数 1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fcast%一次拟合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendendx31fcast%一次拟合预测值for k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendendx41fcast,x%二次拟合预测值%***精度检验p C************////////////////////////////////// k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexelseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses1total = s1total + (x(k5) - xavg)^2;endends1suqare = s1total ./ sizex;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 残差数列x的方差 s1sqrt 为x方差的平方根S1 k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses2total = s2total + (err1(k5) - err1avg)^2; endends2suqare = s2total ./ sizex;%s2suqare 残差数列err1的方差S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C检验值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )elseendendpval = pnum ./ sizex;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围预测步长和数据长度可调整程序参数进行改进运行结果x =1019 1088 1324 1408 1601ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929)ans =x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685)x31fcast =Columns 1 through 31019 1122.89347857097 1264.43142178303 Columns 4 through 61423.80987235488 1603.27758207442 1805.36675232556 x41fcast =Columns 1 through 31019 1118.05685435129 1269.65470492098 Columns 4 through 61429.69153740195 1609.90061644041 1812.82460377782 x =1019 1088 1324 1408 1601Cval =0.139501578334155 pval =1。

分数灰色预测matlab代码详解

分数灰色预测matlab代码详解
分数灰色预测是一种基于灰色系统理论的非线性预测方法，通过对待预测序列的数据进行分形分析，建立分数阶微分方程模型，从而进行预测。

下面我们将详细介绍该方法的matlab代码实现过程。

1. 数据准备
首先，需要准备待预测的时间序列数据，在matlab中可以通过读取文件或手动输入的方式获取数据。

在本文中，我们将使用matlab 自带的load函数读取一个名为data.txt的文本文件中的时间序列数据。

2. 数据预处理
在进行预测之前，需要对数据进行预处理。

这包括去除噪声、平滑处理、归一化等。

在本文中，我们将采用matlab中自带的smooth 函数进行平滑处理，并使用归一化方法将数据缩放到0至1之间。

3. 模型建立
接下来，需要建立分数灰色预测模型。

在matlab中，可以使用greyest函数进行模型参数估计。

在本文中，我们将使用分数阶微分方程模型，因此需要先通过fracdiff函数估计分数阶微分系数。

4. 模型预测
有了模型之后，就可以进行预测了。

在matlab中，可以使用sim 函数进行模型仿真。

在本文中，我们将使用该函数对模型进行预测，并将预测结果可视化。

5. 结果分析
最后，需要对预测结果进行分析。

可以通过计算误差指标、绘制误差曲线等方式进行分析。

在本文中，我们将计算均方误差和平均绝对误差，并绘制预测结果和实际结果的对比图。

综上所述，以上就是分数灰色预测的matlab代码详解。

通过对上述步骤的实现，可以得到较为准确的预测结果，并帮助我们更好地了解该预测方法的原理和应用。

灰色预测MATLAB程序

灰色预测心设尹曲⑴#为原始数列，其1次累＜加生成数列为炉=(孝①宀2\S,其中©=2^°:⑺卫=12…止i-1尋定文沙的灰导数为d(Jt)=玄㈣(Jt)=尤⑴的-工⑴(*-1).令尹为数列壬⑴的邻值生成数列，即尹)(町=加小(町十(1—a)x山(k-1).于是定文GM(1T1)的灰微分方程模型为d(k)+az①(上)=&_即或.严⑹+盘⑴懐)=乩⑴在式(1)中口①的称为灰导数’熬称为发展系数'弧称为白化背景值，b称为灰作用量。

将时刻表庄=23…用代入(O式有j<0)(2)-az⑴(2)=工®⑶—俺叫巧=»于是GMIL)樫型可表示为r=现在问题归结为求巧h在值。

用一元绒性回归，即最小二垂進求它们的估计值住=[]卜护跖护F奕厢上回归分析中求诂计值是用软件计算的，有标淮程博求解，如山訥甜等。

GM(1.1)的白化型对于的(1-1)的获微分方程⑴，如果将解导教矿悶的时報=%…屮观対连续叢里"则工⑴衩为时间i函敕卅®,于是-<'W耐应于导敕重级必％),白化背杲值刃(時对应于导數申⑴。

于是GM(1，1)的换微分方嗨对应于的白微分方程为写®4曲％「)=也⑵GAI(1>1)换色预刪的步叢1-數堀的椅噓弓处理为了保证©M(B1)屋複方达的可行性・需要対已却皴堀锁必要的检峻处Ho 设療皓数攥列为了-计算埶列的级比如果所有的级比都落在可容覆盖区间盂-內・则數摒列X糾可咲建立G*ICL-1)複型且可以避行页色预测。

否则，丙軌据懺适当的叢换处理，如平移銮换：取C使得敕培列严⑹二工蚀盘)+匚用二12…”的级比都落在可啓禎盖内。

(1)残差檢验：计算相对薙差Z 建立GM (L T 1)複型不妬设少弋以m 叫唠霸足上面的要求，以它芮議堀列建立GM(1>1)型蛊(仍(i)+血C1\A)=b ・用回归分祈求得目上的估计值"于是相应的白化模型为气^十小卄工解为工叱)=0)①—勺中1-色-⑶ 应Q于是停到预测值壬⑴(上+1)=0叫1)一勺>加+仝血二12…卫一1=aa伙而相应地得到预«=x co \t +1)=x 0)(t+l)-x a)(i)3i =1,2,-?n-l ?如果对所有的^<0.1・则认为达到鞭嵩的要求：否则，若耐所有的|^)1<0^，则认対达到一般要求©(2)级比偏差値桧验：计算能)=1-呂学©如果对所有的|，则认为达列较高的要求孑吾则若对斫有的，则认为达到一般要求O灰色预测计算实例^…;=:=-■■■■昏例北方某城市1986—1992年道路交通噪声平均声级数据见表6序号年吶寺表拆市近年来交通噪声数据[眶(应)]二諾;二319S872.4第—爭:级比检验建立丢通噪屛均声级数锯时间序列如下：4198972.1j 1990?1.4 619?17201199771.6艸=(•严①卫购(2)厂卅⑺) =(711,72.4.71.4,72.1.71.4,7UQ.71.6)些(1)求级比k(k)忠防护住T)2=(几⑵山⑶.…也⑺)g=(0.982JJ.0042J.0098-0.9917J.0056)(2)级比判断由于所有的X.(10e[0.982J.009S],k=2,3.6故可以用双0)作满意的GM(1,1)建模’第二步：GM(1,1)建模(1)对原始数据X®作一次累加，即卞⑴=(71.L143.5215.9.288359.4.431.4,503)(2)构造数据矩阵B及数据向量Y-2)—H 弋3/>1⑶讦算1T心求解得F'⑴=(工倒〔1〉_-)e 弋Q f+-1*^+1)=0＜l,U)--)£-t +-=-3092^--^+31000⑶求生咸数列值歸型齊看:n令“is 那血由上面的碉醯数可甲得，其中取菱由龙⑴(i}=恥壮曲5加得丁I —"炉閃=进悶-进德-尊(71儿72.4.72.2:72.1:71.9:71.7,71.6)^}=(s"a ＞亍⑴⑵，…，网⑺A＜第三步;模型检验•＞模型的各种检验指标值的计算结果见表工 •t*表7GM(1检验表＜序号年俯原始值模型值残差相对误差级比偏差•>1 19S6 71.1 71.1<219S7 72.4 72.4 -0.0057 0.01%0.0023 <3 19S S 72.4 72.2 0.163S 0.23%0.0203 •>4 19S9 72.1 72.1 0.0329 0.05%-O.(K H8 •>5199071.4 71.9 -0-49S4 0.7%-0.0074 <61991 72.0 71.7 0.21599 037%0.0107<71992 71.6 71.6 0.037S0.05%-0.0032于是得到目=山的餡,立欖型7-B)'1B TY=(dt0.0023 72.6573dt+0.002ix (1>=72.657^心经验证・该模型的精度较高.可进行预测和预报计算的Matlab 程序如下:仃坝测和预报n=length(x); z=0;%取输入数据的样本量for i=1:nz=z+x(i,:)be(i,:)=z; %计算累加值，并将值赋予矩阵beend for i=2:n %对y(i-1,:)=x(i,:)%对原始数列平行移位 endfor i=1：n-1%计算数据矩阵B 的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:)); clCjdearxO=[71H 72.472A 72J71477m c n.lengthtxO)；*'b%注意这里为列帖lamda =xD(l ：n-1),A0(2：n)%计算级比range =minmaxflamda f )%计算级比的范阖 X1=cumsum(xO)；%累加运算B=['0,5*(xl(l ；n ^l)+xl(2:n))t ones(n -1,1)]TY 二甸(2：町；口=B\Y%拟合参数u(l>=a .u(2)=bx=dsolve (+a 'x =b\f x(0)-xO^J ；%求徴分方程的特号解x =subs(xJ*a\,b r /xO ，Mu(l)P u(2)t xO(l)|)i%代入荷计痹擞值和初蜡值yucel =subs %求巳知数擁的扳测位y-vpa(x,6)奄其中的石表示显不白位数字yuce=[x0(l)T diff(yucel)]%羔分运算，还原数据 epsiIon=-yuce%计算战羞作用：求累加数列、求ab 的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.172.472.472.171.472.071.6]; format long ;%设置计算精度if length(x(:,1))==1%对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x endM.I-JTVorhlllst 模型endfor j=1:n-1%计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1%构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y;%计算参数矩阵即ab的值for i=1:n+1%计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha( 2,:)/alpha(1,:);%显示输出预测值的累加数列endvar(1,:)=ago(1,: )for i=1:n%显示输出预测值%如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:);%估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:);%计算残差endc=std(error)/std(x);%调用统计工具箱的标准差函数计算后验差的比值cago alpha var%显示输出预测值的累加数列%显示输出参数数列%显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求ab的值、求预测方程clc,clearx0=[71.172.472.472.171.472.071.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)=bx=dsolve('Dx+a*x=b','x(0)=x0');%求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)})%代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]);%求已知数据的预测值y=vpa(x,6)%其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)]%差分运算，还原数据。

灰色理论预测模型及GM(1,1)matlab程序

灰色理论预测模型及GM(1,1)matlab程序灰色预测方法简介灰色预测是一种对含有不确定因素的系统进行预测的方法。

灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

其用等时距观测到的反应预测对象特征的一系列数量值构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

通过对原始数据的整理寻找数的规律，分为三类：a、累加生成：通过数列间各时刻数据的依个累加得到新的数据与数列。

累加前数列为原始数列，累加后为生成数列。

b、累减生成：前后两个数据之差，累加生成的逆运算。

累减生成可将累加生成还原成非生成数列。

c、映射生成：累加、累减以外的生成方式。

建模步骤a、建模机理b、把原始数据加工成生成数；c、对残差（模型计算值与实际值之差）修订后，建立差分微分方程模型；d、基于关联度收敛的分析；e、gm模型所得数据须经过逆生成还原后才能用。

f、采用“五步建模（系统定性分析、因素分析、初步量化、动态量化、优化）”法，建立一种差分微分方程模型gm(1,1）预测模型。

GM(1,1)程序：% 本程序主要用来计算根据灰色理论建立的模型的预测值。

% 应用的数学模型是GM(1,1)。

% 原始数据的处理方法是一次累加法。

clear;clc;% load ('data.txt');% y=data';y=[3 4 5 4 7 7];n=length(y);yy=ones(n,1);yy(1)=y(1);for i=2:nyy(i)=yy(i-1)+y(i);endB=ones(n-1,2);for i=1:(n-1)B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1;endBT=B';for j=1:n-1YN(j)=y(j+1);endYN=YN';A=inv(BT*B)*BT*YN;a=A(1);u=A(2);t=u/a;t_test=input('请输入需要预测个数：');i=1:t_test+n;yys(i+1)=(y(1)-t).*exp(-a.*i)+t;yys(1)=y(1);for j=n+t_test:-1:2ys(j)=yys(j)-yys(j-1);endx=1:n;xs=2:n+t_test;yn=ys(2:n+t_test);plot(x,y,'^r',xs,yn,'*-b');det=0;for i=2:ndet=det+abs(yn(i)-y(i));enddet=det/(n-1);disp(['百分绝对误差为：',num2str(det),'%']); disp(['预测值为：',num2str(ys(n+1:n+t_test))]);。

基于MATLAB的灰色预测模型

现有给定数据（见表1），每一组编程的同学先熟悉灰色预测)1,1(GM 的算法，根据算法过程（分步完成）编写相应的程序。

要求：
1.程序能够由编程者的需求不同预测出给定数据后任意n 年的值；
2.程序能够求出给定数据与预测数据的残差；
3.程序能够画出初始数据与预测数据的对比图，注意图的格式，黑白；
4.程序能够求出数据矩阵B 及数据向量Y ；
5程序能够得出求预测值的方程即：
c me k x nk +=+)1(^
其中，c n m 、、均为需用程序求解出来的系数；
6.程序中重要的程序语句用‘%’注明语义；
表1 给定数据。

灰色预测模型的matlab运行代码(讲解)

1灰色预测模型 GM(1,1)的matlab 运行代码例由1990—2001年中国蔬菜产量，建立模型预测 2002年中国蔬菜产量，并对预测结果作检验。

分析建模：给定原始时间1990—2001年资料序列X (0)(k),对X (0)(k) 生成1-AGO(累加)序列X (1)(k)及Y n。

见下表K 1 2 3 4 5 6 7 8 9 10 11 12 X ⑼19519 ,19578 ,19637,19695,16602, 25723, 30379, 34473, 38485, 40514, 42400, 48337 X (1)19519, 39097, 58734, 264605,307005,355342Y n 19578 19637 40514 42400 48337其中 X ⑴(k):x (0)(i)；i 1Y n = [X (0)(2),X (0)(3), ,X (0) (12)]TX (0)(k)的 GM(1,1),得到z (1^2) 10.5 X⑴(1) X (。

⑵1 -29308.0 z ⑴(3) 1 0.5 X ⑴(2) X "(3)1 -48915.5 z (1)(4) 10.5 X ⑴(3) X (1^4)1 -68581.5 z ⑴(5) 1 0.5 X(4) X (1^5)1 -86730.0 z (1)(6) 10.5 X ⑴(5) X⑴⑹ 1 -107892.5 z ⑴(7) 1 0.5 X ⑴(6) X ("(7)1 -135943.5 z ⑴(8) 1 0.5 X ⑴(7) X "(8)1 -168369.5 z ⑴(9) 1 0.5 X ⑴(8) X ⑴(9)1-204848.5 z ⑴(10) 1 0.5 X ⑴(9) X ⑴(10) 1 244348.0 z ⑴(11) 10.5 X ⑴(10) X ("(11) 1 -331173.5 z ⑴(12) 10.5 X ⑴(11)X ("(12) 1-5236.21 1 1 1 1 1 1 1 1 1 将B 和Y n 代入辨识算式,有:对上述 ka(B TB) 1B T?Y n b得灰色GM(1,1)模型为0.1062105 13999.9(1)灰微分方程 X (0) (k)-0.1062105 Z ⑴(k)=13999.9刃⑼亿 1) 15248.968e 0.1062105t(5)残差检验：(6)后验差检验:12[X ( 0)(i) X (0)]2112L , X (0)丄 X (0)(i) n 1 12 i 1 C<0.35好；C<0.5合格；00.6不合格。

灰色预测模型实验以及例题分析实验报告

Columns 19 through 20
1.0691 1.2371
预测数据与实际数据的比较图如下：
2.P3-2（P34）
结果如下：
a=-0.0624
b =156.6162
预测数据与实际数据的比较图如下：
3.预测与会代表人数”（P35）
结果如下：
G=
Columns 1 through 9
0.1810 0.1980 0.1735 0.1520 0.1331 0.1166 0.1021 0.0895 0.0784
实验名称
灰色预测模型实验
实验目的
用MATLAB编程验证程序P3-1（P32）及P3-2（P34），并用其解决“预测与会代表人数”（P35）问题。
学会使用matlab编程做灰色预测实验。
实验内容（算法、程序、步骤和方法）
实验内容（算法、程序、步骤和方法）
实验内容（算பைடு நூலகம்、程序、步骤和方法）
1.P3-1（P32）
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)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据
end
t1=1:4;
t2=1:14;
G
plot(t1,A,'o',t2,G) %原始数据与预测数据的比较

灰色预测及MATLAB实现

的未来预测值。
（3）对累加生成数据做均值生成 B 矩阵与常数项向量Yn ，即
0.5(x(1) (1) x(1) (2))
B

0.5(
x(1)
(2)

x(1)
(3))

，Yn

(x(0) (2), x(0) (3),
, x(0) (n))T

0.5(x(1) (n 1) x(1) (n))
3.2 灰色预测的MATLAB程序
3.2.1 典型程序结构
（1）对原始数据进行累加。
矩阵处理， MATLAB的长
处
（2）构造累加矩阵B与常数向量。
（3）求解灰参数。
（4）将参数代入预测模型进行数据预测。
【例】某公司1999-2008年利润为（元/年）：[89677 99215 109655 120333 135823 159878 182321 209407 246619 300670]，预测该公司未来几年的利润情况。
已知本届会议的回执情况（表1），往几届会议代表回执和与会情况（表2），根据这些数据预测本届与会代表。
表1 回执中对住房的要求
要求男女
合住1 154 78
合住2 104 48
合住3 32 17
独住1 107 59
独住2 68 28
独住3 41 19
表2 以往几届代表的回执参会情况表
届次
第一届
第二届
dx
由于aˆ 是通过最小二乘法求出的近似值，因此 xˆ(1) (t 1)事近似表达
式，与原序列区分，多了一个“帽子”。
（6）对函数表达式 xˆ(1) (t 1)及 xˆ(1) (t)进行离散，将二者作差以便还

灰色预测例题matlab求解(精)

建立GM（1，1）模型对产品销售额预测祁诗阳冯晓凯申静某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1预测模型，并预测2005年的产品销售额。

年份1999 2000 2001 2002 2003 2004销售额2.673.13 3.25 3.36 3.56 3.72(亿元有题目知构造累加生成序列对作紧邻均值生成编程如下：x=[2.67 5.8 9.05 12.41 15.97 19.69];z(1=x(1;for i=2:6z(i=0.5*(x(i+x(i-1;endformat long gz结果如下：z =Columns 1 through 42.67 4.235 7.425 10.73Columns 5 through 614.19 17.83因此于是构造B矩阵和Y矩阵如下：对参数进行最小二乘估计，采用matlab编程完成解答如下：B=[[-4.235 -7.425 -10.73 -14.19 -17.83]',ones(5,1];Y=[3.13 3.25 3.36 3.56 3.72]';format long ga=inv(B'*B*B'*Y结果如下：a =-0.04396098154759662.92561659879905即=-0.044，u=2.96 =-66.55则GM(1,1白化方程为预测模型为：1、关联度检验法：采用matlab编程得到模拟序列for i=1:6X(i=69.22*exp(0.044*(i-1-66.55;endformat long gx(1=X(1;for i=2:6x(i=X(i-X(i-1;endX结果如下：x =Columns 1 through 42.673.11367860537808 3.25373920141375 3.40010005288617 Columns 5 through 63.55304456012134 3.71286887145915因此模拟序列为求模拟序列和原始序列的相关度初始化，即将该序列所有数据分别除以第一个数据。

灰色预测系统基于GM(1-1)的matlab程序

function GM1_1(X0)%format long ;X0=input('请输入实测数据');%实测值[m,n]=size(X0);X1=cumsum(X0); %累加X2=[];for i=1:n-1X2(i,:)=X1(i)+X1(i+1);endB=-0.5.*X2 ;t=ones(n-1,1);B=[B,t] ; % 求B矩阵YN=X0(2:end) ;Pt=YN./X1(1:(length(X0)-1)) %对原始数据序列X0进行准光滑性检验， %序列X0的光滑比P(t)=X0(t)/X1(t-1)A=inv(B.'*B)*B.'*YN.' ;a=A(1)u=A(2)c=u/a ;b=X0(1)-c ;X=[num2str(b),'exp','(',num2str(-a),'k',')',num2str(c)];strcat('X(k+1)=',X)%syms k;for t=1:length(X0)k(1,t)=t-1;endkY_k_1=b*exp(-a*k)+c;for j=1:length(k)-1Y(1,j)=Y_k_1(j+1)-Y_k_1(j);endXY=[Y_k_1(1),Y] %预测值CA=abs(XY-X0) ; %残差数列Theta=CA %残差检验绝对误差序列XD_Theta= CA ./ X0 %残差检验相对误差序列AV=mean(CA); % 残差数列平均值R_k=(min(Theta)+0.5*max(Theta))./(Theta+0.5*max(Theta)) ;% P=0.5R=sum(R_k)/length(R_k) %关联度Temp0=(CA-AV).^2 ;Temp1=sum(Temp0)/length(CA);S2=sqrt(Temp1) ; %绝对误差序列的标准差%----------AV_0=mean(X0); % 原始序列平均值Temp_0=(X0-AV_0).^2 ;Temp_1=sum(Temp_0)/length(CA);S1=sqrt(Temp_1) ; %原始序列的标准差TempC=S2/S1*100; %方差比C=strcat(num2str(TempC),'%') %后验差检验 %方差比%----------SS=0.675*S1 ;Delta=abs(CA-AV) ;TempN=find(Delta<=SS);N1=length(TempN);N2=length(CA);TempP=N1/N2*100;P=strcat(num2str(TempP),'%') %后验差检验 %计算小误差概率m=input('请输入预测期数:');for g=1:(length(X0)+m)v(1,g)=g-1;endvm=b*exp(-a*v)+c;for j=1:length(v)-1l(1,j)=m(j+1)-m(j);endxyz=[m(1),l]%预测值disp(['小误差概率为：',num2str(P)]);disp(['后验方差比为：',num2str(C)]);disp(['预测值为：',num2str(xyz)]);。

灰色预测及MATLAB实现

3.1灰色预测基础知识
什么是灰色预测？
灰色预测是就灰色系统所做的预测。所谓灰色系统是介于白色系统和黑箱系统之间的过渡系统，其具体的含义是：如果某一系统的全部信息已知为白色系统，全部信息未知为黑箱系统，部分信息已知，部分信息未知，那么这一系统就是灰色系统。一般地说，社会系统、经济系统、生态系统都是灰色系统。例如物价系统，导致物价上涨的因素很多，但已知的却不多，因此对物价这一灰色系统的预测可以用灰色预测方法。灰色系统理论认为对既含有已知信息又含有未知或非确定信息的系统进行预测，就是对在一定方位内变化的、与时间有关的灰色过程的预测。尽管过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此这一数据集合具备潜在的规律，灰色预测就是利用这种规律建立灰色模型对灰色系统进行预测。
ˆ （4）用最小二乘法求解灰参数 a ( a, ) ( B B )
T T
1
B Yn 。
T
ˆ （5）将灰参数 a 代入
ˆ x dy
(1)
dx
(1)
ax
(1)
，求解得
(1)
dt
(t 1) ( x
(0)

a
)e
at

a
dx ˆ ˆ (1) 由于 a 是通过最小二乘法求出的近似值，因此 x (t 1) 事近似表达
(1)
ˆ 序列，得到近似数据序列 x ˆ x
(0)
(0)
ˆ (t 1) x
(1)
ˆ (t 1) x
(t )
（7）建立灰色预测模型进行检验，步骤如下：
① 计算 x
(0)
与x
(0)
(t ) 之间的残差和相对误差 (t ) x

灰色理论灰色预测模型和灰色关联度分析matlab通用代码(精)

%该程序用于灰色关联分析，其中原始数据的第一行为参考序列，1至15行为正相关序列，16至17为负相关序列clc,clearload x.txt %把原始数据存放在纯文本文件x.txt 中%如果全为正相关序列，则将两个循环替换为下列代码%for i=1:size(x,1%x(i,=x(i,/x(i,1;%endfor i=1:15x(i,=x(i,:/x(i,1; %标准化数据endfor i=16:17x(i,:=x(i,1./x(i,:; %标准化数据enddata=x;n=size(data,1;ck=data(1,:;%分离参考序列bj=data(2:n,:;m1=size(bj,1;for j=1:m1t(j,:=bj(j,:-ck;endjc1=min(min(abs(t';jc2=max(max(abs(t';rho=0.5;%灰色关联度为0.5ksi=(jc1+rho*jc2./(abs(t+rho*jc2;r=sum(ksi'/size(ksi,2;r %灰色关联度向量[rs,rind]=sort(r,'descend' %对关联度进行降序排序%该函数用于灰色预测模型，其中x0为列向量，alpha一般取0.5，将第一个数据视为序号为0,k从0开始的序号矩阵function y=huiseyuce(x0,alpha,kn=length(x0;x1=cumsum(x0;for i=2:nz1(i=alpha*x1(i+(1-alpha*x1(i-1;endz1=z1';B=[-z1(2:n,ones(n-1,1];Y=x0(2:n;ab=B\Y;y1=(x0(1-ab(2/ab(1*exp(-ab(1*k+ab(2/ab(1;%产生预测累加生成序列y=[x0(1 diff(y1]%产生灰色预测数据。

分数灰色预测matlab代码详解

分数灰色预测matlab代码详解分数灰色预测（Fractional Grey Model，FGM）是以函数的非线性分数阶微分方程为基础的一种新型时间序列预测方法。

在FGM中，将已知的数值序列$X(k),k=1,2,\cdots,n$ 作为被预测序列，通过构造二阶单位根非线性微分方程，建立起包含一阶导数、二阶导数、初值以及终值的模型，进而计算模型参数并根据模型实现对未来数值的预测。

该方法对中小样本序列具有良好的预测效果，且其原理简单、易于理解，因此在实际预测中得到了广泛应用。

以下是FGM的matlab代码和具体解释。

首先，导入数据，例如一组销售数据：```matlabY = [42 59 72 88 114 138 157 186 212];```接下来，可以通过以下函数将数据转换为一阶、二阶累加数据：```matlabX = cumsum(Y); % 一阶累加n = length(X);for i = 1:n-1DX(i) = X(i+1)-X(i); % 二阶累加end```然后，可以定义分数阶微分方程。

FGM的分数阶微分方程为：$$_{0}D_{t}^{\alpha}X(t)=a_{1}\cdot X(t)+a_{2}$$其中 $a_{1}$、$a_{2}$ 为待求参数，$\alpha$ 为分数阶微分阶次。

根据变换公式，$a_{1}$ 和 $a_{2}$ 可以通过以下函数求得：```matlabsyms a1 a2eq1 = DX(1) - a1*X(1) - a2;eq2 = DX(2) - a1*X(2) - a2;eq3 = DX(3) - a1*X(3) - a2;eq4 = DX(4) - a1*X(4) - a2;eq5 = DX(5) - a1*X(5) - a2;eq6 = DX(6) - a1*X(6) - a2;eq7 = DX(7) - a1*X(7) - a2;eq8 = DX(8) - a1*X(8) - a2;eqn = DX(n-1) - a1*X(n-1) - a2;[eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eqn][A,B] = equationsToMatrix([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eqn],[a1,a2]); X = linsolve(A,B);a1 = X(1)a2 = X(2)```这里使用了符号计算工具箱，利用上述方程组解出 $a_{1}$ 和 $a_{2}$。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

灰色预测专设工⑼他QA—叫吋)为原始数列.其1次累❖加生成数列为恥=妙①曲⑵，…卅何),其中X° 仇)二工* ° (0.址=1=2= -:n5-1卷定义卫的决导数为d(k) = *町(上)=x 叫咼-x cl)(Jt-l).令为数列工①的邻值生成数列.即却(去)=^(*) + (1- a)x0)(t-lX于是定义GM (L 1)的灰微分方程模型为d(k)-血⑴住)=K即或严＞(£) + “尹⑻=人⑴在式(1)中』。

＞(灼称为灰导数，我称为发展系数, 弧称为白化背景值，b称为灰作用量乜将时刻表殳二2「3「/代入(1)式有V!1「—ay=代⑶ B =Ib*- :X闵0)-Z，：］(K)1于是G\I <1»1）複至可表示为Y = Bu.現在问题归结为求sb 在值。

用一元线性回归・即最小二秦法求它们的活计值为注二实陌上回归分析中求估计值是用软件计尊的・有标准程序求解，iOmaClab 等。

GM <1» 1>的白化晏対于G\I <1> 1）的灰微分方程（1） >如果将灰导数打（Q 的时刻视为连绫变里"则x°）视为时问（函数卅⑺,于是*〉（Q 対血于导数里级心2 >白化背臬值申的对应于导数卅⑴。

否朋，対数摄做适当的克换处理■如平移叢换:取C 使得鞍据列严伙）=工⑴伙）+ G 上=1,2,…,的级比都華住可吝禎盖内。

心⑴⑴ + o?i> (r)二◎ dr<2)GM mi ）质色预测的步骤1 •教摇的枪绘与处連为了ftilGAl （1,1）建複方法的可行性，亲要为已知期S 做必要的检蛉处理。

设原始教据列为了逛=（乂°（1）*6（2）严炉00； >计算数列的级比如果所有的级比都落在可容覆盖区间 • fc =A-2,3"・如果対所有的|p 伙)|<0・1 -则认为达到较高的要求,否则若旳所有的|。

(鋼<0・2 >则认为送到一般要求。

2.逢立G\I C1F 1）樓型不妨设少=(由(1)/®(2〉,…¥切)彌足上面的要求，以它为埶据列建立G 、I CU 1)複塑J°)仏)十血⑴(上)=b.用回归分析求得场b 的佔计值，于是相应的白化僅型为絹为dthk *）© = （2 ⑴-？上 YE +三 a a于是得到预则值护住十1)=(少(1)-2)穴+色,上=12…1,aa从而相应地得到预测值：x <0) (^ +1) = x (1)(^ +l)-x o )(t)5 k = 1,2,…,M-1,s.检验预测值(1)死差检验二计尊相对残差如果对所有的 ^<0.1，认为达到技高的要求,否则，若対所有的|f(/t)|<02 >则认为达到一骰要求。

(2)级比偏差值检验；计算皿)=1上0% l+0.5a 律）,•:•例北方某城市1986〜1992年道路交通噪声平均芦级数据见表6 第一步:级比检验建立交通噪声平均声级数粥时间序列如F ：x (0) = (2),…工⑼⑺)-(71.1, 72A 724, 72L71A 72.0,71.6)⑴求级比/住)心少(上-1) • 2 = (A(2)/(3),…,班7» ‘丿 =(0.982.1.1.0042.1.0098. 0.9917.1.0056) (2)级比判断由于所有的入(k)曰0.982,1.0098], k=2,3、，7,故可以用x(0)作满意厉GM (b 1)建模。

第二步:GM (1, 1)建模(1) 对原始数据X®作一次累加，即 X 。

)=(71.1443.5215.9.28&359.4.431.4.503) (2) 构造数据矩阵B 及数据向量Y•:♦表6市近年来交通噪声数据[dB(A)]序号年份旳z1 1996 71.12 1997 72.43 1988 72.4 4 1989 72.15 1990 71.46 1991 72.07 1992 71.6-駅“弋 6 ) + x (l < 7 )(3) ifMif= ^r ^)-1B r y =(0-0023)<4*®7V72.6573dr+ 0.0023 丿〉=72.6573 dt❖求解得沙⑴= (x"°)(l)-0w + £.aar^(Jt +l) = (x^Xl)--)^ +?=-30929&7陀％ +31O QO(5)求生成教列值: i w (上+1)令k=HS4J6由上面的南冏五鈕数可算得，其中取cX 由丘⑴(1) =沁玖得？ 11讣(上)=尹(上)-炉(上-1),G1-1,72.4: 72.1 72.1, 719： 71.7, 71 6)曲= 3(1)，旳(2)，…，旳(7)>=•:・第三步:模型检验•:•模型的各种检验指标値的计算结果见裘7. •:・表7 G 、I(1,1濃型检验表・：・序号年份原始値模型值残差相对误差级比偏差•:• 1 1986 71.1 71.1 0 0❖ 21987 72.4 72.4 -0.0057 0.01%0.0023=b r 1 12> 0 3>•:・3198872.472.20.16380.23%0.0203• 4198972.172.10.03290.05%-0.0018❖5199071.471.9-0.49840.7%-0.0074•:- 6199172.071.70.26990.37%0.0107❖7199271.671.60.0378 0.05%-0.0032•:•经验证，该模型的精度较高，可进行预测和预报。

计算的M 耐程序如下■蚁和预报lamda = x0(l ：n^l)./id )(2：n) % 计算级比 range = minnax( lamda*) %计算锁比的范国 xl =cumsum(xO)； % 累加运算B = f -0^5* (xl(l ;n - 1) +xl(2：n))t ones(n-lj)]； V = xO(2：n) * u = B\Y % 拟合 ^»u(l)^atU (2)=bx=dsolve(^Dx + a*x = b\f x(Q)=池')< 绻戏徴裁方程的挣号餐X = subfi(xJ ,a ,Jb ,/xO ,| J u( 1)+ u(2) ,xO( I) I )； % 代入牯计赛魏值和初始值 yucel = Bubs (x 「t/* [0： n - 1 ]): %求巳他第据的预測值y = vpa(x t 6) %其中曲6表示显示&位教字 yuce = [xO(l)f diff(yucel)] * 差分运勲还原擬据 delta - abs( epsilon % 讣算祀对课產rbo = l-(l-0.5*u(l))41 心 F ⑴厂“京 Y^rhuJst 模型作用：求累加数列、求ab 的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]; format long ; if len gth(x(:,1))==1 x=x'; end% 设置计算精度 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换n=len gth(x); z=0; for i=1: n z=z+x(i,:); be(i,:)=z; %取输入数据的样本量%计算累加值，并将值赋予矩阵 beendfor i=2: n %对原始数列平行移位y(i-1,:)=x(i,:);endfor i=1: n-1%计算数据矩阵 B 的第一列数据 c(i,:)=-0.5*(be(i,:)+be(i+1,:));xO=[71 J 72,4 72 A 72 Jn =length(xO)；71 A 72 .0 71 #6]\%注意疑里为列向epsilon - xO" — yuce % 计算残差屈丫 t 计鼻奴比備垦扰・讥1endfor j=1:n-1 e(j,:)=1;%计算数据矩阵B 的第二列数据endfor i=1:n-1B(i,1)=c(i,:);B(i,2)=e(i,:);%构造数据矩阵Bendalpha=inv(B'*B)*B'*y; for i=1:n+1 值%计算参数矩阵即a b 的值%计算数据估计值的累加数列，如改为n+1 为n+m 可预测后m-1 个ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha( 2,:)/alpha(1,:); %显示输出预测值的累加数列endvar(1,:)=ago(1,:) for i=1:n %显示输出预测值%如改n 为n+m-1 ，可预测后m-1 个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x);c%调用统计工具箱的标准差函数计算后验差的比值Lxago alpha var %显示输出预测值的累加数列%显示输出参数数列%显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b 的值、求预测方程clc,clearx0=[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)=bx=dsolve( 'Dx+a*x=b' , 'x(0)=x0' ); %求微分方程的符号解x=subs(x,{ 'a' , 'b' , 'x0' },{u(1),u(2),x0(1)}) %代入估计参数值和初始值yuce1=subs(x, 't' ,[0:n-1]); %求已知数据的预测值y=vpa(x,6) %其中的6 表示显示6 位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据。