グレー予測GM(1,1)モデル実装IML&R
関連記事ではGM(1,1)モデルについて言及しているが,劉思峰の『灰色システム』の第6章の詳細はないような気がする.次にマトリクス言語アルゴリズムを用いて実現し,コードは校正後,中国衛生統計における2つの文章につながった.
/* 2006 */
data a1;
INPUT t year xt@@;/* */
yt+xt;/* */
index=1;
zt=-(yt+LAG(yt))/2;/* B */
DATAlINES;
1 1990 24395 2 1991 25286
3 1992 26901 4 1993 27339
5 1994 27871 6 1995 28721
7 1996 29728 8 1997 30067
9 1998 30791 10 1999 31284
11 2000 33716 12 2001 34558
;
PROC IML;/* IML */
USE a1;/* SAS a1*/
READ ALL VAR{zt index}INTO B WHERE(zt ^= .);/* a1 zt index ( B*/
READ ALL VAR{xt}INTO Yn WHERE(zt ^= .);/* a1 xt Yn*/
ahat=INV(B`*B)*B`*Yn;/* */
ahatt=ahat`; /* */
na={a u};
CREATE a2 FROM ahatt[COLNAME=na];/* SAS a2*/
APPEND FROM ahatt;/* */
print b yn ahat ahatt;
QUIT;/* IML */
DATA a3;
SET a2;
index=1;
RUN;
DATA a4;
SET a1;
IF _N_ =1;
xt0=xt;
KEEP xt0 index;
RUN;
DATA a5;
MERGE a1 a3 a4;
BY index;
IF _N_ =1 THEN xp=xt;
ELSE DO
yt1=(xt0-u/a)*EXP(-a*(t-1))+u/a;/* */
yt0=(xt0-u/a)*EXP(-a*(t-2))+u/a;/* */
xp=yt1-yt0;/* */
END;
error=xp-xt;/* */
rerror=error/xt*100;/* */
/* DROP yt index zt yt1 yt0 xt0;*/
PROC PRINT DATA=a5;/* */
RUN;/* */
/* 2008.12 */
/* ： data al jbi=lag (x0)／x0 ；
data al data pan */
data a1;
INPUT t year xt@@;/* */
yt+xt;/* */
index=1;
jbi=lag(xt)/xt;
zt=-(yt+LAG(yt))/2;/* B */
DATAlINES;
1 1990 24395 2 1991 25286
3 1992 26901 4 1993 27339
5 1994 27871 6 1995 28721
7 1996 29728 8 1997 30067
9 1998 30791 10 1999 31284
11 2000 33716 12 2001 34558
;
data pan;/* */
set a1;
if 0.1353< =jbi< =7.389 then good=1;
else good=.; /* */
title 'panduanmoxing';
proc print data=pan;/* */
/* ： ， */
proc means data=a5 std mean noprint;/* */
var xt error;
output out=a5_2 std=sl s2 mean=x_e_;
data a5_3;/* C */
set a5_2;
c=s2/sl;
if 0.65 < c then jdu=0;
else if 0.5<c< =0.65 then jdu=3;
else if 0.35<c<0.5 then jdu=2;
else jdu=1;
/*drop sl s2 _type_ _freq_;*/
title 'houyanchabi';
proc print data=a5_3;run;/* C */
/* ： ， ， */
data a6;/* */
input t year @@;
datalines;
14 2006 15 2007 16 2008 17 2009 18 2010
;
data a7;/* */
merge a3 a4;
array t(6)(12 13 14 15 16 17);/* */
do i=2 to 6;
x1k1=(xt0-u/a)* exp(-a*t(i))+u/a;
x1k0=(xt0-u/a)* exp(-a*t(i-1))+u/a;
xp=x1k1-x1k0;
output;
end;
/*drop tl t2 t3 t4 t5 t6 a b x01 i x1k1 x1k0 index;*/
data a8;
merge a6 a7;
title 'yuce';
proc print data=a8;/* */
run;/* */
# R GM(1,1)
gm11<-function(x0,t){ #x0 ，t
x1<-cumsum(x0) # 1-AG0
b<-numeric(length(x0)-1)
n<-length(x0)-1
for(i in 1:n){ # x1
b[i]<--(x1[i]+x1[i+1])/2
b} # b， x1
D<-numeric(length(x0)-1)
D[]<-1
B<-cbind(b,D)
BT<-t(B)#
M<-solve(BT%*%B)
YN<-numeric(length(x0)-1)
YN<-x0[2:length(x0)]
alpha<-M%*%BT%*%YN # alpha
alpha2<-matrix(alpha,ncol=1)
a<-alpha2[1]
u<-alpha2[2]
cat("GM(1,1) ：",'
'," -a=",-a," "," u=",u,'
','
') # a,u
y<-numeric(length(c(1: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("x(1) ：",'
',y,'
')
xy<-numeric(length(y))
xy[1]<-y[1]
for(o in 2:t){ # x0
xy[o]<-y[o]-y[o-1]
}
cat("x(0) ：",'
',xy,'
','
')
# e
e<-numeric(length(x0))
for(l in 1:length(x0)){
e[l]<-x0[l]-xy[l] #
}
cat(" ：",'
',e,'
')
#
e2<-numeric(length(x0))
for(s in 1:length(x0)){
e2[s]<-(abs(e[s])/x0[s]) #
}
cat(" ：",'
',e2,'
','
')
cat(" =",sum(e^2),'
')
cat(" =",sum(e2)/(length(e2)-1)*100,"%",'
')
cat(" =",(1-(sum(e2)/(length(e2)-1)))*100,"%",'
','
')
#
avge<-mean(abs(e));esum<-sum((abs(e)-avge)^2);evar=esum/(length(e)-1);se=sqrt(evar) # se
avgx0<-mean(x0);x0sum<-sum((x0-avgx0)^2);x0var=x0sum/(length(x0));sx=sqrt(x0var) # x0 sx
cv<-se/sx #
cat(" :",'
',"C =",cv,'
')# ， 。
if(cv < 0.35){
cat("C <0.35, GM(1,1) ： ",'
','
')
}else{
if(cv<0.5){
cat("C [0.35,0.5), GM(1,1) ： ",'
','
')
}else{
if(cv<0.65){
cat("C [0.5,0.65), GM(1,1) ： ",'
','
')
}else{
cat("C >=0.65, GM(1,1) ： ",'
','
')
}
}
}
# x0 x0
plot(xy,col='blue',type='b',pch=16,xlab=' ',ylab=' ')
points(x0,col='red',type='b',pch=4)
legend('topleft',c(' ',' '),pch=c(16,4),lty=l,col=c('blue','red'))
}
a<-c(1.95,2.23,2.4,2.15,1.8,1.95)
gm11(a,length(a)+6)