概述算法：灰色预测模型用于对原始数据（≥4个）做中短期预测，其中，GM(1,1)模型适用于具有较强的指数规律的序列，只能描述单调的变化过程，而GM(2,1)模型适用于非单调的摆动发展序列或具有饱和的S形序列。
GM(1,1)编程步骤：
1.建立时间序列
2.检验数据是否符合要求
3.计算一次累加生成序列
4.计算邻均值等权数列
5.构造矩阵B、Y
6.计算发展系数a和灰作用量b
7.计算模型拟合值
8.模型精度评定（后验差检验）
①计算残差
②计算标准差
③计算后验差比值、小误差概率
④查表定级
GM(2,1)编程步骤与GM(1,1)类似。
下面就一起来看看如何将优雅的数学语言转换成matlab语言吧。
GM(1,1)源代码
代码语言：javascript复制clear;clc;
% 建立时间序列【输入】
x0 = [15.9 15.4 18.1 21.3 20.1 22.0 22.6 21.4]';
% 需要预测几期数据【输入】,预测数据见x0_hat变量
count = 2;
% 检验数据是否符合要求
n1 = length(x0);
lmd = x0(1:end-1)./x0(2:end);
flag = min(lmd)>exp(-2/(n1 1)) & max(lmd)<exp(2/(n1 1));
if ~flag
error('数据级比不符合要求');
end
% 计算一次累加生成序列
x1 = cumsum(x0);
% 计算邻均值等权数列
z1 = 0.5 * ( x1(1:end-1) x1(2:end) );
% 构造矩阵B、Y
B = [-z1,ones(n1-1,1)];
Y = x0(2:n1);
% 计算发展系数a和灰作用量b
u = (B'*B)(B'*Y);
a = u(1);
b = u(2);
% 计算模型拟合值
k = (1:n1-1 count)';
x0_hat = [x0(1);(1-exp(a))*(x0(1)-b/a)*exp(-a*k)];
disp('预测数据：')
x0_hat(n1 1:end)
% 模型精度评定
e = x0(1:n1)-x0_hat(1:n1);
s1 = std(x0);
s2 = std(e);
c = s2/s1;
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95 && c<=0.35
disp('精度等级:1级(好)');
elseif p>=0.80 && c<=0.5
disp('精度等级:2级(合格)');
elseif p>=0.7 && c<=0.65
disp('精度等级:3级(勉强)');
else
disp('精度等级:4级(不合格)');
end
% 绘图说明
plot(1:n1,x0,'ro','LineWidth',2.5);
hold on
plot(1:n1 count,x0_hat,'bo','LineWidth',2.5);
plot(1:n1 count,x0_hat,'b-','LineWidth',2.5);
hold off
legend('实测值','预测值','FontSize',18);
GM(2,1)代码
代码语言：javascript复制clear;clc;
% 建立时间序列【输入】
x0 = [5.6 4.2 3.3 2.5 3.1 4.4 5.8]';
n1 = length(x0);
% 需要预测几期数据【输入】,预测数据见x0_hat变量
count = 2;
% 计算一次累加生成序列
x1 = cumsum(x0);
% 计算一次累减生成序列
alpx0 = x0(2:end)-x0(1:end-1);
% 计算邻均值等权数列
z1 = 0.5 * ( x1(1:end-1) x1(2:end) );
% 构造矩阵B、Y
B = [ -x0(2:end),-z1,ones(n1-1,1)];
Y = alpx0;
% 计算发展系数a和灰作用量b
u = (B'*B)(B'*Y);
% 计算模型拟合值
syms x(t)
x=dsolve(diff(x,2) u(1)*diff(x) u(2)*x==u(3),x(0)==x1(1),x(n1-1)==x1(n1));
xt=vpa(x,2);
x1_hat=subs(x,t,0:n1-1 count);
x1_hat=(double(x1_hat))';
x0_hat = [x0(1);diff(x1_hat)];
disp('预测数据：')
x0_hat(n1 1:end)
% 模型精度评定
e = x0(1:n1)-x0_hat(1:n1);
s1 = std(x0);
s2 = std(e);
c = s2/s1;
p = length(find(e-mean(e)<0.6745*s1))/n1;
if p>=0.95 && c<=0.35
disp('精度等级:1级(好)');
elseif p>=0.80 && c<=0.5
disp('精度等级:2级(合格)');
elseif p>=0.7 && c<=0.65
disp('精度等级:3级(勉强)');
else
disp('精度等级:4级(不合格)');
end
% 绘图说明
plot(1:n1,x0,'ro','LineWidth',2.5);
hold on
plot(1:n1 count,x0_hat,'bo','LineWidth',2.5);
plot(1:n1 count,x0_hat,'b-','LineWidth',2.5);
hold off
legend('实测值','预测值','FontSize',18);
通过学习相关算法并将算法转变为实际的编程语言是练习编程的一种重要途径，这不仅可以提升理论认知，还能提高实践动手能力。鉴于此，matlab爱好者公众号计划推出【编程算法】系列，将逐一介绍各类算法在matlab中实现，与大家一起来在算法的海洋里畅游。
若您对算法感兴趣，并有一定的matlab编程基础，欢迎将所学算法整理成文推送给我们。
