gm_1_1.m

灰色gm(1,1)模型

function [pred,info]=GM_1_1(X,k) 
%********************************************************
%灰色GM(1,1)模型
%Build the calculating dieplate for the typical gray model. 
%Designed by tong.   18 March,2009.   
%--------------------------------------------------------
%输入：
%X：输入值，需要预测的时间序列，n维向量，X=[x_0_(1),x_0_(2),...,x_0_(n)],要求每个值都大于零；
%k：X序列后的需要预测值的个数，k>=1；
%--------------------------------------------------------
%输出：
%pred=prediction，预测值向量，维数由k值确定；
%info:结构体，包含以下六个信息：
%info.err_comp_ave=error_comparative_average，平均相对误差值，数值；
%info.XX:预测还原后的X值，n-1维向量，不包含x_0_(1)对应的值,XX=[x_^0_(2),x_^0_(3),...,x_^0_(n)]；
%info.err_comp=error_comparative,相对误差值，n-1维向量，与XX的每个分量对应；
%info.err_abs=error_absolute,绝对误差值，n-1维向量，与XX的每个分量对应；
%info.equ=equation,白化方程的时间响应序列：x_1_(k)=(x_0_(1)-u/a)*exp(-a*(k-1))+u/a；
%info.au=au=[a u]'，参数a和u；
%--------------------------------------------------------
%Example: 刘思峰 谢乃明 等《灰色系统理论及其应用》第四版 第105页 例5.4.1
%X=[60.7 73.8 86.2 100.4 123.3];
%[pred,info]=GM_1_1(X,1);
%********************************************************


if nargout>2,error('Too many output argument.');end 
if nargin==1,k=1;x_orig=X; 
elseif nargin==0|nargin>2 
   error('Wrong number of input arguments.'); 
end 
%--------------------------------------------------------
%判断序列X是否为正数序列
if X>0;
else
    error('X must be larger than zero!');
end 
%--------------------------------------------------------
x_orig=X; 
predict=k; 
%--------------------------------------------------------
%1-AGO process 
x=cumsum(x_orig); 


%compute the coefficient(a and u)------------------------ 
n=length(x_orig); 
%first generate the matrix B 
for i=1:(n-1); 
B(i)=-(x(i)+x(i+1))/2; 
end 
B=[B' ones(n-1,1)]; 
%then generate the matrix Y 
for i=1:(n-1); 
y(i)=x_orig(i+1); 
end 
Y=y'; 
%get the coefficient. a=au(1) u=au(2) 
au=(inv(B'*B))*(B'*Y); %最小二乘估计系数，au=[a u]'
%-------------------------------------------------------- 
%change the grey model to symbolic expression 
%白化方程的时间响应序列：x_1_(k)=(x_0_(1)-u/a)*exp(-a*(k-1))+u/a
%eq=costr1+costr2*exp(costr3*(t-1))

coef1=au(2)/au(1); 
coef2=x_orig(1)-coef1; 
coef3=0-au(1); 
costr1=num2str(coef1); 
costr2=num2str(abs(coef2)); 
costr3=num2str(coef3); 
eq=strcat('x_1_(k)=(',costr2,')*e^(',costr3,'*(k-1))','+','(',costr1,')'); 

%-------------------------------------------------------- 

%comparison of calculated and observed value 

for t=1:n+predict 
   mcv(t)=coef1+coef2*exp(coef3*(t-1));  %mcv：时间响应序列得出的结果，还需要进行还原
end 
x_mcv0=diff(mcv); %x_mcv0：计算出来的还原值，包括预测值
x_mcve=[x_orig(1) x_mcv0]; %x_mcve：在最前面补上第一个原始值后的还原值，包含predict个预测值
x_mcv=diff(mcv(1:end-predict)); %x_mcv：计算出来的还原值，不包括第一个数据和后面预测值,x_1_(2),...,x_1_(n)
x_orig_n=x_orig(2:end); %x_orig_n:不包含第一个数据的原始数据，x(2),...,x(n)
x_a_error=x_orig_n-x_mcv; %x_a_error：绝对误差值(残差)，n-1维向量，不包含第一个数据的误差，因为第一个数据误差为0
x_c_error=abs(x_a_error./x_orig_n); %x_c_error：相对误差值,同上，没讨论第一个数据
x_error=mean(x_c_error); %平均相对误差，同上，没讨论第一个数据，也不应该讨论


disp('**************************');
if x_error>0.2 
   disp('model disqualification!'); 
elseif x_error>0.1 
   disp('model check out!'); 
else 
   disp('model is perfect!'); 
end 
disp('**************************');


%predicting model and plot gragh 
plot(1:n,x_orig,'b*',1:n+predict,x_mcve,'rs');
p=x_mcve(end-predict+1:end); %p：仅仅给出的是预测值的结果
xlabel('CURVE OF GREY MODEL ANALYSIS'); 
title('GM(1,1)'); 
grid on 

%-------------------------------------------------------- 
pred=p;

info.err_comp_ave=x_error; 
info.XX=x_mcv;
info.err_comp=x_c_error;
info.err_abs=x_a_error;
info.equ=eq; 
info.au=au;