灰色预测[GM(1_1)]MATLAB程序+改进+GM(1,m)多变量的+自己整理并运行过的

灰色预测[GM(1,1) ]MATLAB程序
% 本程序主要用来计算根据灰色理论建立的模型的预测值。

% 应用的数学模型是GM(1,1)。

% 原始数据的处理方法是一次累加法。

X0=input('请输入数据');%输入数据请用如例所示形式：[48.7 57.17 68.76 92.15]
n=length(X0);
X1=ones(n,1);
X1(1)=X0(1);
for i=2:n
X1(i)=X1(i-1)+X0(i);
end
B=ones(n-1,2);
for i=1:(n-1)
B(i,1)=-(X1(i)+X1(i+1))/2;
B(i,2)=1;
end
BT=B';
for j=1:n-1
Yn(j)=X0(j+1);
end
Yn=Yn';
A=inv(BT*B)*BT*Yn;
a=A(1);
u=A(2);
t=u/a;
t_test=input('请输入需要预测个数：');
i=1:t_test+n;
X1s(i+1)=(X0(1)-t).*exp(-a.*i)+t;
X1s(1)=X0(1);
for j=n+t_test:-1:2
X0s(j)=X1s(j)-X1s(j-1);
end
x=1:n;
xs=2:n+t_test;
X0ss=X0s(2:n+t_test);
plot(x,X0,'^r',xs,X0ss,'*-b');
det=0;
for i=2:n
det=det+abs(X0s(i)-X0(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
disp(['预测值为：',num2str(X0s(n+1:n+t_test))]);
改进程序1
（优化C值，c为微分方程的解的常数）见论文“刘发全,职承杰.灰色预测GM(1,1)模型的一点改
进”
% 本程序主要用来计算根据灰色理论建立的模型的预测值。

% 应用的数学模型是 GM(1,1)。

% 原始数据的处理方法是一次累加法。

clear
jm=xlsread('测试温度');
X0=jm;%input('请输入数据 ');%输入数据请用如例所示形式：[48.7 57.17 68.76 92.15] n=length(X0);
X1=ones(n,1);
X1(1)=X0(1);
for i=2:n
X1(i)=X1(i-1)+X0(i);
end
B=ones(n-1,2);
for i=1:(n-1)
B(i,1)=-(X1(i)+X1(i+1))/2;
B(i,2)=1;
end
BT=B';
for j=1:n-1
Yn(j)=X0(j+1);
end
Yn=Yn';
A=inv(BT*B)*BT*Yn;
a=A(1);
u=A(2);
for i=1:n-1
h1=0;
h1=h1+exp(-a.*i)/X0(i+1);
h2=0;
h2=h2+(exp(-a.*i)/X0(i+1))^2;
end
t=u/a;
h3=(X0(1)-t)*(1-exp(a)^(-1))/(X0(1)^2);
h4=((1-exp(a))^(-1))/(X0(1))^2;
%c=(h1+h3)/(h2+h4);
c=h1/h2;
t_test=10;%input('请输入需要预测个数：');
i=1:t_test+n;
X0s(i)=c*exp(-a*i);
x=1:n;
xs=2:n+t_test;
X0ss=X0s(2:n+t_test);
plot(x,X0,'^r',xs,X0ss,'*-b');
det=0;
for i=2:n
det=det+abs(X0s(i)-X0(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
disp(['预测值为： ',num2str(X0s(n+1:n+t_test))]);
改进程序2
对
中的μ进行调节，不在是取均值（μ=2）。

直到误差达到最小
clear
jm=xlsread('测试温度');
X0=jm;%input('请输入数据');%输入数据请用如例所示形式：[48.7 57.17 68.76 92.15] for p=0:0.01:1;
n=length(X0);
X1=ones(n,1);
X1(1)=X0(1);
for i=2:n
X1(i)=X1(i-1)+X0(i);
end
B=ones(n-1,2);
for i=1:(n-1)
B(i,1)=-(p*X1(i)+(1-p)*X1(i+1));
B(i,2)=1;
end
BT=B';
for j=1:n-1
Yn(j)=X0(j+1);
end
Yn=Yn';
A=(inv(BT*B)*BT)*Yn;
a=A(1);
u=A(2);
t=u/a;
for i=1:n-1
h1=0;
h1=h1+X0(i+1)*exp(-a*i);
h2=0;
h2=h2+exp(-2*a*i);
end
h3=(X0(1)-t)*(1-exp(a)^(-1));
h4=(1-exp(a))^(-2);
%c=((X0(1)-t)*(1-exp(a))+h1*((1-exp(a))^2))/(1+h2*((1-exp(a))^2)); c=(h1+h3)/(h2+h4);
t_test=10;%input('请输入需要预测个数：');
i=1:t_test+n;
X0s(i)=c*exp(-a*i);
x=1:n;
xs=2:n+t_test;
X0ss=X0s(2:n+t_test);
figure
plot(x,X0,'^r',xs,X0ss,'*-b');
det=0;
for i=2:n
det=det+abs(X0s(i)-X0(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
disp(['预测值为：',num2str(X0s(n+1:n+t_test))]);
end
多变量GM（1，m）
clear all
clc
%输入待预测时刻k 及原始序列X0
X0=[48.1 166.2; 509.3 175.9; 532.9 200.1;544.5 258.1;554.2 293.5;576.7 334.8;617.3 376.4;618.2 405.9;623.5 421.3;614.7 390.1;621.0 389.7];
%对原始序列X0 累加生成序列X1
[n,m]=size(X0);
t_test=input('请输入需要预测个数：');
k=n+t_test;
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);
X0s(1:m,1)=S';
for k=2:n+t_test
X0s(1:m,k)=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);
end
X0ss=X0s(1:m,n+1:n+t_test)。