R语言#灰色预测模型

#灰色预测模型GM(1,1)
#数据使用文献数据：基于灰色模型的中国职业性耳鼻喉口腔疾病发病预测研究
#已知2011-2018年全国职业性耳鼻喉口腔疾病发病人数，使用GM(1,1)模型预测发病数
x0<-c(532,639,716,880,1097,1276,1608,1528)#原始发病数
#目前R语言中暂无灰色模型，需手工实现R语言的灰色模型函数
#第一步：建模可行性分析，级比检验
exp(-2/(8+1));exp(2/(8+1))#级比区间(0.80,1.25)
len=length(x0)
a<-as.numeric()#循环前需建立一个目标a空值
for (i in 2:len){
  a[i]<-x0[i-1] / x0[i]
  a}
a#序列所有级比落入(0.80,1.25)区间内,可进行建模

#第二步：GM(1,1)建模
gm<-function(x0,t){    #x0为输入序列,t为预测个数
  x1<-cumsum(x0) #（1）累加生成序列1-AGO序列(累加生成序列)
  b<-numeric(length (x0)-1)#length(x0)计算x0的长度，numeric(x0)生成指定长度的0向量，这里生成比x0长度少1的0向量，记作b
  n<-length(x0)-1    ## n 为length(x0)-1长度,因为需要生成MEAN（紧邻均值）序列,其长度少1
  for(i in 1:n){    #（2）循环语句i从1到n循环，生成x1的紧邻均值序列
    b[i]<--(x1[i] + x1[i+1])/2
    b}#得序列b,即为x1的紧邻均值生成序列
  #(3)建立灰微分方程与白化微分方程并求解
  D<-numeric(length(b)) # 向量D初始化
  D[]<-1#向量D的元素全部赋值1，即单位向量
  B<-cbind(b,D)    #cbind(,)以列方式将向量b和D合并成矩阵B
  BT<-t(B)    #t(),将矩阵B的转置
  M<-solve(BT%*%B)#求BT*B的逆(类似于实数的倒数)，这里M是BT*B的逆矩阵，%*%为矩阵乘法符号
  YN<-numeric(length(x0)-1)    # 向量 yN 初始化
  YN<-x0[2:length(x0)] #将原始向量x0除第一个外的其余元素赋与yN
  alpha<-M%*%BT%*%YN #最小二乘法计算微分方程参数,a和u赋予向量alpha
  alpha2<-matrix(alpha, ncol =1)#将结果变成1列
  #得到方程的两个系数
  a<-alpha2[1]    #提取参数a
  u<-alpha2[2]    #提取参数u
  # 输出参数估计值
  cat("GM(1,1)参数估计值：",'\n',"发展系数a=",a,"  ","灰色作用量u=",u,'\n','\n') #利用最小二乘法求得参数估计值a,u
  #（3）残差检验与后验差比值检验
  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("原始数列：",'\n',x0,'\n')
  cat("1-AGO的预测值：",'\n',y,'\n')
  #数据累减还原
  xy<-numeric(length(y))    #向量xy初始化
  xy[1]<-y[1] #向量xy的第1个与y的第1个值不变,初始化
  for(o in 2:t)#运用后减运算还原得模型输入序列x0预测序列
  {
    xy[o]<-y[o]-y[o-1]
  }    #，xy为还原值
  cat("预测值：",'\n',xy,'\n','\n')  
  #计算残差e
  xy<-round(xy,4)#将xy保留4位小数
  m<-length(x0)#m长度
  e<-numeric(length(x0)) #残差向量 e 初始化
  for(L in 1 :m){
    e[L]<-x0[L]-xy[L]
  }    #循环语句计算得残差序列e
  cat("残差：",'\n',e,'\n','\n')
  #计算相对误差
  e<-round(e,4)#e保留4位小数
  q<-numeric(length(x0))    #相对误差向量初始化
  for(L in 1:m){
    q[L]<-(abs(e[L])/x0[L]) #abs为取绝对值
  }    #循环语句计算相对误差向量q
  cat("相对误差：",'\n',q* 100,'\n','\n')
  cat("残差平方和=",sum(e^2),'\n')
  cat("平均相对误差=",sum(q)/(length(q)-1)* 100,"%",'\n')
  cat("相对精度=",(1-(sum(q)/(length(q)-1)))* 100,"%",'\n','\n')
  
  #后验差比值检验
  q<-round(q,4)#小数位
  se<-sd(e)    #计算残差向量e标准差
  sx<-sd(x0)    #计算原数列x0标准差
  C<-se/sx  #得后验差比值C（方差比）
  cat("后验差比值检验:",'\n',"C值=",C,'\n')#对后验差比值进行检验，与一般标准进行比较判断预测结果好坏。
  #计算小误差概率
  P<-sum(abs(e)-mean(abs(e))<0.6745*sx)/length(e)
  cat("小误差概率:",'\n',"P值=",P,'\n')
  
  if((P>=0.95)&(C<=0.35))cat("预测效果好")
  else if((P>=0.8)&(C<0.5))cat("预测效果合格")
  else if((P>=0.7)&(C<0.65))cat("预测效果勉强")
  else cat("预测效果不合格")
  #画出输入序列x0的预测序列及x0的比较图像
  plot(xy,col='blue',type='b',pch=16,xlab='时间序列',ylab='职业性耳鼻喉口腔疾病发病人数')
  points(x0,col='red',type='b',pch=18)
  legend("bottomright",c('预测值','原始值'),cex=0.8,pch=c(16,18),col=c('blue','red'))
}
gm(x0,8)
gm(x0,8+2)
gm(x0,8+5)