Matlab灰色预测模型GM(1,1)代码

function c7fun73
X0=[2.874 3.278 3.307 3.39 3.679];
AU=c7fun73(X0);a=AU(1);u=AU(2);m2=length(X0);
for k=1:1:m2-1
xx1(k+1)=(X0(1)-u/a)*exp(-a*k)+u/a;
end
s=0;xx0(1)=X0(1);
for jj=2:1:m2;
xx0(jj)=xx1(jj)-xx1(jj-1);
end
disp('GM(1,1)对数列进行预测结果');xx0
disp('数列1原始观测数据');X0
disp('a');AU(1)
disp('u');AU(2)
function au=c7fun73(X0)
m=length(X0);
s1=0;
for jj=1:1:m;X1(jj)=s1+X0(jj); s1=X1(jj); end
for ii=1:1:m-1; B(ii)=-(X1(ii)+X1(ii+1))/2; end
B=[B(:),ones(m-1,1)];
y=X0([2:m])';
au=inv((B'*B))*B'*y;



% 用MATLAB的灰色预测GM（1，1）模型
%
%程序中的变量定义；alpha是包含 值的矩阵；ago是预测后累加值矩阵；var是预测值矩阵；error是残差矩阵；c是后验差比值
function my_gm_test()
clc
data=exprnd(5,1,10)%原始10个数据
[ago alpha var error c]=gm1(data)
plot(data)
hold on
plot(var','-r*')
function [ago alpha var error c]=gm1(x); %定义函数gm1(x)
%清屏，以使结果独立显示
format long; %设置计算精度
if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换
x=x';
end
n=length(x); %取输入数据的样本量
z=0;
for i=1:n %计算累加值，并将值赋予矩阵be
z=z+x(i,:);
be(i,:)=z;
end
for i=2:n %对原始数列平行移位
y(i-1,:)=x(i,:);
end
for i=1:n-1 %计算数据矩阵B的第一列数据
c(i,:)=-0.5*(be(i,:)+be(i+1,:));
end
for j=1:n-1 %计算数据矩阵B的第二列数据
e(j,:)=1;
end
for i=1:n-1 %构造数据矩阵B
B(i,1)=c(i,:);
B(i,2)=e(i,:);
end
alpha=inv(B'*B)*B'*y; %计算参数 矩阵
for i=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值
ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);
end
var(1,:)=ago(1,:)
for i=1:n %如改n为n+m-1，可预测后m-1个值
var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值
end
for i=1:n
error(i,:)=var(i,:)-x(i,:); %计算残差
end
c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c
ago %显示输出预测值的累加数列
alpha %显示输出参数 数列
var %显示输出预测值
error %显示输出误差
c %显示后验差的比值c



function []=greymodel(y)
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
yy(i)=yy(i-1)+y(i);
end
B=ones(n-1,2);
for i=1:(n-1)
B(i,1)=-(yy(i)+yy(i+1))/2;
B(i,2)=1;
end
BT=B';
for j=1:n-

1
YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
u=A(2);
t=u/a;
i=1:n+2;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+2:-1:2
ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+2;
yn=ys(2:n+2);
plot(x,y,'^r',xs,yn,'*-b');
det=0;
for i=2:n
det=det+abs(yn(i)-y(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
disp(['下个拟合值为 ',num2str(ys(n+1))]);
disp(['再下个拟合值为',num2str(ys(n+2))]);



x(0)={x(0)(1),x(0)(2),……,x(0)(N)}=
{ 724.57, 746.62, 778.27, 800.8, 827.75,871.1, 912.37, 954.28, 995.01, 1037.2}