灰色预测模型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)] %差分运算，还原数据。

GM(1,1)灰色预测模型IntroductionInitial给定原始序列：x(0) =(x(0)(1), x(0)(2), x(0)(3)…, x(0)(n))Step 1一次AGO(1-AGO)生成序列，以弱化原始序列的随机性和波动性：x(1) =(x(1)(1), x(1)(2), x(1)(3)…, x(1)(n)) Matlab Programclearsyms a b;c=[a b]';fid=fopen('.\Grey Model\test.txt');x0=fscanf(fid,'%f');x0=x0';fclose(fid);x1=cumsum(x0); %原始数据累加n=length(x0);for i=1:(n-1)z(i)=(x1(i)+x1(i+1))/2; %生成累加矩阵end%计算待定参数的值Y=x0;Y(1)=[];Y=Y';B=[-z;ones(1,n-1)];B=B';c=inv(B'*B)*B'*Y;c=c';a=c(1);b=c(2);%预测后续数据%预测之后10个时间单位的数据xx1=[];xx1(1)=x0(1);for i=2:(n+10)xx1(i)=(x0(1)-b/a)/exp(a*(i-1))+b/a; endxx0=[];xx0(1)=x0(1);Step 2(1) dx (1)dt+ax (1)(t )=u ，式中a, u 为待定系数。

灰微分方程模型为：x (0)(k )+az (1)(k )=u ，z 为背景值z (1)(k )=1/2(x (1)(k )+x (1)(k −1))(2) 构造矩阵B 和数据向量Y nY n =Ba ̂Y n =[ x (0)(2)x (0)(3)⋮x (0)(n )] , B =[ −1/2(x (1)(1)+x (1)(2)),−1/2(x (1)(2)+x (1)(3)),⋮−1/2(x (1)(n −1)+x (1)(n )), 1 1 ⋮ 1]a ̂=(au)=(B T B)−1B T Y nStep 3模型响应函数x ̂(1)(k +1)=(x (0)(1)−u a )e −ak +u ax ̂(0)(k +1)=x ̂(1)(k +1)−x ̂(1)(k )Step 4检验和判断GM(1,1)模型的精度 (1) 残差检验for i=2:(n+10)xx0(i)=xx1(i)-xx1(i-1); end%关联度检验 for i=1:ne(i)=abs(x0(i)-xx0(i)); endmmax=max(e); for i=1:nee(i)=0.5*mmax/(e(i)+0.5*mmax); endr=sum(ee)/n; %后验差检验x0bar=sum(x0)/n; s1=0; for i=1:ns1=s1+(x0(i)-x0bar)^2; ends1=sqrt(s1/n); s2=0;ebar=sum(e)/n; for i=1:ns2=s2+(e(i)-ebar)^2; ends2=sqrt(s2/n); C=s2/s1; p=0;for i=1:nif abs(e(i)-ebar)<0.6745*s1绝对误差：ε(k)=|x(0)(k)−x̂(0)(k)|相对误差：Φ(k)=ε(k)x(0)(k)(2) 关联度检验分辨率β一般取0.5，此时若关联度大于0.6则认为模型可接受(3) 后验差检验和小误差概率原始序列标准差：S1=√∑[x(0)(i)−x̅(0)]2n绝对误差序列标准差：S2=√∑[ε(i)−ε̅]2n计算方差比：C=S2S1小误差概率：P=P{|ε(i)−ε̅|<0.6745S1}p=p+1;endendp=p/n;Cpif p>0.95&C<0.35disp('预测精度好');else if p>0.8&C<0.5disp('预测合格');else if p>0.7&C<0.65disp('预测勉强合格'); elsedisp('预测不合格'); endendend%原始数据与预测数据进行比较t1=1:n;t2=1:(n+10);xx0plot(t1,x0,'o',t2,xx0)。

灰色预测专设工⑼他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 >白化背臬值申的对应于导数卅⑴。

于是G\I （1，1）的坝徽 分方樂対应于的白微分方程为内・则数堀列X©可以塗互G\I <19 1） 且可以进行页色预测。

否朋，対数摄做适当的克换处理■如平移叢换:取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 -则认为达到较高的要求,否则 若旳所有的|。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

综上所述，以上就是分数灰色预测的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,n）（灰色模型代码）%灰色预测模型GM（1，n）模型的matlab源代码，包括预测模型的建立，以及模型的精度检验指标c,p的计算%假设预测3步，N=3%如在命令窗口键入：%gm=ycgm1n([1.6,1.7,2,1.8,1.9],[2,2.4,3,3.2,3.1],[3,3.1,3.2,3.5 ,2.8],3)function GM=ycgm1n(data1,data2,data3,N) %data1:纵摇，data2:升沉，data3:波浪T=length(data1);PYX1=data1;PYX2=data2;PYX3=data3;%进行数据预处理，这里用初值化X0_1=PYX1./PYX1(1);X0_2=PYX2./PYX2(1);X0_3=PYX3./PYX3(1);%用AGO生成一阶累加生成模块X1_1(1)=X0_1(1);X1_2(1)=X0_2(1);X1_3(1)=X0_3(1);for i=2:TX1_1(i)=X1_1(i-1)+X0_1(i);X1_2(i)=X1_2(i-1)+X0_2(i);X1_3(i)=X1_3(i-1)+X0_3(i);end%构造累加矩阵Bfor i=1:T-1M1(i)=(0.5*(X1_1(i)+X1_1(i+1)));M2(i)=(0.5*(X1_2(i)+X1_2(i+1)));M3(i)=(0.5*(X1_3(i)+X1_3(i+1)));endB1=zeros(T-1,3);for i=1:(T-1)B1(i,1)=-M1(i); %-(X1_1(i)+X1_1(i+1)))/2; B1(i,2)=X1_2(i+1);B1(i,3)=X1_3(i+1);endB2=zeros(T-1,2);for i=1:(T-1)B2(i,1)=-M2(i); %-(X1_2(i)+X1_2(i+1)))/2; B2(i,2)=X1_3(i+1);endB3=zeros(T-1,2);for i=1:(T-1)B3(i,1)=-M3(i); %-(X1_3(i)+X1_3(i+1)))/2; B3(i,2)=1;endsave B1 B1;save B2 B2;save B3 B3;%构造常数项向量Yfor i=2:TY1(i-1)=X0_1(i);Y2(i-1)=X0_2(i);Y3(i-1)=X0_3(i);endHCS1=inv(B1'*B1)*B1'*Y1'; %用最小二乘法求灰参数HCS1 H1=HCS1'; %H1=[a,b2,b3]HCS2=inv(B2'*B2)*B2'*Y2'; %用最小二乘法求灰参数HCS2 H2=HCS2'; %H2=[a,b3]HCS3=inv(B3'*B3)*B3'*Y3'; %用最小二乘法求灰参数HCS3 H3=HCS3'; %H3=[b,a]%计算出X3的累加序列for i=1:T+NYCX13(i)=(X0_3(1)-H3(2)/H3(1))*exp(-1*H3(1)*(i-1))+H3(2)/H3(1);endfor i=2:T+N% K(i)=XR1(i)-XR1(i-1);YCX0_3(i)=YCX13(i)-YCX13(i-1);endYCX0_3(1)=X0_3(1);%对参数作alpha，beta变换H2=H2./(1+0.5*H2(1));%还原计算出X2的预测值YCX0_2(1)=X0_2(1);for i=2:TYCX0_2(i)=H2(2).*X1_3(i)-H2(1).*X1_2(i-1);endYCX12(T)=X1_2(T);for i=T+1:T+NYCX0_2(i)=H2(2).*YCX13(i)-H2(1).*YCX12(i-1);YCX12(i)=YCX0_2(i)+YCX12(i-1);end%对参数作alpha，beta变换H1=H1./(1+0.5*H1(1));%还原计算出X1的预测值YCX0_1(1)=X0_1(1);for i=2:TYCX0_1(i)=H1(2).*X1_2(i)+H1(3).*X1_3(i)-H1(1).*X1_1(i-1);endYCX11(T)=X1_1(T);for i=T+1:T+NYCX0_1(i)=H1(2).*YCX12(i)+H1(3).*YCX13(i)-H1(1).*YCX11(i-1);YCX11(i)=YCX0_1(i)+YCX11(i-1);end%数据还原GM=YCX0_1; %.*PYX1(1);save GM GM;e0(1,T-1)=zeros;for i=1:T-1 %求X1到X5的残差值e0e0(i)=(X0_1(i+1)-YCX0_1(i+1))/X0_1(i+1); %1-YCX0_1(i+1)/X0_1(i+1);endsave e0 e0;e0_average=sum(abs(e0))/length(e0)p=1-e0_average;X_average=mean(X0_1) %求原始数据x0均值s1=std(PYX1) %求原始数据的标准差s2=std(e0)c=s2/s1 %计算方差比c,c<0.35为好end。

灰色理论预测模型及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))]);。

数学建模-灰⾊预测模型GM（1,1）_MATLAB %GM(1,1).m%建⽴符号变量a(发展系数)和b(灰作⽤量)syms a b;c = [a b]';%原始数列 AA = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];%填⼊已有的数据列！n = length(A);%对原始数列 A 做累加得到数列 BB = cumsum(A);%对数列 B 做紧邻均值⽣成for i = 2:nC(i) = (B(i) + B(i - 1))/2;endC(1) = [];%构造数据矩阵B = [-C;ones(1,n-1)];Y = A; Y(1) = []; Y = Y';%使⽤最⼩⼆乘法计算参数 a(发展系数)和b(灰作⽤量)c = inv(B*B')*B*Y;c = c';a = c(1);b = c(2);%预测后续数据F = []; F(1) = A(1);for i = 2:(n+10) %这⾥10代表向后预测的数⽬，如果只预测⼀个的话为1F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;end%对数列 F 累减还原,得到预测出的数据G = []; G(1) = A(1);for i = 2:(n+10) %10同上G(i) = F(i) - F(i-1); %得到预测出来的数据enddisp('预测数据为：');G%模型检验H = G(1:10); %这⾥的10是已有数据的个数%计算残差序列epsilon = A - H;%法⼀：相对残差Q检验%计算相对误差序列delta = abs(epsilon./A);%计算相对误差Qdisp('相对残差Q检验：')Q = mean(delta)%法⼆：⽅差⽐C检验disp('⽅差⽐C检验：')C = std(epsilon, 1)/std(A, 1)%法三：⼩误差概率P检验S1 = std(A, 1);tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);disp('⼩误差概率P检验：')P = length(tmp)/n%绘制曲线图t1 = 1995:2004;%⽤⾃⼰的，如1 2 3 4 5...t2 = 1995:2014;%⽤⾃⼰的，如1 2 3 4 5... plot(t1, A,'ro'); hold on;plot(t2, G, 'g-');xlabel('年份'); ylabel('污⽔量/亿吨');legend('实际污⽔排放量','预测污⽔排放量'); title('长江污⽔排放量增长曲线'); %都⽤⾃⼰的grid on;。

灰色预测作用：求累加数列、求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,:)/alp ha(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 %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求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)] %差分运算，还原数据。

%x=[1019,1088,1324,1408,1601];gm1(x);测试数据%二次拟合预测GM(1,1) 模型function gmcal=gm1(x)if nargin==0x=[1019,1088,1324,1408,1601]end format long g sizex=length(x); %求数组长度k=0;for y1=x k=k+1; if k>1x1(k)=x1(k-1)+x(k);% 累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1 维数减1，用于计算B yn1(k-1)=x(k);elsex1(k)=x(k);end end %x1,z1,k,yn1 sizez1=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,rightb ra,'+ ',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)%输出时间响应方程for 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);endend endx31fcast %一次拟合预测值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);endendend%二次拟合预测值%***精度检验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);%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 =运行结果x =ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929) ans = x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685) x31fcast =Columns 1 through 3Columns 4 through 61429.691537401951609.90061644041 1812.824603777821019 1088 1324 1408 16011019 1088 13241408 160110191122.89347857097 1264.43142178303 Columns 4 through 61423.80987235488 1603.27758207442 1805.36675232556x41fcast =Columns 1 through 310191118.05685435129 1269.65470492098Cval =0.139501578334155 pval =1。

GM(,n)(灰色模型代码)%灰色预测模型GM（1，n）模型的matlab源代码，包括预测模型的建立，以及模型的精度检验指标c,p的计算%假设预测3步，N=3%如在命令窗口键入：%gm=ycgm1n([1.6,1.7,2,1.8,1.9],[2,2.4,3,3.2,3.1],[3,3.1,3.2,3.5,2.8],3)function GM=ycgm1n(data1,data2,data3,N) %data1:纵摇，data2:升沉，data3:波浪T=length(data1);PYX1=data1;PYX2=data2;PYX3=data3;%进行数据预处理，这里用初值化X0_1=PYX1./PYX1(1);X0_2=PYX2./PYX2(1);X0_3=PYX3./PYX3(1);%用AGO生成一阶累加生成模块X1_1(1)=X0_1(1);X1_2(1)=X0_2(1);X1_3(1)=X0_3(1);for i=2:TX1_1(i)=X1_1(i-1)+X0_1(i);X1_2(i)=X1_2(i-1)+X0_2(i);X1_3(i)=X1_3(i-1)+X0_3(i);end%构造累加矩阵Bfor i=1:T-1M1(i)=(0.5*(X1_1(i)+X1_1(i+1)));M2(i)=(0.5*(X1_2(i)+X1_2(i+1)));M3(i)=(0.5*(X1_3(i)+X1_3(i+1)));endB1=zeros(T-1,3);for i=1:(T-1)B1(i,1)=-M1(i); %-(X1_1(i)+X1_1(i+1)))/2;B1(i,2)=X1_2(i+1);B1(i,3)=X1_3(i+1);endB2=zeros(T-1,2);for i=1:(T-1)B2(i,1)=-M2(i); %-(X1_2(i)+X1_2(i+1)))/2;B2(i,2)=X1_3(i+1);endB3=zeros(T-1,2);for i=1:(T-1)B3(i,1)=-M3(i); %-(X1_3(i)+X1_3(i+1)))/2;B3(i,2)=1;endsave B1 B1;save B2 B2;save B3 B3;%构造常数项向量Yfor i=2:TY1(i-1)=X0_1(i);Y2(i-1)=X0_2(i);Y3(i-1)=X0_3(i);endHCS1=inv(B1&#39;*B1)*B1&#39;*Y1&#39;; %用最小二乘法求灰参数HCS1H1=HCS1&#39;; %H1=[a,b2,b3]HCS2=inv(B2&#39;*B2)*B2&#39;*Y2&#39;; %用最小二乘法求灰参数HCS2H2=HCS2&#39;; %H2=[a,b3]HCS3=inv(B3&#39;*B3)*B3&#39;*Y3&#39;; %用最小二乘法求灰参数HCS3H3=HCS3&#39;; %H3=[b,a]%计算出X3的累加序列for i=1:T+NYCX13(i)=(X0_3(1)-H3(2)/H3(1))*exp(-1*H3(1)*(i-1))+H3(2)/H3(1);endfor i=2:T+N% K(i)=XR1(i)-XR1(i-1);YCX0_3(i)=YCX13(i)-YCX13(i-1);endYCX0_3(1)=X0_3(1);%对参数作alpha，beta变换H2=H2./(1+0.5*H2(1));%还原计算出X2的预测值YCX0_2(1)=X0_2(1);for i=2:TYCX0_2(i)=H2(2).*X1_3(i)-H2(1).*X1_2(i-1);endYCX12(T)=X1_2(T);for i=T+1:T+NYCX0_2(i)=H2(2).*YCX13(i)-H2(1).*YCX12(i-1);YCX12(i)=YCX0_2(i)+YCX12(i-1);end%对参数作alpha，beta变换H1=H1./(1+0.5*H1(1));%还原计算出X1的预测值YCX0_1(1)=X0_1(1);for i=2:TYCX0_1(i)=H1(2).*X1_2(i)+H1(3).*X1_3(i)-H1(1).*X1_1(i-1);endYCX11(T)=X1_1(T);for i=T+1:T+NYCX0_1(i)=H1(2).*YCX12(i)+H1(3).*YCX13(i)-H1(1).*YCX11(i-1);YCX11(i)=YCX0_1(i)+YCX11(i-1);end%数据还原GM=YCX0_1; %.*PYX1(1);save GM GM;e0(1,T-1)=zeros;for i=1:T-1 %求X1到X5的残差值e0e0(i)=(X0_1(i+1)-YCX0_1(i+1))/X0_1(i+1); %1-YCX0_1(i+1)/X0_1(i+1); endsave e0 e0;e0_average=sum(abs(e0))/length(e0)p=1-e0_average;X_average=mean(X0_1) %求原始数据x0均值s1=std(PYX1) %求原始数据的标准差s2=std(e0)c=s2/s1 %计算方差比c,c&lt;0.35为好end。

灰色预测模型matlab程序精确版(推荐完整)编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（灰色预测模型matlab程序精确版(推荐完整)）的内容能够给您的工作和学习带来便利。

同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。

本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快业绩进步，以下为灰色预测模型matlab程序精确版(推荐完整)的全部内容。

灰色预测模型matlab程序精确版（推荐完整)编辑整理：张嬗雒老师尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布到文库，发布之前我们对文中内容进行仔细校对,但是难免会有疏漏的地方,但是我们任然希望灰色预测模型matlab程序精确版(推荐完整) 这篇文档能够给您的工作和学习带来便利.同时我们也真诚的希望收到您的建议和反馈到下面的留言区,这将是我们进步的源泉,前进的动力。

本文可编辑可修改，如果觉得对您有帮助请下载收藏以便随时查阅，最后祝您生活愉快业绩进步,以下为〈灰色预测模型matlab程序精确版(推荐完整）〉这篇文档的全部内容。

灰色预测模型matlab程序灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性%下面程序是灰色模型GM（1,1)程序二次拟合和等维新陈代谢改进预测程序，matla b6.5 ,使用本程序请注明,程序存储为gm1。

m%x = [5999，5903，5848,5700,7884］；gm1(x); 测试数据%二次拟合预测GM（1,1）模型function gmcal=gm1（x）sizexd2 = size（x,2）;%求数组长度k=0;for y1=xk=k+1；if k〉1x1(k）=x1(k—1)+x（k);％累加生成z1（k-1)=-0。

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代码详解分数灰色预测（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}$。