看到很少有用R语言作GM(1,N)模型的，自己根据GM(1,1)的代码修改了一下，出现了一些问题。GM(1,2)的时候跑《灰色系统理论及其应用》里例题的数据是和书上吻合的，但是拓展到GM(1,3)的情况时，预测结果误差非常大，我觉得不应该是这样，但是不知道是为什么？有大神可以帮忙看看我的代码哪里出现问题了吗？十分感谢！
附件的excel中sheet2是本次数据。
- #GM(1,2)
- gm12<-function(x10,x20,t){
- x11 <- cumsum(x10)
- x21 <- cumsum(x20)
- b <- numeric(length(x10)-1)
- n <- length(x10)-1
- for(i in 1:n){
- b[i]<--(x11[i]+x11[i+1])/2
- b}
- B <- cbind(b,x21[-1])
- BT <- t(B)
- M <- solve(BT%*%B)
- YN <- numeric(length(x10)-1)
- YN <- x10[2:length(x10)]
- alpha <- M%*%BT%*%YN
- alpha2 <- matrix(alpha,ncol=1)
- a <- alpha2[1]
- u1 <- alpha2[2]
- cat("GM(1,N)参数估计值：",'\n',"系统发展系数-a=",-a,"",
- "驱动系数u1=",u1,'\n','\n')
- y <- numeric(length(c(1:t)))
- y[1]<-x10[1]
- for(w in 1:(t-1)){
- y[w+1] <- (x10[1]-u1/a*x21[w+1])*exp(-a*w)+u1/a*x21[w+1]
- }
- cat("x(11)模拟值",'\n',y,'\n')
- xy <- numeric(length(y))
- xy[1] <- y[1]
- for(o in 2:t){
- xy[o] <- y[o] - y[o-1]
- }
- cat("x(10)模拟值",'\n',xy,'\n','\n')
- #计算残差e
- e <- numeric(length(x10))
- for(l in 1:length(x10)){
- e[l] <- x10[l]-xy[l]
- }
- cat("残差：",'\n',e,'\n')
- #计算相对误差
- e2 <- numeric(length(x10))
- for(s in 1:length(x10)){
- e2[s] <- (abs(e[s])/x10[s]) #μÃÏà¶ÔÎó2î
- }
- 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)
- avgx10<-mean(x10);x10sum<-sum((x10-avgx10)^2);x10var=x10sum/(length(x10));sx=sqrt(x10var)
- cv<-se/sx
- cat("后验差比值检验:",'\n',"C值=",cv,'\n')
- if(cv < 0.35){
- cat("C值<0.35, GM(1,N)预测精度等级为：好",'\n','\n')
- }else{
- if(cv<0.5){
- cat("C值属于[0.35,0.5), GM(1,N)模型预测精度等级为：合格",'\n','\n')
- }else{
- if(cv<0.65){
- cat("C值属于[0.5,0.65), GM(1,N)模型预测精度等级为：勉强合格",'\n','\n')
- }else{
- cat("C值>=0.65, GM(1,N)模型预测精度等级为：不合格",'\n','\n')
- }
- }
- }
- #画出输入序列x10的预测序列及x10的比较图像
- plot(xy,col='blue',type='b',pch=16,xlab='时间/年',ylab='全国人均每次国内旅游消费/元')
- points(x10,col='red',type='b',pch=4)
- legend('topleft',c('预测值','原始值'),pch=c(16,4),lty=l,col=c('blue','red'))
- }
- a<-c(2.874,3.278,3.307,3.39,3.679)
- b<-c(7.04,7.645,8.075,8.53,8.774)
- gm12(x20,x30,length(x20))
*可以看到误差非常大
|
熟悉论坛请点击新手指南
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
