灰色模型GM(1-1)&matlab实现

年份全国人口年份全国人口19789.6259199612.238919799.7542199712.362619809.8705199812.4761198110.0072199912.5786198210.1654200012.6743198310.3008200112.7627198410.4357200212.8453198510.5851200312.9227198610.7507200412.9988198710.93200513.0756198811.1026200613.1448198911.2704200713.2129199011.4333200813.2802199111.5823200913.345199211.7171201013.4091199311.8517201113.4735199411.985201213.5404199512.1121单位：亿请根据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),,()]Xxxxxn设(1)X为(0)X的一次累加序列：即：(1)(0)(1)(1)(1)(1)(0)()(-1)(),2,,xxxkxkxkkn得(1)(1)(1)(1)(1)[(1),(2),(3),,()]Xxxxxn利用(1)X计算GM(1,1)模型参数a、u，令ˆ[,]Taau,则有1ˆ()TTNaBBBY其中(1)(1)(1)(1)(1)(1)1((1)(2))121((2)(3))121((1)())12xxxxBxnxn(0)(0)(0)[(2),(3),,()]TNYxxxn由此获得：(1)(1)ˆ(1)((1))akuuxkxeaa于是：(0)(1)(1)ˆˆˆ(1)(1)()1xkxkxk（）或(0)(0)ˆ(1)-((1))akuxkaxea（2）注：（1）式是根据(0)X和(1)X的关系的到的（2）式是利用数学求导还原得到的至于用哪个，最好看相对误差5.2模型求解：求解得：六、模型检验6.1残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验。设模拟值的残差序列为(0)()ek，则(0)(0)(0)ˆ()()()ekxkxk令()k为残差相对值，即残差百分比为(0)(0)(0)ˆ()()()%()xkxkkxk令为平均残差，11()nttn.一般要求20%t，最好是10%t，符合要求.6.2后验差检验后验差检验，即对残差分布的统计特性进行检验.检验步骤如下：1、计算原始时间数列0(0)(0)(0)(1),(2),,()Xxxxn的均值和方差2(0)(0)2(0)11111(),()nnttxxtSxtxnn2、计算残差数列(0)(0)(0)(0)(1),(2),,()eeeen的均值e和方差22s2(0)2(0)21111(),()nntteetSetenn其中(0)(0)(0)ˆ()()(),1,2,,etxtxttn为残差数列.3、计算后验差比值21CSS4、计算小误差频率(0)1()0.6745PPeteS令0S=0.67451S，(0)()|()|tete，即0()PPtS.若对给定的00C，当0CC时，称模型为方差比合格模型；若对给定的00P，当0PP时，称模型为小残差概率合格模型.PC模型精度0.950.35优0.800.5合格0.700.65勉强合格0.700.65不合格表3后验差检验判别参照表6.3关联度检验关联度是用来定量描述各变化过程之间的差别.关联系数越大，说明预测值和实际值越接近.设(0)(0)(0)(0)ˆˆˆˆ()(1),(2),,()Xtxxxn(0)(0)(0)(0)()(1),(2),,()Xtxxxn序列关联系数定义为(0)(0)(0)(0)(0)(0)(0)(0)ˆˆmin()()max()(),0ˆˆ()()max()()1,0xtxtxtxtttxtxtxtxtt式中,(0)(0)ˆ()()xtxt为第t个点(0)x和(0)ˆx的绝对误差，()t为第t个数据的关联系数，称为分辨率，即取定的最大差百分比，0,一般取0.5.(0)()xt和(0)ˆ()xt的关联度为11ntrtn精度等级关联度均方差比值小误差概率好（1级）0.900.350.95合格（2级）0.800.500.80勉强（3级）0.700.650.70不合格（4级）0.700.650.70表2精度检验等级function[A,X0,X]=GM(x0)%GM(1,1)模型%x0原始数据，是单行或者单列数据%这里把x转成行向量处理%A是参数矩阵，第一个为参数a，第二个为参数u%X0矩阵有4列%第1列：原始数据%第2列：预测数据%第3列：残差%第4列：相对误差%X相对误差（第4列）的平均值X0=[];A=[];X=[];[m,n]=size(x0);ifm~=1&&n~=1disp('*************************');disp('数据有误');disp('*************************');elseifn==1x0=x0';endn=size(x0,2);x1(1)=x0(1);B=[];fori=2:nx1(i)=x1(i-1)+x0(i);%对原始数据进行一次累加，从而得到x1b=-0.5*(x1(i-1)+x1(i));b1=[b,1];B=[B;b1];%利用x1，得到B矩阵endYN=x0';YN(1)=[];%得到YN矩阵，注：x0是行向量A=(inv(B'*B))*(B'*YN);%inv(X)求方阵的逆，另外det(x)是求方阵的行列式%A矩阵有两个数，参数a和参数u%现在求利用求得的参数a、u进行预测（即把k值带入）X1(1)=x0(1);a=A(1);%参数au=A(2);%参数ufori=2:nX1(i)=(x0(1)-u/a)*exp(-1*a*(i-1))+u/a;%累加的预测数据endX0(1)=x0(1);fori=2:nX0(i)=X1(i)-X1(i-1);%预测，这是有X0与X1的关系求的%X0(i)=-1*a*(x0(1)-u/a)*exp(-1*a*(i-1));%这是有数学推导求的%至于那一个好，就看平均误差吧endX0=[x0'X0'];%实际值，预测值%残差检验fori=1:nX0(i,3)=abs(X0(i,2)-X0(i,1));%残差数列X0(i,4)=X0(i,3)/X0(i,1);%相对误差endX=mean(X0(:,4));%相对误差的平均值end