关于R语言实现GM(1,N)模型中出现的问题

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

看到很少有用R语言作GM(1,N)模型的，自己根据GM(1,1)的代码修改了一下，出现了一些问题。GM(1,2)的时候跑《灰色系统理论及其应用》里例题的数据是和书上吻合的，但是拓展到GM(1,3)的情况时，预测结果误差非常大，我觉得不应该是这样，但是不知道是为什么？有大神可以帮忙看看我的代码哪里出现问题了吗？十分感谢！

附件的excel中sheet2是本次数据。

1. #GM(1,2)
   
2. 
   
3. gm12<-function(x10,x20,t){
   
4.   x11 <- cumsum(x10)
   
5.   x21 <- cumsum(x20)
   
6.   b <- numeric(length(x10)-1)
   
7.   n <- length(x10)-1
   
8.   for(i in 1:n){
   
9.     b[i]<--(x11[i]+x11[i+1])/2
   
10.     b}
    
11.   B <- cbind(b,x21[-1])
    
12.   BT <- t(B)
    
13.   M <- solve(BT%*%B)
    
14.   YN <- numeric(length(x10)-1)
    
15.   YN <- x10[2:length(x10)]
    
16.   alpha <- M%*%BT%*%YN  
    
17.   alpha2 <- matrix(alpha,ncol=1)
    
18.   a <- alpha2[1]
    
19.   u1 <- alpha2[2]
    
20.   cat("GM(1,N)参数估计值：",'\n',"系统发展系数-a=",-a,"  ",
    
21.       "驱动系数u1=",u1,'\n','\n')
    
22.   y <- numeric(length(c(1:t)))
    
23.   y[1]<-x10[1]
    
24.   for(w in 1:(t-1)){  
    
25.     y[w+1] <- (x10[1]-u1/a*x21[w+1])*exp(-a*w)+u1/a*x21[w+1]
    
26.   }
    
27.   cat("x(11)模拟值",'\n',y,'\n')
    
28.   xy <- numeric(length(y))
    
29.   xy[1] <- y[1]
    
30.   for(o in 2:t){
    
31.     xy[o] <- y[o] - y[o-1]
    
32.   }
    
33.   cat("x(10)模拟值",'\n',xy,'\n','\n')                       
    
34.   
    
35.   #计算残差e
    
36.   e <- numeric(length(x10))
    
37.   for(l in 1:length(x10)){
    
38.     e[l] <- x10[l]-xy[l]
    
39.   }
    
40.   cat("残差：",'\n',e,'\n')
    
41.   #计算相对误差
    
42.   e2 <- numeric(length(x10))
    
43.   for(s in 1:length(x10)){
    
44.     e2[s] <- (abs(e[s])/x10[s]) #μÃÏà¶ÔÎó2î
    
45.   }
    
46.   cat("相对残差",'\n',e2,'\n','\n')
    
47.   cat("残差平方和=",sum(e^2),'\n')
    
48.   cat("平均相对误差=",sum(e2)/(length(e2)-1)*100,"%",'\n')
    
49.   cat("相对精度=",(1-(sum(e2)/(length(e2)-1)))*100,"%",'\n','\n')
    
50.   
    
51.   #后验差比值检验
    
52.   avge<-mean(abs(e));esum<-sum((abs(e)-avge)^2);evar=esum/(length(e)-1);se=sqrt(evar)  
    
53.   avgx10<-mean(x10);x10sum<-sum((x10-avgx10)^2);x10var=x10sum/(length(x10));sx=sqrt(x10var)  
    
54.   cv<-se/sx  
    
55.   cat("后验差比值检验:",'\n',"C值=",cv,'\n')
    
56.   if(cv < 0.35){     
    
57.     cat("C值<0.35, GM(1,N)预测精度等级为：好",'\n','\n')
    
58.   }else{
    
59.     if(cv<0.5){
    
60.       cat("C值属于[0.35,0.5), GM(1,N)模型预测精度等级为：合格",'\n','\n')
    
61.     }else{
    
62.       if(cv<0.65){
    
63.         cat("C值属于[0.5,0.65), GM(1,N)模型预测精度等级为：勉强合格",'\n','\n')
    
64.       }else{
    
65.         cat("C值>=0.65, GM(1,N)模型预测精度等级为：不合格",'\n','\n')
    
66.       }
    
67.     }
    
68.   }
    
69.   #画出输入序列x10的预测序列及x10的比较图像
    
70.   plot(xy,col='blue',type='b',pch=16,xlab='时间/年',ylab='全国人均每次国内旅游消费/元')
    
71.   points(x10,col='red',type='b',pch=4)
    
72.   legend('topleft',c('预测值','原始值'),pch=c(16,4),lty=l,col=c('blue','red'))
    
73. }
    
74. 
    
75. 
    
76. a<-c(2.874,3.278,3.307,3.39,3.679)
    
77. b<-c(7.04,7.645,8.075,8.53,8.774)
    
78. gm12(x20,x30,length(x20))

复制代码

1. #GM(1,3)

复制代码

*可以看到误差非常大

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

1. gm13<-function(x10,x20,x30,t){ #x10为系统特征数据序列，xi0为相关因素序列，t为预测个数
   
2.   x11 <- cumsum(x10) #一次累加生成序列1-AG0序列
   
3.   x21 <- cumsum(x20)
   
4.   x31 <- cumsum(x30)
   
5.   b <- numeric(length(x10)-1)
   
6.   n <- length(x10)-1
   
7.   for(i in 1:n){ #生成x10的紧邻均值生成序列
   
8.     b[i]<--(x11[i]+x11[i+1])/2
   
9.     b} #得序列b，即为x10的紧邻均值生成序列
   
10.   B <- cbind(b,x21[-1],x31[-1])
    
11.   BT <- t(B)#做逆矩阵
    
12.   M <- solve(BT%*%B)
    
13.   YN <- numeric(length(x10)-1)
    
14.   YN <- x10[2:length(x10)]
    
15.   alpha <- M%*%BT%*%YN  #模型的最小二乘估计参数列满足alpha尖
    
16.   alpha2 <- matrix(alpha,ncol=1)
    
17.   a <- alpha2[1]
    
18.   u1 <- alpha2[2]
    
19.   u2 <- alpha2[3]
    
20.   cat("GM(1,N)参数估计值：",'\n',"系统发展系数-a=",-a,"  ",
    
21.       "驱动系数u1=",u1,'\n','\n',
    
22.       "驱动系数u2=",u2,'\n','\n')#利用最小二乘法求得参数估计值a,u1,u2
    
23.   y <- numeric(length(c(1:t)))
    
24.   y[1]<-x10[1]
    
25.   for(w in 1:(t-1)){  #将a,u1,u2的估计值代入时间响应序列函数计算x10拟合序列y
    
26.     y[w+1] <- (x10[1]-u1/a*x21[w+1]-u2/a*x31[w+1])*exp(-a*w)+u1/a*x21[w+1]+u2/a*x31[w+1]
    
27.   }
    
28.   cat("x(11)的模拟值：",'\n',y,'\n')
    
29.   xy <- numeric(length(y))
    
30.   xy[1] <- y[1]
    
31.   for(o in 2:t){ #运用后减运算还原得模型输入序列10预测序列
    
32.     xy[o] <- y[o] - y[o-1]
    
33.   }
    
34.   cat("x(10)的模拟值：",'\n',xy,'\n','\n')                       
    
35.   
    
36.   #计算残差e
    
37.   e <- numeric(length(x10))
    
38.   for(l in 1:length(x10)){
    
39.     e[l] <- x10[l]-xy[l] #得残差
    
40.   }
    
41.   cat("残差：",'\n',e,'\n')
    
42.   #计算相对误差
    
43.   e2 <- numeric(length(x10))
    
44.   for(s in 1:length(x10)){
    
45.     e2[s] <- (abs(e[s])/x10[s]) #得相对误差
    
46.   }
    
47.   cat("相对残差：",'\n',e2,'\n','\n')
    
48.   cat("残差平方和=",sum(e^2),'\n')
    
49.   cat("平均相对误差=",sum(e2)/(length(e2)-1)*100,"%",'\n')
    
50.   cat("相对精度=",(1-(sum(e2)/(length(e2)-1)))*100,"%",'\n','\n')
    
51.   
    
52.   #后验差比值检验
    
53.   avge<-mean(abs(e));esum<-sum((abs(e)-avge)^2);evar=esum/(length(e)-1);se=sqrt(evar)  #计算残差的方差se
    
54.   avgx10<-mean(x10);x10sum<-sum((x10-avgx10)^2);x10var=x10sum/(length(x10));sx=sqrt(x10var)  #计算原序列x10的方差sx
    
55.   cv<-se/sx  #得验差比值
    
56.   cat("后验差比值检验:",'\n',"C值=",cv,'\n')#对后验差比值进行检验，与一般标准进行比较判断预测结果好坏。
    
57.   if(cv < 0.35){     
    
58.     cat("C值<0.35, GM(1,N)预测精度等级为：好",'\n','\n')
    
59.   }else{
    
60.     if(cv<0.5){
    
61.       cat("C值属于[0.35,0.5), GM(1,N)模型预测精度等级为：合格",'\n','\n')
    
62.     }else{
    
63.       if(cv<0.65){
    
64.         cat("C值属于[0.5,0.65), GM(1,N)模型预测精度等级为：勉强合格",'\n','\n')
    
65.       }else{
    
66.         cat("C值>=0.65, GM(1,N)模型预测精度等级为：不合格",'\n','\n')
    
67.       }
    
68.     }
    
69.   }
    
70.   #画出输入序列x10的预测序列及x10的比较图像
    
71.   plot(xy,col='blue',type='b',pch=16,xlab='时间/年',ylab='全国人均每次国内旅游消费/元')
    
72.   points(x10,col='red',type='b',pch=4)
    
73.   legend('topleft',c('预测值','原始值'),pch=c(16,4),lty=l,col=c('blue','red'))
    
74. }
    
75. 
    
76. 
    
77. 
    
78. getwd()
    
79. 
    
80. #缺失值处理
    
81. data=read.xlsx("data1.xls",2)
    
82. data$X2=as.numeric(as.character(data$X2))
    
83. sub=which(is.na(data$X2))
    
84. print (sub)
    
85. print(data[1,])
    
86. value=data$X1[1]-(data$X2[3]-data$X2[2])
    
87. data$X2[sub]=value
    
88. 
    
89. x10 = data$Y..
    
90. x20 = data$X1
    
91. x30 = data$X2
    
92. 
    
93. gm13(x10,x20,x30,length(x10))

复制代码

2L是GM(1,3)的代码，不知道为什么1L没显示

xc793997241 发表于 2019-10-16 14:33
2L是GM(1,3)的代码，不知道为什么1L没显示

请问误差较大的原因找到了吗