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

Python039

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-1

for(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)%*%Yn

a<-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)

这里写了四个函数，方便在Matlab里面调用，分别是GM（1，1），残差GM（1，1），新陈代谢GM（1，1），Verhust自己写得难免有所疏忽，需要的朋友自己找本书本来试验一下。。

Gm（1，1）

function [px0,ab,rel]=gm11(x0,number)

%[px0,ab,rel]=gm11(x0,number)

%px0为预测数列，rel为平均相对误差,rel为平均相对误差(为百分比)

%默认的number参数为原数组大小

if nargin==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换

number=max(size(x0))

end

n=max(size(x0))%取输入数据的样本量

x1=zeros(size(x0))

for k=1:n

for i=1:k

x1(k)=x1(k)+x0(i)%计算累加值，并将值赋予矩阵be

end

end

z=zeros(size(x0))

for k=2:n

z(k)=0.5*(x1(k)+x1(k-1))%计算数据矩阵B的第一列数据

end

y=x0'

y(1)=[]

b(:,1)=-z'

b(:,2)=1

b(1,:)=[]

ab=inv(b'*b)*b'*y %计算参数 矩阵

a=ab(1)

b=ab(2)

px0(1)=x0(1)

%求还原值系列

for k=1:number-1

px0(k+1)=(1-exp(a)) * ( x0(1)-b/a ) * exp(-a*k)

end

temp=px0(1:n)

x0

temp=(temp-x0)./x0 %相对误差

temp(1)=[]%删除第一个为零的误差

temp=abs(temp)

rel=sum(temp)/(n-1)*100

残差Gm（1，1）

function [px0,ab,rel]=ccgm11(x0,number)

%[px0,ab,rel]=gm11(x0,number)

%px0为残差预测数列，ab为求得的系数,rel为平均相对误差(为百分比)

%默认的number参数为原数组大小

if nargin==1

number=max(size(x0))

end

n=max(size(x0)) %数组大小..

[px0,ab,rel]=gm11(x0,number)

wucha=x0-px0(1:n)

i=n

%求后面的同号的数目.

while(wucha(i)*wucha(i-1)>0 &i>=2)

i=i-1

end

start=i

length=n-i+1

new=wucha(start:n)

if length>=4

pwucha=gm11(new)

px0(start:n)=px0(start:n)+pwucha

clear wucha

wucha=px0-x0

wucha=wucha./x0 %相对误差

wucha=abs(wucha)

rel=sum(wucha)/(n-1)*100

end

verhust

function [px0,ab,rel]=verhust(x1,number)

%[px0,ab,rel]=verhust(x0,number)

%px0为预测数列，rel为平均相对误差,rel为平均相对误差(为百分比)

%默认的number参数为原数组大小

if nargin==1

number=max(size(x1))

end

n=max(size(x1))

x0(1)=x1(1)

for k=2:n

x0(k)=x1(k)-x1(k-1)

z(k)=0.5*(x1(k)+x1(k-1))

end

x0

z

B=[-(z(2:n))' (z(2:n).^2)']

B

Y=(x0(2:n))'

Y

ab=inv(B'*B)*B'*Y

a=ab(1)b=ab(2)

for k=1:number

px0(k)=(a*x1(1))/(b*x1(1)+(a-b*x1(1)).*exp(a*(k-1)))

end

temp=px0(1:n)

x1

temp=(temp-x1)./x1 %相对误差

temp(1)=[]%删除第一个为零的误差

temp=abs(temp)

rel=sum(temp)/(n-1)*100

新陈代谢Gm（1，1）

function [px0,ab,rel]=xcdxgm11(x0,number,step)

%[px0,ab,rel]=xcdxgm11(x0,number,step)

%x0为原系列,number为要预测的数目,step为基本步长

%px0为预测数列，rel为平均相对误差,rel为平均相对误差(为百分比)

%默认的number参数为原数组大小

%模型假设预测的数据和原始数据都要大于等于5

if nargin==1

number=max(size(x0))

step=max(size(x0))

end

if nargin==2

step=max(size(x0))

end

n=max(size(x0))

if n<step | n<5

error('此模型要求至少有五个原始数据,并且原始数据个数要大于新陈代谢的步长.')

end

[px0,ab,rel]=gm11(x0,n)

last=n

x0

px0

while last<number

begin=last-step+1

temp=px0(begin:last)

temp=gm11(temp,step+1)

last=last+1

px0(last)=temp(step+1)

end