灰色系统新陈代谢模型源码(matlab)

灰⾊预测模型及MATLAB实例下⾯将主要从三⽅⾯进⾏⼤致讲解，灰⾊预测概念及原理、灰⾊预测的分类及求解步骤、灰⾊预测的实例讲解。

⼀、灰⾊预测概念及原理：1.概述：关于所谓的“颜⾊”预测或者检测等，⼤致分为三⾊：⿊、⽩、灰，在此以预测为例阐述。

其中，⽩⾊预测是指系统的内部特征完全已知，系统信息完全充分；⿊⾊预测指系统的内部特征⼀⽆所知，只能通过观测其与外界的联系来进⾏研究；灰⾊预测则是介于⿊、⽩两者之间的⼀种预测，⼀部分已知，⼀部分未知，系统因素间有不确定的关系。

细致度⽐较：⽩>⿊>灰。

2.原理：灰⾊预测是通过计算各因素之间的关联度，鉴别系统各因素之间发展趋势的相异程度。

其核⼼体系是灰⾊模型（Grey Model，GM），即对原始数据做累加⽣成（或者累减、均值等⽅法）⽣成近似的指数规律在进⾏建模的⽅法。

⼆、灰⾊预测的分类及求解步骤：1.GM（1,1）与GM（2,1）、DGM、Verhulst模型的分类⽐较：预测模型适⽤场景涉及的序列GM（1,1）模型⼀阶微分⽅程，只含有1个变量的灰⾊模型。

适⽤于有较强指数规律的序列。

累加序列均值序列GM（2,1）模型适⽤于预测预测具有饱和的S形序列或者单调的摆动发展序列缺陷。

累加序列累减序列均值序列DGM模型累加序列累减序列Verhulst模型累加序列均值序列2.求解步骤思维导图：其中预测过程可能会涉及以下三种序列、⽩化微分⽅程、以及⼀系列检验，由于⼤致都相同，仅仅是某些使⽤累加和累减，⽽另外⼀些则使⽤累加、累减和均值三个序列的差别⽽已。

于是下⾯笔者将对其进⾏归纳总结再进⾏绘制思维导图，帮助读者理解。

（1）原始序列（参考数据列）：（2）1次累加序列（1-AGO）：（3）1次累减序列（1-IAGO ）：（也就是原始序列中，后⼀项依次减去前⼀项的值，例如，[x(2)-x(1),x(3-x(2),...,x(n)-x(n-1))]。

）（4）均值⽣成序列：（这是对累加序列"（前⼀项+后⼀项）/2"得出的结果。

%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。

灰色预测心设尹曲⑴#为原始数列，其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 实现预备知识（1）灰色系统白色系统是指系统内部特征是完全已知的;黑色系统是指系统内部信息完全未知的;而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知,另一部分信息未知或不确定。

(2)灰色预测 灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行 预测。

尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。

灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。

目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的ＧＭ（1，1)模型。

它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。

经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。

因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。

1 灰色系统的模型GM(1,１）１．1 GM(1，1)的一般形式设有变量X (0)={X (0)(i ）,i ＝1,2,．..，n}为某一预测对象的非负单调原始数据列，为建立灰色预测模型:首先对X（0)进行一次累加（１—ＡGO ， Acum ｕl ａt ｅｄ Ge ｎeｒa ｔｉng Opera ｔo ｒ)生成一次累加序列: X (1)={X（1）(k ),k =1，2,…,n}其中X (1）(k )＝∑=ki 1X (0）(i)＝X (1)(ｋ-１）+ X （0）(ｋ) (1)对X（1)可建立下述白化形式的微分方程：dtdX )1(十)1(aX =u (2)即G Ｍ（１，1)模型。

上述白化微分方程的解为(离散响应): ∧X（1)(k +１)＝(X (０)(1)－a u ）ake -+au (3）或∧X （1）（k )=(X (0）(1)-a u ))1(--k a e +au （４)式中：k 为时间序列，可取年、季或月。

matlab灰色预测模型函数Matlab是一种广泛应用于科学计算和工程领域的软件工具，它提供了许多函数和工具箱，用于数据分析和建模。

其中一个重要的函数是灰色预测模型函数，它可以用来预测和分析时间序列数据。

本文将介绍灰色预测模型函数的原理和应用，并通过一个示例来演示其使用方法。

灰色预测模型是一种基于灰色系统理论的预测方法，它适用于样本数据较少、不完整或不规律的情况。

在灰色预测模型中，数据被分为两类：发展数据和规律数据。

发展数据是指拥有较完整信息的数据，规律数据是指缺乏完整信息的数据。

通过对规律数据进行处理和建模，可以预测未来数据的趋势和变化。

在Matlab中，可以使用灰色预测模型函数进行数据预测和分析。

该函数可以通过输入历史数据和需要预测的时间步长来生成预测结果。

具体的使用方法如下：1. 导入数据：首先需要导入需要预测的时间序列数据。

可以使用Matlab中的数据导入函数来读取数据文件或手动输入数据。

2. 数据预处理：对导入的数据进行预处理，包括去除异常值、平滑数据等。

可以使用Matlab中的数据处理函数来完成这些操作。

3. 构建灰色模型：使用灰色预测模型函数来构建灰色模型。

该函数需要输入历史数据和需要预测的时间步长。

4. 模型评估：对构建的灰色模型进行评估，包括计算预测误差、拟合度等指标。

可以使用Matlab中的统计函数来完成这些计算。

5. 预测结果：根据构建的灰色模型，可以生成未来数据的预测结果。

可以使用Matlab中的预测函数来完成这一步骤。

下面通过一个示例来演示灰色预测模型函数的使用方法。

假设我们有一个销售数据的时间序列，我们希望预测未来三个月的销售额。

首先，我们需要导入销售数据，并进行数据预处理。

然后，我们可以使用灰色预测模型函数来构建灰色模型。

最后，我们可以通过预测函数来生成未来三个月的销售额预测结果。

在实际操作中，我们可能还需要对模型进行参数调优和模型选择。

可以使用交叉验证等方法来选择最优的参数和模型。

function GM1_1(X0)%format long ;[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) ;P_t=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.5 R=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),'%') %后验差检验%计算小误差概率调用例子：X0=[2.874,3.278,3.337,3.39,3.679];GM1_1(X0)参考资料：我自己想的，算法不是很优，希望你自己改进吧。

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的使用，对于进行灰色系统建模具有重要意义。

二、数据预处理在进行灰色系统建模之前，需要对原始数据进行预处理，以提高后续建模的准确性和可靠性。

数据预处理包括数据清洗、数据平滑和数据归一化等步骤。

Matlab 提供了丰富的数据处理函数和工具，可以快速、灵活地完成这些操作。

1. 数据清洗数据清洗是指删除或修正含有噪声、异常值或缺失值的数据。

Matlab中可以使用滤波函数、插值函数和替换函数等方法对数据进行清洗。

例如，可以使用median函数对数据进行中值滤波，去除噪声干扰。

另外，使用interp1函数进行数据插值，可以填补缺失值，使数据更加完整。

2. 数据平滑数据平滑是指通过降低数据的波动性，使其更具有连续性和稳定性。

Matlab中常用的数据平滑方法包括移动平均法、指数平滑法和小波平滑法等。

移动平均法通过计算滑动窗口内数据的平均值，来平滑原始数据。

指数平滑法则采用指数加权平均的方式，对数据进行平滑处理。

小波平滑法则利用小波分析的方法，对数据进行平滑处理。

3. 数据归一化数据归一化是指将不同量纲或取值范围的数据，转换为统一的尺度。

常用的归一化方法包括最小-最大归一化和Z-score归一化等。

最小-最大归一化将数据线性映射到[0,1]的范围内，使数据具有统一的尺度和可比性。

Z-score归一化则通过计算数据与均值的偏差，除以标准差，将数据标准化为均值为0，标准差为1的分布。

灰色关联度matlab源程序（完整版）近几天一直在写算法，其实网上可以下到这些算法的源程序的，但是为了搞懂，最搞清楚，还是自己一个一个的看了，写了，作为自身的积累，而且自己的的矩阵计算类库也迅速得到补充，以后关于算法方面，基本的矩阵运算不用再重复写了，挺好的，是种积累，下面把灰关联的matlab程序与大家分享。

灰色关联度分析法是将研究对象及影响因素的因子值视为一条线上的点,与待识别对象及影响因素的因子值所绘制的曲线进行比较,比较它们之间的贴近度,并分别量化,计算出研究对象与待识别对象各影响因素之间的贴近程度的关联度,通过比较各关联度的大小来判断待识别对象对研究对象的影响程度。

简言之，灰色关联度分析的意义是指在系统发展过程中，如果两个因素变化的态势是一致的，即同步变化程度较高，则可以认为两者关联较大；反之，则两者关联度较小。

因此，灰色关联度分析对于一个系统发展变化态势提供了量化的度量，非常适合动态(Dynamic)的历程分析。

灰色关联度可分成“局部性灰色关联度”与“整体性灰色关联度”两类。

主要的差别在于局部性灰色关联度有一参考序列，而整体性灰色关联度是任一序列均可为参考序列。

关联度分析是基于灰色系统的灰色过程, 进行因素间时间序列的比较来确定哪些是影响大的主导因素, 是一种动态过程的研究。

关联度计算的预处理，一般初值化或者均值化，根据我的实际需要，本程序中使用的是比较序列与参考序列组成的矩阵除以参考序列的列均值等到的，当然也可以是其他方法。

%注意：由于需要，均值化方法采用各组值除以样本的各列平均值clear;clc;yangben=[47.924375 25.168125 827.4105438 330.08875 1045.164375 261.37437516.3372 6.62 940.2824 709.2752 962.1284 84.87455.69666667 30.80333333 885.21 275.8066667 1052.42 435.81]; %样本数据fangzhen=[36.27 14.59 836.15 420.41 1011.83 189.5464.73 35.63 755.45 331.32 978.5 257.8742.44 23.07 846 348.05 1025.4 296.6959.34 39.7 794.31 334.63 1016.4 317.2752.91 17.14 821.79 306.92 1141.94 122.044.21 4.86 1815.52 2584.68 963.61 0.006.01 2.43 1791.61 2338.17 1278.08 30.873.01 1.58 1220.54 956.14 1244.75 3.9125.65 7.42 790.17 328.88 1026.01 92.82115.80 27 926.5 350.93 1079.49 544.3812.63 8.75 1055.50 1379.00 875.10 1.65]; %待判数据[rows,cols]=size(fangzhen);p=0.5; %分辨系数[m,n]=size(yangben);R=[];for irow=1:rowsyy=fangzhen(irow,:);data=[yy;yangben];data_gyh1=mean(yangben)for i=1:m+1for j=1:ndata_gyh(i,j)=data(i,j)/data_gyh1(j);endendfor i=2:m+1for j=1:nDij(i-1,j)=abs(data_gyh(1,j)-data_gyh(i,j));endendDijmax=max(max(Dij));Dijmin=min(min(Dij));for i=1:mfor j=1:nLij(i,j)=(Dijmin+p*Dijmax)/(Dij(i,j)+p*Dijmax); endendLijRowSum=sum(Lij');for i=1:mRij(i)=LijRowSum(i)/n;endR=[R;Rij];endRmatlab求灰色关联度矩阵源代码2010-12-11 22:57 function greyrelationaldegree(X,c)%GRAYRELATIONALDEGREE this function is used for calculating the gery %relation between squence%rememeber that the first column of the input matrix is the desicion %attribution squences.what we want to calculate is the grey ralational degree between%it and other attributions%X is the squence matrix, c is the parameter used in the function%in most of the time, the value of c is 0.5firstrow = X(1,:);reci_firstrow = 1./firstrow;reci_convert = diag(reci_firstrow);initialMIRROR = X*reci_convert;% find the initial value mirror of the sequce matrixA = initialMIRROR'[nrow,ncolumn] = size(A);for (i=2:nrow)C = A(i,:)-A(1,:)D=abs(C);eval(['B' num2str(i) '=D']);amax = max(eval(['B' num2str(i)]))amin = min(eval(['B' num2str(i)]))maxarray(i-1)=amaxminarray(i-1)=aminend %find the difference squence and the max value and min value of each squencemaxmax = max(maxarray)minmin = min(minarray)for(i=2:nrow)for(j=1:ncolumn)eval(['greyrelationdegree' num2str(i)'(j)=(minmin+c*maxmax)/(B' num2str(i) '(j)+c*maxmax)'])endend % calculate the greyralational degree of each datafor(i=2:nrow)eval(['greyrelatioanaldegree_value' num2str(i) '= mean (greyrelationdegree' num2str(i) ')' ])end基于matlab灰色关联度计算的实现2006年07月28日星期五上午 11:06 function r=incident_degree(x0,x1)%compute the incident degree for grey model.%Designed by NIXIUHUI,DalianFisherUniversity.%17 August,2004,Last modified by NXH at 21 August,2004%数据初值化处理x0_initial=x0./x0(1);temp=size(x1);b=repmat(x1(:,1),[1 temp(2)]);x1_initial=x1./b;%分辨系数选择K=0.1;disp('The grey interconnect degree is: ');x0_ext=repmat(x0_initial,[temp(1) 1]);contrast_mat=abs(x0_ext-x1_initial);delta_min=min(min(contrast_mat));%delta_min在数据初值化后实际为零delta_max=max(max(contrast_mat));a=delta_min+K*delta_max;incidence_coefficient=a./(contrast_mat+K*delta_max);%得到关联系数r=(sum(incidence_coefficient'))'/temp(2); %得到邓氏面积关联度我们根据图1的步骤和图2的数据进行编程实现，程序如下：%清除存空间等clear;close all;clc;%载入源数据 %其实这里可以载入execl表格的n=15; %参与评价的人数m=4; %参与评价的指标个数X_0=zeros(n,m); % 数据矩阵X_2=zeros(n,m); %偏差结果的求取矩阵X_3=zeros(n,m); % 相关系数计算矩阵a1_0=[13 18 17 18 17 17 18 17 13 17 18 13 18 13 18];a2_0=[18 18 17 17 18 13 17 13 18 13 17 13 13 17 17];a3_0=[48.67 43.33 43.56 41.89 39.47 43.44 37.97 41.14 39.67 39.83 34.11 40.58 34.19 30.75 21.22];a4_0=[10 10.7 3 5.4 5.4 0.7 4.2 0.5 9.3 0.85 2.9 5.45 4.2 2.7 6]; %指标数X_1=[a1_0',a2_0',a3_0',a4_0']; %最后使用到的数据矩阵%1 寻找参考列x0=[max(a1_0),max(a2_0),max(a3_0),max(a4_0)]; %取每列的最大值（指标的最大值）%2 计算偏差结果i=1;while(i~=m+1) %为什么这个地方会出问题呢for j=1:1:nX_2(j,i)=abs(X_1(j,i)-x0(i));end;i=i+1;end%3 确定偏差的最值error_min=min(min(X_2));error_max=max(max(X_2));%4 计算相关系数i=1;p=0.5;while(i~=m+1)for j=1:1:nX_3(j,i)=(error_min+p*error_max)/(X_2(j,i)+p*error_max); end;i=i+1;end%X_3 %可以在此观察关联矩阵%5 计算各个学生的关连序a=zeros(1,n);for j=1:1:nfor i=1:1:ma(j)=a(j)+X_3(j,i); %%%%其实可以直接用sumend;a(j)=a(j)/m; %%%%%%%%%可以改进%%%%%%%%%%123下一页%end%a %在此可以观测各个学生的序%改进:如果各个指标的所占权重不一样的话,可以添加相应的权系数%6 排序b=a';[c,s]=sort(b);for i=1:1:nd(i)=i;endd=d';result=[d b c s]%7 将结果显示出来figure(1);plot(a);figure(2)bar(a); %柱状图最后所得到的结果如图3到图5所示。

单位：亿请根据1978-2012年的全国人口数据，用GM(1,1)来预测未来人口数量，能预测50年内的数据吗？100年呢？为什么？摘要：关键字：人口预测、GM(1,1)一、问题的重述二、问题的分析三、模型假设1、假设本问题所使用的数据均真实有效，具有统计分析价值2、不考虑战争、灾害、疾病对人口数目的影响；3、...4、...四、符号说明与名词解释4.1符号说明见模型假设。

4.2名词解释...五、模型的建立与求解5.1模型建立在灰色系统理论中，称抽象的逆过程为灰色模型，也称GM。

它是根据关联度、生成数灰导数，灰微分等观点和一系列数学方法建立起来的连续型的微分方程。

通常GM表示为GM(n,h)。

当n=h=1时即构成了单变量一阶灰色预测模型。

设原始时间序列为：(0)(0)(0)(0)(0)[(1),(2),(3),,()]X x x x x n=L设(1)X为(0)X的一次累加序列：即：(1)(0)(1)(1)(1)(1)(0)()(-1)(),2,, x xx k x k x k k n ⎧=⎪⎨⎪=+=⎩L得(1)(1)(1)(1)(1)[(1),(2),(3),,()] X x x x x n=L利用(1)X 计算GM(1,1)模型参数a 、u ，令ˆ[,]T aa u =, 则有1ˆ()T T N aB B B Y -= 其中(1)(1)(1)(1)(1)(1)1((1)(2))121((2)(3))121((1)())12x x x x B x n x n ⎡⎤-+⎢⎥⎢⎥⎢⎥-+⎢⎥=⎢⎥⎢⎥⎢⎥--+⎢⎥⎣⎦L L (0)(0)(0)[(2),(3),,()]T N Y x x x n =L由此获得：(1)(1)ˆ(1)((1))ak u uxk x e a a -+=-+ 于是：(0)(1)(1)ˆˆˆ(1)(1)()1xk x k x k +=+-（）或(0)(0)ˆ(1)-((1))ak uxk a x e a-+=-（2）注：（1）式是根据(0)X 和(1)X 的关系的到的（2）式是利用数学求导还原得到的至于用哪个，最好看相对误差5.2模型求解：求解得：六、模型检验6.1残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验。

function g(x); %定义函数gm(x)for j=1:10clc; %清屏,以使计算结果独立显示format long; %设置计算精度if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';end;n=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; %计算参数α、μ矩阵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,:)=var(i,:)-x(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值cago %显示输出预测值的累加数列alpha %显示输出参数α、μ数列var %显示输出预测值error %显示输出误差c %显示后验差的比值cvar=fix(var);for i=2:11x(i-1)=var(i);end。