🔆 文章首发于我的个人博客：欢迎大佬们来逛逛

数学建模：灰色预测模型

灰色预测

三个基本方法：

累加数列：计算一阶累加生成数列

$$x^{(1)}(k)=\sum _{i=1}^k{x}^{(0)}(i),k=1,2,\cdots ,n,$$

---

累减数列：计算一阶累减生成数列

$$x^{(0)}(k)=x^{(1)}(k)-x^{(1)}(k-1),k=2,3,\cdots ,n,$$

---

加权累加：计算一阶等权邻接生成数，图片描述有误，此处计算的是一次累加的加权邻值生成

$$z^{(0)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1),k=2,3,\cdots ,n,$$

---

算法步骤

1. 进行级比检验，检查是否满足建立微分方程的前提条件。

$$\lambda (k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)}$$

1. 对原数据做一次累加，计算加权邻值生成数
2. 构造数据矩阵 $B$ ，与数据向量 $Y$

$$B=\left[\begin{array}{cccc}-\frac{1}{2}\left(x^{(1)}\right(1)+x^{(1)}(& 2)& & \\ -\frac{1}{2}\left(x^{(1)}\right(2)+x^{(1)}(& 3)& )& 1\\ & ⋮& & ⋮\\ -\frac{1}{2}\left(x^{(1)}\right(n-1)+x^{(1)}(& n)& )& \end{array}\right],Y=\left[\begin{array}{cc}{x}^{(0)}(& 2)\\ x^{(0)}(& 3)\\ ⋮& \\ x^{(0)}(& n)\end{array}\right]$$

1. 计算 $a$ 与 $b$ 的值

$$\hat{u}=(\hat{a},\hat{b}{)}^T=(B^T\cdot B{)}^{-1}{B}^{T}Y$$

1. 构建模型

$$x^{(1)}(t)=(x^{(0)}(1)-\frac{b}{a})e^{-a(t-1)}+\frac{b}{a}.$$

1. 计算生成模型值 ${\hat{x}}^{(1)}(k)$ 和模型还原值 ${\hat{x}}^{(0)}(k)$ 并且带入预测

$${\hat{x}}^{(0)}(k)={\hat{x}}^{(1)}(k)-{\hat{x}}^{(1)}(k-1)$$

1. 检验预测值

---

代码实现

%95至04年数据
clc;clear;
data = [174 179 183 189 207 234 220.5 256 270 285];
n = length(data);
 
%% 级比检验通过
check = [];
for k = 2:n
    lambda(k) = data(k-1)/data(k);
    if (exp(-2/(n+1))<lambda(k))&&(lambda(k)<exp(2/(n+1)))
        check(end+1) = 1;
    else check(end+1) = 0;
    end
end 
 
%% 计算累加数列
X1 = cumsum(data);

%% 计算加权
for i=2:n
    z(i) = 0.5*(X1(i-1)+X1(i));
end
 
%% 数据矩阵B及数据向量Y
Y = data(2:n)';
B = [-z(2:n)',ones(n-1,1)];
u = (B'*B)\B'*Y;
% u = B\Y; 表示B的逆 乘以 Y
a = u(1,1);
b = u(2,1);
 
%% 构造模型并且带入预测值
% 生成预测一次累加数列
f_X1 = [];
f_X0 = [];
for k=1:n-1
    f_X1(1)=data(1);
    f_X1(k+1) = (data(1)-b/a)*exp(-a*k) + b/a;
end
% 前缀和反推原始数据
for k=2:n
    f_X0(1)=data(1);
    f_X0(k)=f_X1(k)-f_X1(k-1);
end
 
%% 残差检验 与 级比偏差值检验
for k=1:n-1
    sigma(k)=abs((data(k)-f_X0(k))/data(k));
    rho(k+1)=abs(1-((1-0.5*a)*lambda(k+1))/(1+0.5*a));
end

%% 预测下n个值
test = input('nums:');
nums = 5;
n=n+test;
f_f_X1 = [];
f_f_X0 = [];
for k=1:n-1
    f_f_X1(1)=data(1);
    f_f_X1(k+1) = (data(1)-b/a)*exp(-a*k) + b/a;
end
for k=2:n
    f_f_X0(1)=data(1);
    f_f_X0(k)=f_f_X1(k)-f_f_X1(k-1);
end

%% 绘图
xAxis = 1995:2004;
xAxisPredict = 1995:1995+n-1; 
h = plot(xAxis,data,'o',xAxisPredict,f_f_X0,'-');
set(gca, 'XScale', 'log', 'YScale', 'log');
set(h,'LineWidth',1.5);