灰色预测

理论提出者
　　首先是由华中科技大学的邓聚龙教授提出的理论。理论概念
　　通过少量的、不完全的信息，建立灰色微分预测模型，对事物发展规律作出模糊性的长期描述（模糊预测领域中理论、方法较为完善的预测学分支）。关联度
　　提出系统的关联度分析方法，是对系统发展态势的量化比较分析。关联度的一般表达式为：
　　n
　　ri=1/n∑ xi(k)
　　i=1
　　ri 是曲线xi对参考曲线x0的关联度。生成数
　　通过对原始数据的整理寻找数的规律，分为三类：
　　a、累加生成：通过数列间各时刻数据的依个累加得到新的数据与数列。累加前数列为原始数列，累加后为生成数列。
　　b、累减生成：前后两个数据之差，累加生成的逆运算。累减生成可将累加生成还原成非生成数列。
　　c、映射生成：累加、累减以外的生成方式。基本关系式
　　记x(0)为原始数列
　　x(0)=( x(0)(k)xk=1,2,…,n)=(x(0)(1),x(0)(2),…,x(0)(n))
　　记x(1)为生成数列
　　x(1)=( x(1)(k)xk=1,2,…,n)=(x(1)(1),x(1)(2),…,x(1)(n))
　　如果x(0) 与x(1)之间满足下列关系，即
　　k
　　x(1)(k)= ∑x(0)(i)
　　i=a
　　称为一次累加生成。建立模型
　　a、建模机理
　　b、 把原始数据加工成生成数；
　　c、 对残差(模型计算值与实际值之差)修订后，建立差分微分方程模型；
　　d、 基于关联度收敛的分析；
　　e、 gm模型所得数据须经过逆生成还原后才能用。
　　f、采用“五步建模（系统定性分析、因素分析、初步量化、动态量化、优化）”法，建立一种差分微分方程模型gm(1,1)预测模型。基本算式为：
　　令 x(0)=(x(0)(1),x(0)(2),…,x(0)(n))
　　作一次累加生成， k
　　x(1)(k)= ∑x(0)(m)
　　m=1
　　有 x(1)=(x(1)(1),x(1)(2),…,x(1)(n))
　　=(x(0)(1),x(1)(1)+x(0)(2),…,x(1)(n-1)+x(0)(n))
　　x(1)可建立白化方程：dx(1)/dt+ax(1)=u 即gm(1,1).
　　该方程的解为: x(1)(k+1)=(x(1)(1)-u/a)e-ak+u/a预测方法：
　　a、数列预测
　　b、灾变预测
　　c、季节灾变预测
　　d、拓扑预测
　　e、系统综合预测
　　f、模糊预测特点
　　对于一个模糊系统来说，传统的预测方法就会失去作用。处理模糊预测问题的数学方法是模糊数学。模糊数学的基础是模糊集合论，而模糊集合是普通集合的扩展。美国学者l.a.zadeh教授建立的模糊集合论，为模糊预测理论与方法的研究奠定了理论基础。它用简捷有力的方法处理复杂系统，在某种程度上弥补了经典数学与统计数学的不足。应用
　　在预测应用上，如气象预报、地震预报、病虫害预报等，国内学者做出了许多有益的研究。基于混沌理论的分析预测
　　混沌理论是近年来长足发展的一门学科。混沌向世界规律运动的假定性提出挑战。一方面，它告诉我们，宇宙远比我们想得到的要怪异，它使许多传统的科学方法受到怀疑。另一方面，混沌认为许多无规则的事物实际上可能是简单规律的结果。混沌展现给我们的是一些新的规律。混沌理论
　　遵从简单规律的系统会以令人惊讶的复杂方式表现其行为。混沌是隐秘形式的秩序。
　　混沌系统是指敏感地依赖于初始条件的内在变化系统，对外来变化的敏感性本身并不意味着混沌。混沌理论最令人兴奋的是：一个非常简单的决定论系统能够产生异常复杂的输出结果。给定一个简单规则和初始条件，系统将产生复杂连续系列，这一点类似“无中生有”。研究应用实例
　　美国科学家帕卡德和他的同事基于混沌和生物进化理论，借助计算机，致力于用图形来描述金融市场的混沌现象。帕卡德认为，世界上有大量不同的随机现象，他所研究的是大体只需几个变量就能描述系统行为的一种混沌现象。他试图建立一种学习算法，对进化模型进行处理。而对于众多的模型，帕卡德采用一种称为遗传算法的方法处理数据。它用类似生物繁殖中突变和杂交现象的方法来改变模型。这种方法的核心是，计算机不断设定新的假设环境，从而使学习算法更具有适应性。认为一个好的学习算法不仅能建立适应模型，它还能时刻观测数据的变化。所谓“学习算法”是一种特别的程序，他擅长对大量的、各种各样的模型进行比较研究，找出哪个模型最适用于分析目前和未来的数据。matlab程序
　　%下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,matlab6.5 ,使用本程序请注明，程序存储为gm1.m
　　%x = [5999,5903,5848,5700,7884];gm1(x); 测试数据
　　%二次拟合预测GM(1,1)模型
　　function gmcal=gm1(x)
　　ynyz = [
　　124 113 154.1666667 240.0909091 224.3333333
　　180.9 155.0909091 223.4166667 274.8333333 327.25
　　33.33333333 92.5 57.5 382.5 101.3333333
　　60 69.75 62.66666667 136.6666667 202.75
　　48.66666667 60.33333333 49.75 97.25 139.6666667
　　111.4 98.18181818 118.75 168.8 243.5454545
　　280.25 296.125 316.7777778 354 458.8888889
　　62.66666667 61.33333333 95.75 39.66666667 79
　　223.1 205 235.3636364 303.8 355
　　148.25 162 161.1666667 243 361.6
　　650 522.8 629.7142857 915.6666667 1270.4
　　780.3333333 726.6666667 924 1483.666667 1656.666667
　　140 104.75 168.2 322 337.3333333
　　292.6666667 344.3333333 454.3333333 416.3333333 943.5
　　168.3333333 161.5 339 415.6666667 262.8
　　290.3333333 270.25 317 548 494
　　313.25 543 565.6666667 696.75 713
　　308 328.4285714 379 432.1428571 521.8333333
　　1723.25 2060 2107 3643.666667 2416.6
　　285.75 508 364.25 526.25 1142
　　386.3333333 511.2 558.1666667 549.5 674.5
　　348.2857143 412.3461538 502.96 694.88 747.56
　　4091.914286 4137.055556 4823.685714 6175.722222 7706.305556
　　520 619 547.5 786 1142
　　643.2352941 489.8095238 699.6666667 1133.25 1298.533333
　　964.25 972 1220.809524 1408.565217 2286.909091
　　210.6 234.8 304 428.5 544.8571429
　　79.75 110.6666667 107.2 283 415
　　265.7333333 286.5714286 400.8125 415 425.6875
　　1873.285714 1623 1854.222222 2042.9 2665.5
　　219.175 311.7 406.8928571 480.75 642.0769231
　　88.2 171 233 178.8 438
　　33.33333333 183 223 269 733
　　105.4 122.8125 130.5625 265.5 235.9444444
　　134.4 104.3333333 247.25 293 385.1666667
　　134.8 147 216.8571429 186.6666667 246.25
　　33.33333333 44.25 116.75 96 54.6
　　124.2 147 253 544 530.75
　　114.5 109.5 150.5 219.3333333 300.8571429
　　57.66666667 51.6 89.6 166.6666667 253
　　959.75 1120.25 2188 3439.5 4430
　　2461.857143 3007.714286 4046.571429 4783 6518.5
　　1871.2 1785.833333 3962.25 2489 7272
　　2396.833333 2322.714286 2950.142857 3408 4164.25
　　5154.25 5092 8304.333333 19573 13888.33333
　　2737 2675.428571 3227 3299.75 6860.4
　　1500.75 1813.857143 2180.166667 2722 3583.857143
　　1417 1205.8 1414.5 1211.777778 1276.9
　　255 278.6363636 305.1666667 386.4615385 596.0769231
　　20 33.33333333 100 132 217.5
　　894 614.3333333 545 425.3333333 563
　　285.8 308.4615385 443.75 511.2 710.1
　　299.6 316 443 893.3333333 1023.666667
　　240 278.5714286 361.6666667 704 848.6
　　50 61.33333333 100 273 841
　　119.2 192.25 251 445 715
　　159.6 142.4 225.5 396.3333333 357.5
　　33.33333333 63.66666667 71 177.5 137
　　256.6666667 257 375.6666667 664.5 677
　　215.1666667 194 329.6666667 361.1428571 385
　　296.6666667 265 244.25 249 240.6666667
　　292.5 355.75 415.75 407.2 574.4
　　317.5 306.2 408.5 389.2 543.6
　　315 409.4 504.8 472.6666667 711.4
　　84.33333333 116.5 153 128 253
　　632 1013 855.3333333 1089.333333 1361.333333
　　227.25 275.75 348 361 520
　　225.3333333 260.4 264.5 305 426.5
　　237.5 351.8 378 375.6666667 550.8333333
　　77 106.6 104.3333333 166.8 346.8
　　148 443 189.3333333 379.5 604.5
　　474.5 1182 1250 1019 1263.5
　　];
　　x = ynyz(72,:);
　　sizexd2 = size(x,2);
　　%求数组长度
　　k=0;
　　for y1=x
　　k=k+1;
　　if k>1
　　x1(k)=x1(k-1)+x(k);
　　%累加生成
　　z1(k-1)=-0.5*(x1(k)+x1(k-1));
　　%z1维数减1，用于计算B
　　yn1(k-1)=x(k);
　　else
　　x1(k)=x(k);
　　end
　　end
　　%x1,z1,k,yn1
　　sizez1=size(z1,2);
　　%size(yn1);
　　z2 = z1';
　　z3 = ones(1,sizez1)';
　　YN = yn1'; %转置
　　%YN
　　B=[z2 z3];
　　au0=inv(B'*B)*B'*YN;
　　au = au0';
　　%B,au0,au
　　afor = 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,rightbra
　　strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)
　　%输出时间响应方程
　　%******************************************************
　　%二次拟合
　　k2 = 0;
　　for y2 = x1
　　k2 = k2 + 1;
　　if k2 > k
　　else
　　ze1(k2) = exp(-(k2-1)*afor);
　　end
　　end
　　%ze1
　　sizeze1 = size(ze1,2);
　　z4 = ones(1,sizeze1)';
　　G=[ze1' z4];
　　X1 = x1';
　　au20=inv(G'*G)*G'*X1;
　　au2 = au20';
　　%z4,X1,G,au20
　　Aval = au2(1);
　　Bval = au2(2);
　　%Aval,Bval
　　%输出预测的 A,B的值
　　strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(Bval),rightbra)
　　%输出时间响应方程
　　nfinal = sizexd2-1 + 1;
　　%决定预测的步骤数5 这个步骤可以通过函数传入
　　%nfinal = sizexd2 - 1 + 1;
　　%预测的步骤数 1
　　for k3=1:nfinal
　　x3fcast(k3) = constant1*exp(afor1*k3)+ua;
　　end
　　%x3fcast
　　%一次拟合累加值
　　for k31=nfinal:-1:0
　　if k31>1
　　x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);
　　else
　　if k31>0
　　x31fcast(k31+1) = x3fcast(k31)-x(1);
　　else
　　x31fcast(k31+1) = x(1);
　　end
　　end
　　end
　　x31fcast
　　%一次拟合预测值
　　for k4=1:nfinal
　　x4fcast(k4) = Aval*exp(afor1*k4)+Bval;
　　end
　　%x4fcast
　　for k41=nfinal:-1:0
　　if k41>1
　　x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);
　　else
　　if k41>0
　　x41fcast(k41+1) = x4fcast(k41)-x(1);
　　else
　　x41fcast(k41+1) = x(1);
　　end
　　end
　　end
　　x41fcast,x
　　%二次拟合预测值
　　%***精度检验p C************//////////////////////////////////
　　k5 = 0;
　　for y5 = x
　　k5 = k5 + 1;
　　if k5 > sizexd2
　　else
　　err1(k5) = x(k5) - x41fcast(k5);
　　end
　　end
　　%err1
　　%绝对误差
　　xavg = mean(x);
　　%xavg
　　%x平均值
　　err1avg = mean(err1);
　　%err1avg
　　%err1平均值
　　k5 = 0;
　　s1total = 0 ;
　　for y5 = x
　　k5 = k5 + 1;
　　if k5 > sizexd2
　　else
　　s1total = s1total + (x(k5) - xavg)^2;
　　end
　　end
　　s1suqare = s1total ./ sizexd2;
　　s1sqrt = sqrt(s1suqare);
　　%s1suqare,s1sqrt
　　%s1suqare 残差数列x的方差 s1sqrt 为x方差的平方根S1
　　k5 = 0;
　　s2total = 0 ;
　　for y5 = x
　　k5 = k5 + 1;
　　if k5 > sizexd2
　　else
　　s2total = s2total + (err1(k5) - err1avg)^2;
　　end
　　end
　　s2suqare = s2total ./ sizexd2;
　　%s2suqare 残差数列err1的方差S2
　　Cval = sqrt(s2suqare ./ s1suqare);
　　Cval
　　%nnn = 0.6745 * s1sqrt
　　%Cval C检验值
　　k5 = 0;
　　pnum = 0 ;
　　for y5 = x
　　k5 = k5 + 1;
　　if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrt
　　pnum = pnum + 1;
　　%ppp = abs( err1(k5) - err1avg )
　　else
　　end
　　end
　　pval = pnum ./ sizexd2;
　　pval
　　%p检验值
　　%arr1 = x41fcast(1:6)
　　%预测结果为区间范围 预测步长和数据长度可调整程序参数进行改进