灰色预测程序

灰色预测建模步骤

灰色预测1.理论灰色理论认为系统的行为现象尽管是朦胧的，数据是复杂的，但它毕竟是有序的，是有整体功能的。

灰数的生成，就是从杂乱中寻找出规律。

同时，灰色理论建立的是生成数据模型，不是原始数据模型，因此，灰色预测是一种对含有不确定因素的系统进行预测的方法。

灰色预测是一种对含有不确定因素的系统进行预测的方法。

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

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

2.分类①灰色时间序列预测；即用观察到的反映预测对象特征的时间序列来构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

②畸变预测；即通过灰色模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。

③系统预测；通过对系统行为特征指标建立一组相互关联的灰色预测模型，预测系统中众多变量间的相互协调关系的变化。

④拓扑预测；将原始数据作曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点3.关联度提出系统的关联度分析方法，是对系统发展态势的量化比较分析。

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

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

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

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

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

关系式记x(0)为原始数列x(0)=(x(0)(k)xk=1,2,…，n)=(x(0）⑴，x(0）⑵，…，x(0)(n))记x⑴为生成数列x⑴=(x⑴（k)xk=1,2，…，n)=(x⑴⑴，x⑴⑵，…，x⑴（n))如果x(0) 与x⑴之间满足下列关系，即称为一次累加生成5.建模步骤a、建模机理b、把原始数据加工成生成数；c、对残差（模型计算值与实际值之差）修订后，建立差分微分方程模型；d、基于关联度收敛的分析；e、gm模型所得数据须经过逆生成还原后才能用。

灰色系统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

灰色预测法GM(1,1)总结

灰色预测模型一、灰色预测的概念1. 灰色预测法是一种对含有不确定因素的系统进行预测的方法。

灰色系统是介于白色系统和黑色系统之间的一种系统。

灰色系统内的一部分信息是已知的，另一部分信息时未知的，系统内各因素间具有不确定的关系。

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

尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此可以通过对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

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

二、灰色预测的类型1. 灰色时间序列预测；即用观察到的反映预测对象特征的时间序列来构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

2. 畸变预测；即通过灰色模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。

3. 系统预测；通过对系统行为特征指标建立一组相互关联的灰色预测模型，预测系统中众多变量间的相互协调关系的变化。

4. 拓扑预测；将原始数据作曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点三、GM （1，1）模型的建立 1. 数据处理为了弱化原始时间序列的随机性，在建立灰色预测模型之前，需先对原始时间序列进行数据处理，经过数据处理后的时间序列即称为生成列。

i. 设()()()()()()()()(){},,, (00000)123X X X X X n = 是所要预测的某项指标的原始数据，计算数列的级比()()()(),,,,()00123X t t t n X t λ-==。

如果绝大部分的级比都落在可容覆盖区间(,)2211n n ee-++内，则可以建立GM(1,1)模型且可以进行灰色预测。

灰色预测(由过去预测未来)

end
a
fprintf('参数a=');alpha(1)
fprintf('参数u=');alpha(2)
fprintf('预测未来十年干流6类为');
for i=11:20
a(i-10,:)= var(:,i);
end
lamda=x(1:n-1)./x(2:n);%计算级比
range=minmax(lamda');%级比范围
rho=1-(1-0.5*alpha(1))/(1+0.5*alpha(1))*lamda;%计算级比偏差值，alpha(1)=a
c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c
format long; %设置计算精度
if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换源自 x=x'; end
n=length(x); %取输入数据的样本量
x1=cumsum(x);%累加
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
y=x(2:n);
alpha=B\y;
%预测数据的累加，共n+1个
for i=1:n+10
ago(i)=(x(1)-alpha(2)/alpha(1))*exp(-alpha(1)*(i-1))+alpha(2)/alpha(1);
end
%此处共预测n+1个值
%估计值第一个与实际值相同
var(1,:)=ago(1,:);
%开始第二个估计值

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

灰色预测原理及实例

灰色预测理论
陈文广
青岛理工大学管理学院
灰色系统的定义
白色系统是指一个系统的内部特征是完全已知的，即系统的信息是完全充分的。
黑色系统是指一个系统的内部信息对外界来说是一无所知的，只能通过它与外界的联系来加以观测研究。
灰色系统是指“部分信息已知,部分信息未知”的“小样本”,“ 贫信息”的不确定性系统,它通过对“部分”已知信息的生成、开发去了解、认识现实世界，实现对系统运行行为和演化规律的正确把握和描述.
k
k1
x(r1) (i) x(r1) (i) x(r1) (k)
i1
i1
(2) (x(r) (k)) x(r2) (k)
(i) (x(r) (k)) x(ri) (k)
灰色系统的定义和特点
常用的灰色预测有五种：
（1）数列预测，即用观察到的反映预测对象特征的时间序列来构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。
（2）灾变与异常值预测，即通过灰色模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。
（3）季节灾变与异常值预测，即通过灰色模型预测灾变值发生在一年内某个特定的时区或季节的灾变预测。
(i) (x(r) (k )) i1 (x(r) (k )) i1 (x (r) (k 1)) 式中， 0为0次累减,即无累加,从而有关系式(3 5), 1为1次累减 (i)为i次累减
从而可得下述关系
(1) (x(r) (k)) 0(x(r) (k))0(x(r)C(k 1)) x(r) (k) x(r) (k 1)
数列累加
【例1】设原始数据序列
x ( 0 ) { x ( 0 ) ( 1 ) x ( 0 , ) ( 2 ) ,x ( 0 ) ( N ) } { 6 , 3 , 8 , 1 , 7 } 0

灰色预测法(GM(1-1)模型)

X 3 3.4,3.3,3.5,3.5
商业
X 4 6.7,6.8,5.4,4.7
参考序列分别为 X1, X 2 ，被比较序列为 X 3, X 4,
试求关联度。
回总目录回本章目录
. #;
解答：
以 X1 为参考序列求关联度。
第一步：初始化，即将该序列所有数据分别除以第一个数据。得到：
X1 1,0.9475,0.9235,0.9138
回总目录回本章目录
. #;
• 系统预测
通过对系统行为特征指标建立一组相互关联的灰色预测模型，预测系统中众多变量间的相互协调关系的变化。
• 拓扑预测
将原始数据做曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点。
回总目录回本章目录
dt
其中：α称为发展灰数；μ称为内生控制灰数。
回总目录回本章目录
. #;
设 ˆ 为待估参数向量，ˆ
a
，可利用
最小二乘法求解。解得：
ˆ BT B 1 BT Yn
求解微分方程，即可得预测模型：
Xˆ
1 k
1
X
0 1
a
e ak
a
k 0,1,2..., n
回总目录回本章目录
. #;
二、模型检验灰色预测检验一般有残差检验、关联度检
X 2 1,1.063,1.1227,1.1483
X
3
1,.097,1.0294,1.0294
X 4 1,1.0149,0.805,0.7
回总目录回本章目录
. #;
第二步：求序列差 2 0,0.1155,0.1992,0.2335

python实现灰色预测GM（1,1）模型灰色系统预测灰色预测公式推导

python实现灰⾊预测GM（1,1）模型灰⾊系统预测灰⾊预测公式推导来源公式推导连接关键词：灰⾊预测 python 实现灰⾊预测 GM(1,1)模型灰⾊系统预测灰⾊预测公式推导⼀、前⾔本⽂的⽬的是⽤Python和类对灰⾊预测进⾏封装⼆、原理简述1.灰⾊预测概述灰⾊预测是⽤灰⾊模型GM(1,1)来进⾏定量分析的，通常分为以下⼏类： (1) 灰⾊时间序列预测。

⽤等时距观测到的反映预测对象特征的⼀系列数量（如产量、销量、⼈⼝数量、存款数量、利率等）构造灰⾊预测模型，预测未来某⼀时刻的特征量，或者达到某特征量的时间。

(2) 畸变预测（灾变预测）。

通过模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。

(3) 波形预测，或称为拓扑预测，它是通过灰⾊模型预测事物未来变动的轨迹。

(4) 系统预测，对系统⾏为特征指标建⽴⼀族相互关联的灰⾊预测理论模型，在预测系统整体变化的同时，预测系统各个环节的变化。

上述灰⾊预测⽅法的共同特点是：（1）允许少数据预测；（2）允许对灰因果律事件进⾏预测，例如：灰因⽩果律事件：在粮⾷⽣产预测中，影响粮⾷⽣产的因⼦很多，多到⽆法枚举，故为灰因，然⽽粮⾷产量却是具体的，故为⽩果。

粮⾷预测即为灰因⽩果律事件预测。

⽩因灰果律事件：在开发项⽬前景预测时，开发项⽬的投⼊是具体的，为⽩因，⽽项⽬的效益暂时不很清楚，为灰果。

项⽬前景预测即为灰因⽩果律事件预测。

（3）具有可检验性，包括：建模可⾏性的级⽐检验（事前检验），建模精度检验（模型检验），预测的滚动检验（预测检验）。

2.GM(1,1)模型理论 GM(1,1)模型适合具有较强的指数规律的数列，只能描述单调的变化过程。

已知元素序列数据：做⼀次累加⽣成（1-AGO）序列：其中，令为的紧邻均值⽣成序列：其中，建⽴GM(1,1)的灰微分⽅程模型为：其中，为发展系数，为灰⾊作⽤量。

设为待估参数向量，即，则灰微分⽅程的最⼩⼆乘估计参数列满⾜其中再建⽴灰⾊微分⽅程的⽩化⽅程（也叫影⼦⽅程）：⽩化⽅程的解（也叫时间响应函数）为那么相应的GM(1,1)灰⾊微分⽅程的时间响应序列为：取，则再做累减还原可得即为预测⽅程。

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

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

>> %普通的灰色预测GM1.mclearX=input('请输入原始数据:','s');%原始数据（可以多行，每一行为一类原始数据，即可多类原始数据）X=str2num(X);[m1 m2]=size(X);%m1和m2分别表示X的行数和列数k0=input('请输入所要预测的阶数:');%GM(1,1)模型for i=1:m1n=i;x0=X(i,:);%将原始数据X中的第i行数据赋给x0，即取出一类原始数据disp('1.原始数据:');Y='';for z=1:m2Y=strcat(Y,'(',num2str(x0(z)),')');enddisp(Y);% 1. 利用一次累加（1-AGO）生成新数列E=triu(ones(m2));%E表示元素为1右上三角阵x1=x0*E;%对原始数据进行一次累加（1-AGO）生成新数列x1disp('2.一次累加（1-AGO）生成的数据:');Y='';for z=1:m2Y=strcat(Y,'(',num2str(x1(z)),')');enddisp(Y);% 2. 计算出发展系数a，灰作用量ub1=x1;b1(1)=[];b2=x1;b2(m2)=[];b=-0.5*(b1+b2);%生成B的第一行B=[b;ones(1,m2-1)];%生成BB=B'; %对B进行转置y0=x0;y0(1)=[];%生成Yy0=y0'; %对Y进行转置A=((inv(B'*B))*B')*y0; %根据Y=BA(A为待辨识参数向量)，可知A=[inv（B'B）]B'Ya=A(1);u=A(2);%A=[a u]'% 3. 确立模型且求出模拟值u_a=u/a;for k=0:m2+k0-1x2(k+1)=(x0(1)-u_a)*exp(-k*a)+u_a;end %求出新数列的模拟值x2x3=x2;x3(m2+k0)=[];x4=[0 x3];x5=x2-x4;%利用累减生成法求出原始数据的模拟值x5disp('3.一次累加（1-AGO）生成的数据的模拟值:')Y='';for z=1:m2+k0Y=strcat(Y,'(',num2str(x2(z)),')');enddisp(Y);disp('4.原始数据的预测值:')Y='';for z=1:m2+k0Y=strcat(Y,'(',num2str(x5(z)),')');enddisp(Y);% 4. 模型检验（算出的值到等级参照表中检查其精度等级）%计算后验差比C和残差序列Qx6=x5(1:m2);Q=x0-x6;%Q为残差序列s1=std(Q);%s1为残差序列Q的标准差s2=std(x0);%s2为初始序列x0的标准差C=s1/s2;%后验差比C(越小越好)w1=1:m2;w1=[ones(m2,1) w1'];w2=Q';[bb,bint,r1,rint,stats]=regress(w2,w1);rcoplot(r1,rint)C1=strcat('5.后验差比(均方差比值): C=',num2str(C));disp(C1);if C<=0.35disp(' 由于C<=0.35，则此模型精度等级为1级（好）。

No.15-第6章-灰色预测

xˆu(1) (7) fu (6 1) x(1) (6) 1 min 33.0071 xˆu(1) (8) fu (6 2) x(1) (6) 2 min 36.6595 xˆu(1) (9) fu (6 3) x(1) (6) 3 min 40.3119

(4.9445,5.5828,5.3441,5.2669, 4.5640,3.6524)
试计算其一次累加生成序列 x(1) (7), x(1) (8), x(1) (9) 的最高预测值、最低预测值和基本预测值。
max

max{x(0) (k)}
1n{x(0)(k)}
四、区间预测
最高预测值最低预测值基本预测值
xˆs(1) (7) fs (6 1) x(1) (6) 1 max 34.9375 xˆs(1) (8) fs (6 2) x(1) (6) 2 max 40.5203 xˆs(1) (9) fs (6 3) x(1) (6) 3 max 46.1031
X (0) =47343.2,53641.3,57072.6,59746.9,66581.6
(1)级比检验
σ(k) x(k -1) = σ(2), σ(3), σ(4), σ(5)
x(k)
= 0.8826, 0.9399, 0.9552, 0.8973
n5
级比覆盖区域为
2
2
(e n1, en1 ) (0.716531,1.395612)
级比序列都在其覆盖范围内，可以建立GM(1,1)模型。
【案例6-1】城镇基本医疗参保人数预测
(2)构造累加生成序列 X(1) 47343.2,100984.5,158057.1,217804,284385.6 (3)构造数据矩阵 B 和数据向量 Y 计算 αˆ

灰色预测例题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因此模拟序列为求模拟序列和原始序列的相关度初始化，即将该序列所有数据分别除以第一个数据。

灰色预测一些程序

clc,clearsyms a b;c=[a b]';A=[174 179 183 189 207 234 220.5 256 270 285];%注意这里为列向量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)F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;endG=[];G(1)=A(1);for i=2:(n+10)G(i)=F(i)-F(i-1);endt1=1995:2004;t2=1995:2014;G,a,b %输出预测值、发展系数和灰色作用量plot(t1,A,'o',t2,G)function []=greymodel(y)% 本程序主要用来计算根据灰色理论建立的模型的预测值。

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

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

y=input('请输入数据 ');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;i=1:n+2;yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1);for j=n+2:-1:2ys(j)=yys(j)-yys(j-1);endx=1:n;xs=2:n+2;yn=ys(2:n+2);plot(x,y,'^r',xs,yn,'*-b');det=0;sum1=0;sumpe=0;for i=1:nsumpe=sumpe+y(i);endpe=sumpe/n;for i=1:n;sum1=sum1+(y(i)-pe).^2;ends1=sqrt(sum1/n);sumce=0;for i=2:nsumce=sumce+(y(i)-yn(i)); endce=sumce/(n-1);sum2=0;for i=2:n;sum2=sum2+(y(i)-yn(i)-ce).^2;ends2=sqrt(sum2/(n-1));c=(s2)/(s1);disp(['后验差比值为：',num2str(c)]);if c<0.35disp('系统预测精度好')else if c<0.5disp('系统预测精度合格')else if c<0.65disp('系统预测精度勉强')elsedisp('系统预测精度不合格')endendenddisp(['下个拟合值为 ',num2str(ys(n+1))]); disp(['再下个拟合值为',num2str(ys(n+2))]);。

灰色预测实例

第一题kN m k b pk N m L g f mgp S )()(1∑--=--+--=当m<=N 时f mgp S -=当m>N 时kN m k b pk N m L g f mgp S )()(1∑--=--+--=现在设旅客达到机场概率为p=90%，N=300，f=0.6Ng ，g L b 5.0= 现在km k pk m gg mg S )300(*5.1180*9.0301∑-=----=取m=301 经过计算得到 S=（90.9-2.53*10^(-14)）*g 取m=302经过计算得到 S=（91.8-8.095*10^(-13)）*g取m=307经过计算得到S=（96.3-4.065*10^(-8)）*g取m=311经过计算得到S=（99.9-9.865*10^(-6)）*g取m=318经过计算得到S=（106.2-5.68*10^(-3)）*g取m=325经过计算得到S=（112.5-2.59*10^(-1)）*g取m=332经过计算得到S=（118.8-2.42）*g=116.38*g取m=336经过计算得到S=（122.4-5.42）*g=116.98g取m=337经过计算得到S=（123.3-6.38）*g=116.92g所以航空公司在出售336张票的时候收益最大值为116.98g，由于这只是单方面考虑到肮空公司的利润，在实际中，国内超售可以达到5%，国外一般是2%。

对于拒载的赔偿问题，早已有法律规定是按照里程数进行赔偿，程序 m=337; x=0.9*m-180 y=0;for k=0:1:(m-301)y=y+(m-300-k)*nchoosek(m,k)*0.1^(k)*0.9^(m-k); end 1.5*y第二题首先假设购买打折票的旅客与全票的旅客不到概率是一样的都为pa 为购买打折票未到的人数，b 为购买全票未到的人数，k 为未到达的人数，k=a+b 。

灰色预测

回归分析是应用最广泛的一种办法。

但回归分析要求大样本，只有通过大量的数据才能得到量化的规律，这对很多无法得到或一时缺乏数据的实际问题的解决带来困难。

回归分析还要求样本有较好的分布规律，而很多实际情形并非如此。

例如，我国建国以来经济方面有几次大起大落，难以满足样本有较规律的分布要求。

因此，有了大量的数据也不一定能得到统计规律，甚至即使得到了统计规律，也并非任何情况都可以分析。

另外，回归分析不能分析因素间动态的关联程度，即使是静态，其精度也不高，且常常出现反常现象灰色系统理论提出了一种新的分析方法—关联度分析方法，即根据因素之间发展态势的相似或相异程度来衡量因素间关联的程度，它揭示了事物动态关联的特征与程度。

由于以发展态势为立足点，因此对样本量的多少没有过分的要求，也不需要典型的分布规律，计算量少到甚至可用手算，且不致出现关联度的量化结果与定性分析不一致的情况。

这种方法已应用到农业经济、水利、宏观经济等各方面，都取得了较好的效果。

灰色系统理论建模的主要任务是根据具体灰色系统的行为特征数据，充分开发并利用不多的数据中的显信息和隐信息，寻找因素间或因素本身的数学关系。

通常的办法是采用离散模型，建立一个按时间作逐段分析的模型。

但是，离散模型只能对客观系统的发展做短期分析，适应不了从现在起做较长远的分析、规划、决策的要求。

尽管连续系统的离散近似模型对许多工程应用来讲是有用的，但在某些研究领域中，人们却常常希望使用微分方程模型。

事实上，微分方程的系统描述了我们所希望辨识的系统内部的物理或化学过程的本质。

灰色系统理论首先基于对客观系统的新的认识。

尽管某些系统的信息不够充分，但作为系统必然是有特定功能和有序的，只是其内在规律并未充分外露。

有些随机量、无规则的干扰成分以及杂乱无章的数据列，从灰色系统的观点看，并不认为是不可捉摸的。

相反地，灰色系统理论将随机量看作是在一定范围内变化的灰色量，按适当的办法将原始数据进行处理，将灰色数变换为生成数，从生成数进而得到规律性较强的生成函数。

灰色预测GM（1，1）方法

灰色预测法一、相关知识1、灰色预测通过原始数据的处理和灰色模型的建立，发现、掌握系统发展规律，对系统的未来状态做出科学的定量预测。

2、灰数简介：（1）灰数的定义：是指未明确指定的数，即处在某一范围内的数，灰数是区间数的一种推广。

灰数实际上指在某一个区间或某个一般的数集内取值的不确定数，通常用记号“⊗”表示灰数。

（2）灰数的分类：(Ⅰ)有下界而无上界的灰数[)∞∈⊗,a 或()a ⊗，如大树的重量必大于零，但不可能用一般手段知道其准确的重量,所以其重量为灰数[)∞∈⊗,0。

(Ⅱ)有上界而无下界的灰数(,]a ⊗∈-∞或()a ⊗，如一项投资工程，要有个最高投资限额，一件电器设备要有个承受电压或通过电流的最高临界值。

(Ⅲ)既有下界a 又有上界a 的灰数称为区间灰数，记为[]a a ,∈⊗。

如海豹的重量在20--25公斤之间，某人的身高在1.8-1.9米之间，可分别记为[]25,201∈⊗，[]9.1,8.12∈⊗(Ⅳ)黑数：当()∞∞-∈⊗,或()21,⊗⊗∈⊗，即当⊗的上、下界皆为无穷或上、下界都是灰数时，称⊗为黑数。

(Ⅴ)白数：当[,]a a ⊗∈且a a =时，称⊗为白数。

(3)本征灰数是指不能或暂时还不能找到一个白数作为其“代表”的灰数，比如一般的事前预测值、宇宙的总能量、准确到秒或微妙的“年龄”等都是本征灰数。

非本征灰数是指凭先验信息或某种手段，可以找到一个白数作为其“代表”的灰数。

我们称此白数为相应灰数的白化值，记为⊗~，并用()a ⊗表示以a 为白化值的灰数。

如托人代买一件价格100元左右的衣服，可将100作为预购衣服价格()100⊗的白化数，记为()100100~=⊗。

例：（1）气温不超过36℃，[]36,0∈⊗。

（2）预计某地区今年夏粮产量在100万吨以上，[)∞∈⊗,100；（3）估计某储蓄所年底居民存款总额将达7000万到9000万，[]9000,7000∈⊗；（4）如某人希望至少获得1万元科研经费，并且越多越好，[)∞∈⊗,10000；（5）有的数，从系统的高层次，即宏观层次、整体层次或认识的概括层次上看是白的，可到低层次上，即到系统的微观层次、分部层次或认识的深化层次则可能是灰的。