%%%%%%%%%%%%%%%%把下面函数保存为gmcal.m文件%%%%%%%%%%%
function gmcal=gm1(x)
sizexd2 = size(x,2);
%求数组长度
k=0;
for y1=x
k=k+1;
if k>1
x1(k)=x1(k-1)+x(k);
%累加生成
z1(k-1)=-0.5*(x1(k)+x1(k-1));
%z1维数减1，用于计算B
yn1(k-1)=x(k);
else
x1(k)=x(k);
end
end
%x1,z1,k,yn1
sizez1=size(z1,2);
%size(yn1);
z2 = z1';
z3 = ones(1,sizez1)';
YN = yn1'; %转置
%YN
B=[z2 z3];
au0=inv(B'*B)*B'*YN;
au = au0';
%B,au0,au
afor = au(1);
ufor = au(2);
ua = au(2)./au(1);
%afor,ufor,ua
%输出预测的 a u 和 u/a的值
constant1 = x(1)-ua;
afor1 = -afor;
x1t1 = 'x1(t+1)';
estr = 'exp';
tstr = 't';
leftbra = '(';
rightbra = ')';
%constant1,afor1,x1t1,estr,tstr,leftbra,rightbra
strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)
%输出时间响含信皮掘应方程谈握轮,也就是最终要求的灰色模型
%%%%%%%%%%%%%%%%%%%%%在workspace里输入%%%%%%%%%%%%
x =[5999,5903,5848,5700,7884];gm1(x)
%其中5999,5903,5848,5700,7884可以换成已知的
历史数据，无论几个都可以。
%灰色预测模型
function y=gm11(x,n)
%x为行向量数据
%做一次累加
x1=zeros(size(x));
for i=1:size(x1,2)
x1(i)=sum(x(1:i));
end
%x1的均值数列
z1=zeros(size(x));
for i=1:size(x1,2)-1
z1(i+1)=0.5*x1(i+1)+0.5*x1(i);
end
Yn=x(2:size(x,2))';
B=-z1(2:size(z1,2))';
B(:,2)=1;
u=inv((B'*B))*B'*Yn;
a=u(1);
b=u(2);
%预测
x2=zeros(1,n);
x2(1)=x(1);
for i=1:n-1
x2(1+i)=(x(1)-b/a)*exp(-a*i)+b/a;
end
x2=[0 x2];
y=diff(x2);
用法：
假设数列1 2 3 4 5.5 6 7.5 为已基没知数据，搏磨纳你要预测游手后面3项，那么保存上面的代码后输入命令：
gm11([1 2 3 4 5.5 6 7.5],10)
clear;
clc;
x0=[127627 128453 129227 129988 130756];
n=length(x0);
lamda=x0(1:n-1)./x0(2:n)
range=minmax(lamda)
x1=cumsum(x0)
for i=2:n
z(i)=0.5*(x1(i)+x1(i-1));
end
B=[-z(2:n)',ones(n-1,1)];
Y=x0(2:n)'塌携誉;
u=B\隐裂Y
x=dsolve('Dx+a*x=b','x(0)=x0');
x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});
yuce1=subs(x,'t',[0:n-1]);
digits(6),y=vpa(x) %为提高预测精度，先计算预测值，再显示微分方程的解
yuce=[x0(1),diff(yuce1)]
epsilon=(x0-yuce)./x0 %计算残差
delta=abs(epsilon./x0) %计算相对团段误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda %计算级比偏差值