GM（1，n）的R语言代码，有会编的给一个！！

Python016

GM（1，n）的R语言代码，有会编的给一个！！,第1张

#灰色预测模型GM(1,1)#用法：#假设数列1 2 3 4 5.5 6 7.5 为已知数据，你要预测后面3项，gm11([1 2 3 4 5.5 6 7.5],10) # 10=7+3# 序列输入格式为：x<-c(1,2,3,4,5.5,6,7.5)gm11<-function(x,k){#x为行向量数据#做一次累加n<-length(x)x1<-numeric(n)for(i in 1:n){x1[i]<-sum(x[1:i])}#x1的均值数列z1<-numeric(n)m<-n-1for(j in 1:m){z1[j+1]<-(0.5*x1[j+1]+0.5*x1[j])}Yn=t(t(x[2:n]))B<-matrix(1,nrow=n-1,ncol=2)B[,1]<-t(t(-z1[2:n]))#solve(M)求M的逆#最小二乘法求解参数列u<-solve(t(B)%*%B)%*%t(B)%*%Yna<-u[1]b<-u[2]#预测x2<-numeric(k)x2[1]<-x[1]for(i in 1:k-1){x2[1+i]=(x[1]-b/a)*exp(-a*i)+b/a}x2=c(0,x2)#还原数据y=diff(x2)y}#调用函数x<-c(1,2,3,4,5.5,6,7.5)gm11(x,10)

clearclc

X0=[128453 129227 130765 129988 131448 132129 132802 133450 134091 135404]

pre_num=5  %预测5年

%% 级比检验

n=length(X0)

Xle=exp(-2/(n+1))

Xre=exp(2/(n+1))

lambda=X0(1:end-1)./X0(2:end)

range=minmax(lambda)

if range(2)<Xre && range(1)>Xle

    disp('所有的级比都落在可容覆盖区间，可以建立GM模型')

else

    disp('没有通过级比检验')

end

%% 建模GM（1,1）

X1=cumsum(X0)

Z1=0.5*(X1(2:end)+X1(1:end-1))

Y=X0(2:end)'

B=[-Z1(1:end)' ones(n-1,1)]

u=B\Y %u=inv(B'*B)*B'*Y

a=u(1)

b=u(2)

%% 输出结果

Xpre=[X0(1) ones(1,n-1+pre_num)]

for k=1:n-1+pre_num

    Xpre(k+1)=(X0(1)-b/a)*(exp(-a*k)-exp(-a*(k-1)))

end

err=X0-Xpre(1:n)   %计算残差

epsilon=abs(err)./X0*100 %计算相对误差

rho=1-(1-0.5*a)/(1+0.5*a)*lambda %计算级比偏差值

%% 画图

t1=2002:2011

t2=2002:2011+pre_num

plot(t1,X0,'o',t2,Xpre,'r','LineWidth',2)

xlabel('年份')

ylabel('人口')

legend('原始数据','预测数据', 'Location','SouthEast')

预测的话，应该用接下来的时间，所以应该是预测2014,2015....

程序如下：

new<-data.frame(year=2014)

lm.pred<-predict(z,new,interval="prediction",level=0.95)

lm.pred

解释：第一行表示输入新的点year=2014，注意，即使就一个点，也要采用数据框结构；第二行的函数predict（）给出相应的预测值，参数interval="prediction"表示同时要给出相应的置信区间，参数level=0.95表示相应的概率为0.95.这个参数也可以不写，因为它的缺省值为0.95.

你提到的2013的数据不是预测，而是拟合。我们可以通过得到的模型对原来的year这个变量的数据进行拟合。

程序如下：

fit<-fitted(z)

fit

得到的就是在你得到的模型下2006-2013这8年的拟合值了。

希望能对你有所帮助~