R实现灰色预测

原文链接: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为预测个数
  x1<-cumsum(x0) #一次累加生成序列1-AG0序列(累加序列)
  b<-numeric(length(x0)-1)# 初始化长度为length(x0)-1的整数部分个numeric类型且值为0的数据集b
  n<-length(x0)-1# n 为length(x0)-1长度 因为需要生成MEAN(紧邻均值)生成序列 其长度少1
  for(i in 1:n){ #生成x1的紧邻均值生成序列
    b[i]<--(x1[i]+x1[i+1])/2 
    b} #得序列b,即为x1的紧邻均值生成序列
  D<-numeric(length(x0)-1)
  D[]<-1
  B<-cbind(b,D)#作B矩阵
  BT<-t(B)#B矩阵转置
  M<-solve(BT%*%B)#求BT*B得逆
  YN<-numeric(length(x0)-1)
  YN<-x0[2:length(x0)]
  alpha<-M%*%BT%*%YN  #模型的最小二乘估计参数列满足alpha尖
  alpha2<-matrix(alpha,ncol=1)# 将结果变成一列
  # 得到方程的两个系数
  a<-alpha2[1]
  u<-alpha2[2]
  
  # 下面为结果输出
  # 输出参数估计值及模拟值
  cat("GM(1,1)参数估计值:",'\n',"发展系数-a=",-a,"  ","灰色作用量u=",u,'\n','\n') #利用最小二乘法求得参数估计值a,u
  y<-numeric(length(c(1:t)))# t为给定的预测个数
  y[1]<-x1[1] # 第一个数不变
  for(w in 1:(t-1)){  #将a,u的估计值代入时间响应序列函数计算x1拟合序列y
    y[w+1]<-(x1[1]-u/a)*exp(-a*w)+u/a 
  }
  cat("x(1)的模拟值:",'\n',y,'\n')
  xy<-numeric(length(y))
  xy[1]<-y[1]
  for(o in 2:t){ #运用后减运算还原得模型输入序列x0预测序列
    xy[o]<-y[o]-y[o-1] 
  } 
  cat("x(0)的模拟值:",'\n',xy,'\n','\n')                       
  
  # 计算残差e
  e<-numeric(length(x0))
  for(l in 1:length(x0)){
    e[l]<-x0[l]-xy[l] #得残差序列(未取绝对值)
  }
  cat("绝对残差:",'\n',e,'\n')
  #计算相对误差
  e2<-numeric(length(x0))
  for(s in 1:length(x0)){
    e2[s]<-(abs(e[s])/x0[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')
  
  #后验差比值检验
  avge<-mean(abs(e));esum<-sum((abs(e)-avge)^2);evar=esum/(length(e)-1);se=sqrt(evar)  #计算残差的均方差se
  avgx0<-mean(x0);x0sum<-sum((x0-avgx0)^2);x0var=x0sum/(length(x0));sx=sqrt(x0var)  #计算原序列x0的方差sx
  cv<-se/sx  #得验差比值(方差比)
  cat("后验差比值检验:",'\n',"C值=",cv,'\n')#对后验差比值进行检验,与一般标准进行比较判断预测结果好坏。
  #计算小残差概率
  P<-sum((abs(e)-avge)<0.6745*sx)/length(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)预测精度等级为:好

R实现灰色预测

 

分享到:更多 ()