灰色模型的malab程序

function [y,p,e]=huise_1_1(X,k)                     %灰色模型的malab程序
%Example [y,p]=gm_1_1([200 250 300 350],2) 
%接口描述： X的预测的初始数列，|X|>4，K是指向后进行预测的个数
%命令格式： 程序保存的文件名，eg：huise.m  则命令是： huise([579.8 547.5 527.0 492.3 437.0],5) 
 if nargout>3;
    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_orig=X;
predict=k;                                %AGO 处理，即是对初始数列进行一阶累加
 x=cumsum(x_orig);                         %计算系数(a 和 u)------------------------
 n=length(x_orig);                           %生成矩阵 B
 
for i=1:(n-1);
  B(i)=-(x(i)+x(i+1))/2;
 end
 B=[B' ones(n-1,1)];                            %生成矩阵 Y
 
for i=1:(n-1);
  y(i)=x_orig(i+1);
 end
 Y=y';                                     %计算系数 a=au(1) u=au(2)
 
au=(inv(B'*B))*(B'*Y);            %--------------------------------------------------------
                                             %把huise模型公式转换成符号
 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(costr1,'+',costr2,'e^',costr3,'*(t-1))');            %计算每一个值
 
for t=1:(n+predict)
 mcv(t)=coef1+coef2*exp(coef3*(t-1));
 end
 x_mcv0=diff(mcv);
 x_mcve=[x_orig(1) x_mcv0]                         %输出图形中的各点
 x_mcv=diff(mcv(1:end-predict));
 x_orig_n=x_orig(2:end);
 x_c_error=x_orig_n-x_mcv;
 x_error=mean(abs(x_c_error./x_orig_n));
 
if x_error>0.2                                      %相对误差的均值
disp('model disqualification!');
 elseif x_error>0.1
 disp('model check out');
 else
 disp('model is perfect!');
 end                                    
 plot(1:n,x_orig,'o',1:n+predict,x_mcve);
 p=x_mcve(end-predict+1:end);                      %画出预测模型和初始数列的点
 xlabel('年份(从第一个数据年份起)');
ylabel('产水量（万吨）');
 title('灰度模型    GM(1,1)');
 grid on
 y=eq;
 e=x_error;
 p=x_mcve(end-predict+1:end);