多变量灰色理论系统GM(1,n)

clear all
%输入待预测时刻k及原始序列X0
k=1;
X0=xlsread('mark','sheet1','a:b');
%对原始序列X0累加生成序列X1
[n,m]=size(X0);
for j=1:m c=0;
for i=1:n
c=X0(i,j)+c;
X1(i,j)=c;
end
end
%计算数据矩阵L
for j=1:m
for i=1:n-1
l(i,j)=(X1(i,j)+X1(i+1,j))/2;
end
end
L=[l ones(n-1,1)];
%计算Y及参数估计值
for j=1:m
Y(1:n-1,j)=X0(2:n,j);
a(:,j)=inv(L'*L)*L'*Y(1:n-1,j);
end
a=a';
A=a(1:end,1:end-1);
B=a(1:end,end);
%计算模型的拟合值或预测值
S=X1(1,1:end);
if k==1
Z=S'
elseif k>1
Z=expm(A*(k-1))*S'+inv(A)*(expm(A*(k-1))-eye(size(expm(A*(k-1)))))*B-(expm(A*(k-2))*S'+inv(A)*(expm(A*(k-2))-eye(size(expm(A*(k-2)))))*B)
else disp('输入错误!k不得小于1')
end