分享

R实现灰色预测

 火树银花5 2017-12-04


原文链接:http://www.cnblogs.com/homewch/p/5783073.html


1.简介


预测就是借助于对过去的探讨去推测、了解未来。灰色预测通过原始数据的处理和灰色模型的建立,发现、掌握系统发展规律,对系统的未来状态做出科学的定量预测。对于一个具体的问题,究竟选择什么样的预测模型应以充分的定性分析结论为依据。模型的选择不是一成不变的。一个模型要经过多种检验才能判定其是否合适,是否合格。只有通过检验的模型才能用来进行预测。本章将简要介绍灰数、灰色预测的概念,灰色预测模型的构造、检验、应用,最后对灾变预测的原理作了介绍。

灰色系统理论的产生和发展动态

  • 1982邓聚龙发表第一篇中文论文《灰色控制系统》标志着灰色系统这一学科诞生。

  • 1985灰色系统研究会成立,灰色系统相关研究发展迅速。

  • 1989海洋出版社出版英文版《灰色系统论文集》,同年,英文版国际刊物《灰色系统》杂志正式创刊。目前,国际、国内200多种期刊发表灰色系统论文,许多国际会议把灰色系统列为讨论专题。国际著名检索已检索我国学者的灰色系统论著500多次。灰色系统理论已应用范围已拓展到工业、农业、社会、经济、能源、地质、石油等众多科学领域,成功地解决了生产、生活和科学研究中的大量实际问题,取得了显著成果。


灰色系统与模糊数学、黑箱方法的区别

  • 主要在于对系统内涵与外延处理态度不同;研究对象内涵与外延的性质不同。

  • 灰色系统着重外延明确、内涵不明确的对象,模糊数学着重外延不明确、内涵明确的对象。

  • “黑箱”方法着重系统外部行为数据的处理方法,是因果关系的两户方法,使扬外延而弃内涵的处理方法,而灰色系统方法是外延内涵均注重的方法。


2.GM(1,1)的R语言实现


2.1 R源码

#灰色预测模型G(1,1)GM11<-function(x0,t){ #x0为输入学列,t为预测个数=""><-cumsum(x0) #一次累加生成序列1-ag0序列(累加序列)=""><-numeric(length(x0)-1)# 初始化长度为length(x0)-1的整数部分个numeric类型且值为0的数据集b=""><-length(x0)-1# n="" 为length(x0)-1长度="" 因为需要生成mean(紧邻均值)生成序列="" 其长度少1=""  for(i="" in="" 1:n){="" #生成x1的紧邻均值生成序列=""  =""><--(x1[i]+x1[i+1])  =""  b}="" #得序列b,即为x1的紧邻均值生成序列=""><-numeric(length(x0)-1)><-1><-cbind(b,d)#作b矩阵><-t(b)#b矩阵转置><-solve(bt%*%b)#求bt*b得逆><-numeric(length(x0)-1)><-x0[2:length(x0)]><-m%*%bt%*%yn  #模型的最小二乘估计参数列满足alpha尖=""><-matrix(alpha,ncol=1)# 将结果变成一列=""  #="" 得到方程的两个系数=""><-alpha2[1]><-alpha2[2]  =""  #="" 下面为结果输出=""  #="" 输出参数估计值及模拟值=""  cat('gm(1,1)参数估计值:','\n','发展系数-a=',-a,'  ','灰色作用量u=',u,' \n','\n')="" #利用最小二乘法求得参数估计值a,u=""><-numeric(length(c(1:t)))# t为给定的预测个数=""><-x1[1] #="" 第一个数不变=""  for(w="" in="" 1:(t-1)){=""  #将a,u的估计值代入时间响应序列函数计算x1拟合序列y=""  =""><-(x1[1]-u )*exp(-a*w)+u/a=""  }=""  cat('x(1)的模拟值:','\n',y,'\n')=""><-numeric(length(y))><-y[1]  for(o="" in="" 2:t){="" #运用后减运算还原得模型输入序列x0预测序列=""  =""><-y[o]-y[o-1]  }=""  cat('x(0)的模拟值:','\n',xy,'\n','\n')=""  =""  =""  =""  =""  =""  =""  =""  =""  =""  =""  =""  =""  #="" 计算残差e=""><-numeric(length(x0))  for(l="" in="" 1:length(x0)){=""  =""><-x0[l]-xy[l] #得残差序列(未取绝对值)=""  }=""  cat('绝对残差:','\n',e,'\n')=""  #计算相对误差=""><-numeric(length(x0))  for(s="" in="" 1:length(x0)){=""  =""><-(abs(e[s]) 0[s])="" #得相对误差=""  }=""  cat('相对残差:','\n',e2,'\n','\n')=""  cat('残差平方和=',sum(e^2),' \n')=""  cat('平均相对误差=',sum(e2)/(length(e2)-1)*100,' %','\n')=""  cat('相对精度=',(1-(sum(e2)/(length(e2)-1)))*100,' %','\n','\n')=""  =""  #后验差比值检验=""><><-sum((abs(e)-avge)^2);evar=esum length(e)-1);se="sqrt(evar)"  #计算残差的均方差se=""><><-sum((x0-avgx0)^2);x0var=x0sum length(x0));sx="sqrt(x0var)"  #计算原序列x0的方差sx=""><-se x=""  #得验差比值(方差比)=""  cat('后验差比值检验:','\n','c值=',cv,' \n')#对后验差比值进行检验,与一般标准进行比较判断预测结果好坏。=""  #计算小残差概率=""><><0.6745*sx) ength(e)=""  cat('小残差概率:','\n','p值=',P,' \n')=""  =""  if(cv="">< 0.35="" &&="" p="">0.95){        cat('C<0.35, p="">0.95,GM(1,1)预测精度等级为:好','\n','\n')  }else{    if(cv<0.5 &&="" p="">0.80){      cat('C值属于[0.35,0.5), P>0.80,GM(1,1)模型预测精度等级为:合格','\n','\n')    }else{      if(cv<0.65 &&="" p="">0.70){        cat('C值属于[0.5,0.65), P>0.70,GM(1,1)模型预测精度等级为:勉强合格','\n','\n')      }else{        cat('C值>=0.65, GM(1,1)模型预测精度等级为:不合格','\n','\n')      }    }  }  #画出输入序列x0的预测序列及x0的比较图像  plot(xy,col='blue',type='b',pch=16,xlab='时间序列',ylab='值')  points(x0,col='red',type='b',pch=4)  legend('topleft',c('预测','原始'),title = '预测时序与原始时序对比',pch=c(16,4),lty=l,col=c('blue','red'))}

2.2 GM11 测试

> x<-c(2.67,3.13,3.25,3.36,3.56,3.72)> GM11(x,length(x)+2)GM(1,1)参数估计值: 发展系数-a= 0.04396098    灰色作用量u= 2.925617 x(1)的模拟值: 2.67 5.78087 9.031547 12.42832 15.97774 19.68668 23.56231 27.61211 x(0)的模拟值: 2.67 3.11087 3.250677 3.396768 3.549424 3.708941 3.875626 4.049803 绝对残差: 0 0.01913011 -0.0006772995 -0.03676788 0.010576 0.01105927 相对残差: 0 0.006111858 0.0002083999 0.01094282 0.002970786 0.002972921 残差平方和= 0.001952456 平均相对误差= 0.4641357 % 相对精度= 99.53586 % 后验差比值检验: C值= 0.0407599 小残差概率: P值= 1 C<0.35, p="">0.95,GM(1,1)预测精度等级为:好



    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多