文章目录 前言一、灰色预测模型灰色预测适用情况GM (1,1)模型 二、示例指数规律检验(原始数据级比检验)级比检验的定义GM(1,1) 模型的级比检验 模型求解求解微分方程 模型评价(检验模型对原始数据的拟合程度)残差检验级比偏差检验 三、代码实现----Matlab级比检验代码模型求解代码调用模型求解代码进行预测
前言
通过模型算法，熟练对Matlab的应用。 学习视频链接： https://www.bilibili.com/video/BV1EK41187QF?p=48&vd_source=67471d3a1b4f517b7a7964093e62f7e6
一、灰色预测模型
灰色系统理论在另一篇博客中已经阐述过了，此处不作赘述。
参考链接：【学习笔记】Matlab和python双语言的学习(灰色关联分析法）
灰色预测
所谓灰色预测，就是对既含有已知信息又含有不确定信息的系统进行预测，就是对在一定范围内变化的、与时间有关的灰色过程进行预测。
适用情况
数据是以年份度量的非负数据(如果是月份或者季度数据一般要用时间序列模型)，比如定时求量的题目，即已知一些年份数据，预测下一年的数据，常见有GDP、人口数量、耕地面积、粮食产量等；或者定量求时，已知一些年份数据和某灾变的阈值，预测下次灾变时间。数据能经过准指数规律的检验(除了前两期外，后面至少90%的期数的光滑比要低于0.5)。数据的期数较短且和其他数据之间的关联性不强(小于等于10，也不能太短了，比如只有3期数据)，要是数据期数较长，一般用传统的时间序列模型比较合适。
GM (1,1)模型
G M : G r e y m o d e l GM:Greymodel GM:Greymodel 灰色模型， ( 1 , 1 ) : (1,1){:} (1,1): 只含有一个变量的一阶微分方程模型如何用 G M ( 1 , 1 ) GM(1,1) GM(1,1) 进行灰色预测？ 根据原始的离散非负数据列，通过累加等方式削弱随机性、获得有规律的离散数据列建立相应的微分方程模型，得到离散点处的解再通过累减求得的原始数据的估计值，从而对原始数据预测
二、示例
长江在 1995-2004 年废水排放总量如下，如果不采取保护措施，请对今后水质污染发展趋势做出预测。
指数规律检验(原始数据级比检验)
级比检验的定义
为了确定原始数据是否可使用灰色预测模型，需要对原始数据进行级比检验要使用灰色预测数据首先应具有准指数规律： 累加 r r r 次的序列为： x ( r ) = ( x ( r ) ( 1 ) , x ( r ) ( 2 ) , . . . , x ( r ) ( n ) ) x^{(r)}=\left(x^{(r)}(1),x^{(r)}(2),...,x^{(r)}(n)\right) x(r)=(x(r)(1),x(r)(2),...,x(r)(n))，定义级比 σ ( k ) = x ( r ) ( k ) x ( r ) ( k − 1 ) , k = 2 , 3 , . . . , n \sigma(k)=\frac{x^{(r)}(k)}{x^{(r)}(k-1)},k=2,3,...,n σ(k)=x(r)(k−1)x(r)(k),k=2,3,...,n 对于 ∀ k , σ ( k ) ∈ [ a , b ] \forall k,\sigma(k)\in[a,b] ∀k,σ(k)∈[a,b]，且区间长度 δ = b − a < 0.5 \delta=b-a<0.5 δ=b−a<0.5，则称累加 r r r 次后的序列具有准指数规律
GM(1,1) 模型的级比检验
对于 G M ( 1 , 1 ) GM(1,1) GM(1,1) 模型中，我们只需要判断累加一次后的序列 x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , … , x ( 1 ) ( n ) ) x^{(1)}=\left(x^{(1)}(1),x^{(1)}(2),\ldots,x^{(1)}(n)\right) x(1)=(x(1)(1),x(1)(2),…,x(1)(n)) 是否具有准指数规律根据上述公式：序 列 x ( 1 ) x^{( 1) } x(1) 的级比 σ ( k ) = x ( 1 ) ( k ) x ( 1 ) ( k − 1 ) = x ( 0 ) ( k ) + x ( 1 ) ( k − 1 ) x ( 1 ) ( k − 1 ) = x ( 0 ) ( k ) x ( 1 ) ( k − 1 ) + 1 \sigma ( k) = \frac {x^{( 1) }( k) }{x^{( 1) }( k- 1) }= \frac {x^{( 0) }( k) + x^{( 1) }( k- 1) }{x^{( 1) }( k- 1) }= \frac {x^{( 0) }( k) }{x^{( 1) }( k- 1) }+ 1 σ(k)=x(1)(k−1)x(1)(k)=x(1)(k−1)x(0)(k)+x(1)(k−1)=x(1)(k−1)x(0)(k)+1定义 ρ ( k ) = x ( 0 ) ( k ) x ( 1 ) ( k − 1 ) \rho(k)=\frac{x^{(0)}(k)}{x^{(1)}(k-1)} ρ(k)=x(1)(k−1)x(0)(k) 为原始序列 x ( 0 ) x^{(0)} x(0) 的光滑比，注意到 ρ ( k ) = x ( 0 ) ( k ) x ( 0 ) ( 1 ) + x ( 0 ) ( 2 ) + ⋯ + x ( 0 ) ( k − 1 ) \rho(k)=\frac{x^{(0)}(k)}{x^{(0)}(1)+x^{(0)}(2)+\cdots+x^{(0)}(k-1)} ρ(k)=x(0)(1)+x(0)(2)+⋯+x(0)(k−1)x(0)(k)，假设 x ( 0 ) x^{(0)} x(0) 为非负序列(生活中常见的时间序列几乎都满足非负性)，那么随着 k k k 增加，最终 ρ ( k ) \rho(k) ρ(k) 会逐渐接近0，因此要使得具有 x ( 1 ) x^{(1)} x(1) 具有准指数规律，即 ∀ k \forall k ∀k，区间长度 δ < 0.5 \delta<0.5 δ<0.5，只需要保证 ρ ( k ) ∈ ( 0 , 0.5 ) \rho(k)\in(0,0.5) ρ(k)∈(0,0.5) 即可，此时序列 x ( 1 ) x^{(1)} x(1) 的级比 σ ( k ) ∈ ( 1 , 1.5 ) \sigma(k)\in(1,1.5) σ(k)∈(1,1.5)。实际建模中，我们要计算出 ρ ( k ) ∈ ( 0 , 0.5 ) \rho(k)\in(0,0.5) ρ(k)∈(0,0.5) 的占比，占比越高越好(一般 ρ ( 2 ) \rho(2) ρ(2)和 ρ ( 3 ) \rho(3) ρ(3)可能不符合，重点关注后面的数据)
模型求解
观察发现，没有明显规律，数据也比较少， 而且是以年份度量的，可以考虑用灰色预测那看不出规律怎么办，可以制造规律：
设 x ( 0 ) = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , . . . , x ( 0 ) ( n ) ) x^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)) x(0)=(x(0)(1),x(0)(2),...,x(0)(n))是最初的非负数据列，我们可以对其累加，得到新的数据列 x ( 1 ) x^{(1)} x(1) x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 1 ) ( 2 ) , … , x ( 1 ) ( n ) ) x^{(1)}=\left(x^{(1)}(1),x^{(1)}(2),\ldots,x^{(1)}(n)\right) x(1)=(x(1)(1),x(1)(2),…,x(1)(n)) 其中： x ( 1 ) ( m ) = ∑ i = 1 m x ( 0 ) ( i ) , m = 1 , 2 , . . . , n x^{(1)}(m)=\sum_{\mathrm{i}=1}^{\mathrm{m}}x^{(0)}(i),m=1,2,...,n x(1)(m)=∑i=1mx(0)(i),m=1,2,...,n
观察可知，新序列 x ( 1 ) x^{(1)} x(1) 曲线像一个指数曲线(直线)我们可以用指数曲线的表达式来逼近序列 x ( 1 ) x^{(1)} x(1)，相应可以构建一阶常微分方程来求解拟合指数曲线的函数表达式一阶常微分方程 d x ( 1 ) d t + a x ( 1 ) = u \frac{dx^{(1)}}{dt}+ax^{(1)}=u dtdx(1)+ax(1)=u 要求出 x ( 1 ) x^{(1)} x(1) 的表达式，就需要解出常微分方程，所以要先知道参数 a a a 和 u u u。
求解微分方程
一阶常微分方程 d x ( 1 ) d t + a x ( 1 ) = u \frac{dx^{(1)}}{dt}+ax^{(1)}=u dtdx(1)+ax(1)=u已知，我们的数据是离散的，所以 d x ( 1 ) d t \frac {dx^{(1)}}{dt} dtdx(1) 等同于 Δ x ( 1 ) Δ t = Δ x ( 1 ) t − ( t − 1 ) = Δ x ( 1 ) = x ( 1 ) ( t ) − x ( 1 ) ( t − 1 ) = x ( 0 ) ( t ) \frac {\Delta x^{( 1) }}{\Delta t}= \frac {\Delta x^{(1)}}{t-(t-1)}= \Delta x^{(1)}= x^{(1)}(t) - x^{(1)}(t-1) = x^{(0)}(t) ΔtΔx(1)=t−(t−1)Δx(1)=Δx(1)=x(1)(t)−x(1)(t−1)=x(0)(t) 则微分方程变为 x ( 0 ) ( t ) + a x ( 1 ) ( t ) = u x^{(0)}(t)+ax^{(1)}(t)=u x(0)(t)+ax(1)(t)=u上式为常见的一元线性方程，为了消除数据随机性，定义 z ( 1 ) = ( z ( 1 ) ( 1 ) , z ( 1 ) ( 2 ) , . . . , z ( 1 ) ( n ) ) z^{(1)}=(z^{(1)}(1),z^{(1)}(2),...,z^{(1)}(n)) z(1)=(z(1)(1),z(1)(2),...,z(1)(n)) 其中： z ( 1 ) ( m ) = δ x ( 1 ) ( m ) + ( 1 − δ ) x ( 1 ) ( m − 1 ) , m = 2 , 3 , . . . , n z^{(1)}(m)=\delta x^{(1)}(m)+(1-\delta)x^{(1)}(m-1),m=2,3,...,n z(1)(m)=δx(1)(m)+(1−δ)x(1)(m−1),m=2,3,...,n 且 δ = 0.5 \delta=0.5 δ=0.5，即为前后时刻的均值则微分方程改为 x ( 0 ) ( t ) = − a z ( 1 ) ( t ) + u x^{(0)}(t)=-az^{(1)}(t)+u x(0)(t)=−az(1)(t)+u我们已知 x ( 0 ) ( t ) , z ( 1 ) ( t ) x^{(0)}(t),z^{(1)}(t) x(0)(t),z(1)(t)的数据，结合线性回归的知识，可利用线性回归或者用最小二乘法求解参数把 a ^ \hat{a} a^ 和 u ^ \hat{u} u^ 代入微分方程 d x ( 1 ) d t + a x ( 1 ) = u \frac{dx^{(1)}}{dt}+ax^{(1)}=u dtdx(1)+ax(1)=u，并求解可得 x ^ ( 1 ) ( m + 1 ) = [ x ( 0 ) ( 1 ) − u ^ a ^ ] e − a ^ m + u ^ a ^ , m = 1 , 2 , . . . , n − 1 \hat{x}^{(1)}(m+1)=\left[x^{(0)}(1)-\frac{\hat{u}}{\hat{a}}\right] e^{-\hat{a}m}+\frac{\hat{u}}{\hat{a}},m=1,2,...,n-1 x^(1)(m+1)=[x(0)(1)−a^u^]e−a^m+a^u^,m=1,2,...,n−1由于 x ( 1 ) ( m ) = ∑ i = 1 m x ( 0 ) ( i ) , m = 1 , 2 , . . . , n x^{(1)}(m)=\sum_{\mathrm{i}=1}^{\mathrm{m}}x^{(0)}(i),m=1,2,...,n x(1)(m)=∑i=1mx(0)(i),m=1,2,...,n，所以我们可以得到： x ^ ( 0 ) ( m + 1 ) = x ^ ( 1 ) ( m + 1 ) − x ^ ( 1 ) ( m ) = ( 1 − e a ^ ) [ x ( 0 ) ( 1 ) − u ^ a ^ ] e − a ^ m , m = 1 , 2 , . . . , n − 1 \hat{x}^{(0)}(m+1)=\hat{x}^{(1)}(m+1)-\hat{x}^{(1)}(m)=\left(1-e^{\hat{a}}\right)\left[x^{(0)}(1)-\frac{\hat{u}}{\hat{a}}\right]e^{-\hat{a}m},m=1,2,...,n-1 x^(0)(m+1)=x^(1)(m+1)−x^(1)(m)=(1−ea^)[x(0)(1)−a^u^]e−a^m,m=1,2,...,n−1当 m m m 取 0,1,…,9 时，得到的 x ^ ( 0 ) \hat{x}^{(0)} x^(0) 为拟合值，大于 9 时，得到的为预测值
模型评价(检验模型对原始数据的拟合程度)
残差检验
绝对残差： ϵ ( k ) = x ( 0 ) ( k ) − x ^ ( 0 ) ( k ) , k = 2 , 3 , . . . , n \epsilon(k)=x^{(0)}(k)-\hat{x}^{(0)}(k),k=2,3,...,n ϵ(k)=x(0)(k)−x^(0)(k),k=2,3,...,n相 对 残 差： ϵ r ( k ) = ∣ x ( 0 ) ( k ) − x ^ ( 0 ) ( k ) ∣ x ( 0 ) ( k ) × 100 % , k = 2 , 3 , . . . , n \epsilon _r( k) = \frac {\left | x^{( 0) }( k) - \hat{x} ^{( 0) }( k) \right | }{x^{( 0) }( k) }\times 100\% , k= 2, 3, . . . , n ϵr(k)=x(0)(k)∣x(0)(k)−x^(0)(k)∣×100%,k=2,3,...,n平均相对残差： ϵ ˉ r = 1 n − 1 ∑ k = 2 n ∣ ϵ r ( k ) ∣ \bar{\epsilon}_r=\frac1{n-1}\sum_{k=2}^n|\epsilon_r(k)| ϵˉr=n−11∑k=2n∣ϵr(k)∣ 如果 ϵ ˉ r < 20 % \bar{\epsilon}_r<20\% ϵˉr<20%，则认为 G M ( 1 , 1 ) GM(1,1) GM(1,1) 对原数据的拟合达到一般要求如果 ϵ ˉ r < 10 % \bar{\epsilon}_r<10\% ϵˉr<10%，则认为 G M ( 1 , 1 ) GM(1,1) GM(1,1) 对原数据的拟合效果非常不错
级比偏差检验
首先计算原始数据级比 σ ( k ) = x ( 0 ) ( k ) x ( 0 ) ( k − 1 ) ( k = 2 , 3 , … , n ) \sigma(k)=\frac{x^{(0)}(k)}{x^{(0)}(k-1)}\left(k=2,3,\ldots,n\right) σ(k)=x(0)(k−1)x(0)(k)(k=2,3,…,n)
再根据预测发展系数 ( − a ^ ) (-\hat{a}) (−a^) 计算级比偏差和平均级比偏差 η ( k ) = ∣ 1 − 1 − 0.5 a ^ 1 + 0.5 a ^ 1 σ ( k ) ∣ , η ‾ = ∑ k = 2 n η ( k ) / ( n − 1 ) \eta(k)=\left|1-\frac{1-0.5\hat{a}}{1+0.5\hat{a}}\frac{1}{\sigma(k)}\right|,\:\overline{\eta}=\sum_{k=2}^{n}\eta(k)/(n-1) η(k)= 1−1+0.5a^1−0.5a^σ(k)1 ,η=k=2∑nη(k)/(n−1) 如果 η ˉ < 0.2 \bar{\eta}<0.2 ηˉ<0.2，则认为模型对原数据拟合达到一般要求； η ˉ < 0.1 \bar{\eta}<0.1 ηˉ<0.1，则认为模型对原数据拟合效果非常不错
三、代码实现----Matlab
级比检验代码
clear;clc
year =[1995:1:2004]'; % 横坐标表示年份，写成列向量的形式
x0 = [174,179,183,189,207,234,220.5,256,270,285]'; %原始数据序列，写成列向量的形式n = size(x0,1);
x1 = zeros(n,1);
for i = 1:nx1(i) = sum(x0(1:i,1));
end% 级比检验
rho = zeros(n,1);
for i = 2:nrho(i) = x0(i) / x1(i-1);
end
figure(2)
plot(year(2:n,1),rho(2:n,1),'o-',[year(2),year(n)],[0.5,0.5],'-'); grid on;
text(year(end-1)+0.2,0.55,'临界线') % 在坐标(year(end-1)+0.2,0.55)上添加文本
set(gca,'xtick',year(2:1:end)) % 设置x轴横坐标的间隔为1
xlabel('年份'); ylabel('原始数据的光滑度'); % 给坐标轴加上标签
% 指标1：光滑比小于0.5的数据占比
num1 = sum(rho<0.5)/(n-1)
% 指标2：除去前两个时期外，光滑比小于0.5的数据占比
num2 = sum(rho(3:end)<0.5)/(n-3)
运行结果：
参考标准：指标1一般要大于60%, 指标2要大于90%，所以通过了级比检验。
模型求解代码
function [result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num)n = size(x0,1);x1 = cumsum(x0);z1 = 0.5 * x1(2:end) + 0.5 * x1(1:n-1,1);y = x0(2:end);x = z1;% 最小二乘法求解k = ((n-1)*sum(x.*y)-sum(x)*sum(y))/((n-1)*sum(x.*x)-sum(x)*sum(x));b = (sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/((n-1)*sum(x.*x)-sum(x)*sum(x));a = -k;u = b;x0_hat = zeros(n,1);x0_hat(1) = x0(1);for m = 1: n-1x0_hat(m+1) = (1-exp(a))*(x0(1)-u/a)*exp(-a*m);endresult = zeros(predict_num,1); % 初始化用来保存预测值的向量for i = 1: predict_numresult(i) = (1-exp(a))*(x0(1)-u/a)*exp(-a*(n+i-1)); % 带入公式直接计算end% 计算绝对残差和相对残差absolute_residuals = x0(2:end) - x0_hat(2:end); % 从第二项开始计算绝对残差，因为第一项是相同的relative_residuals = abs(absolute_residuals) ./ x0(2:end); % 计算相对残差% 计算级比和级比偏差class_ratio = x0(2:end) ./ x0(1:end-1) ; % 计算级比eta = abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio)); % 计算级比偏差
end
% 函数作用：使用传统的GM(1,1)模型对数据进行预测
% x0：要预测的原始数据
% predict_num： 向后预测的期数% 输出变量
% result：预测值
% x0_hat：对原始数据的拟合值
% relative_residuals： 对模型进行评价时计算得到的相对残差
% eta： 对模型进行评价时计算得到的级比偏差
调用模型求解代码进行预测
if num1 > 0.6 && num2 > 0.9if n > 7 % 将数据分为训练组和试验组(根据原数据量大小n来取，n小于7则取最后两年为试验组，n大于7则取最后三年为试验组)test_num = 3;elsetest_num = 2;endtrain_x0 = x0(1:end-test_num); % 训练数据disp('训练数据是: ')disp(mat2str(train_x0')) % mat2str可以将矩阵或者向量转换为字符串显示test_x0 = x0(end-test_num+1:end); % 试验数据disp('试验数据是: ')disp(mat2str(test_x0'))% 使用GM（1，1）模型对训练数据进行训练disp('GM(1,1)模型预测')result1 = gm11(train_x0, test_num); %预测后test_num期的结果% 绘制对试验数据进行预测的图形test_year = year(end-test_num+1:end); % 试验组对应的年份figure(3)plot(test_year,test_x0,'o-',test_year,result1,'*-'); grid on;set(gca,'xtick',year(end-test_num+1): 1 :year(end))legend('试验组的真实数据','GM(1,1)预测结果') xlabel('年份'); ylabel('排污总量'); % 给坐标轴加上标签predict_num = input('请输入你要往后面预测的期数： ');% 计算使用传统GM模型的结果[result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num);%% 绘制相对残差和级比偏差的图形figure(4)subplot(2,1,1) % 绘制子图（将图分块）plot(year(2:end), relative_residuals,'*-'); grid on; % 原数据中的各时期和相对残差legend('相对残差'); xlabel('年份');set(gca,'xtick',year(2:1:end)) % 设置x轴横坐标的间隔为1subplot(2,1,2)plot(year(2:end), eta,'o-'); grid on; % 原数据中的各时期和级比偏差legend('级比偏差'); xlabel('年份');set(gca,'xtick',year(2:1:end))%% 残差检验average_relative_residuals = mean(relative_residuals); % 计算平均相对残差 mean函数用来均值disp(strcat('平均相对残差为',num2str(average_relative_residuals)))%% 级比偏差检验average_eta = mean(eta); % 计算平均级比偏差disp(strcat('平均级比偏差为',num2str(average_eta)))%% 绘制最终的预测效果图figure(5) plot(year,x0,'-o', year,x0_hat,'-*m', year(end)+1:year(end)+predict_num,result,'-*b' ); grid on;hold on;plot([year(end),year(end)+1],[x0(end),result(1)],'-*b')legend('原始数据','拟合数据','预测数据') set(gca,'xtick',[year(1):1:year(end)+predict_num]) xlabel('年份'); ylabel('排污总量');
end
运行结果：
可以看出：平均相对残差为：0.025999 ＜ 10 % 级比偏差为 0.047041 ＜ 0.1 拟合效果非常不错
2024/8/25 更新
在求解微分方程时，我们把微分方程转化成了一元线性回归，求解出回归系数后再求解微分方程。
求解回归系数的过程:
我们可以直接调用 Matlab 的 regress 函数，调用时需注意， x x x 需加上常数项。还可以使用最小二乘法求解。最小二乘法的通解形式为： β ^ = ( X T X ) − 1 ( X T Y ) \hat{\beta}=(X^{T}X)^{-1}(X^{T}Y) β^=(XTX)−1(XTY)或者直接套用公式： k = ( n − 1 ) ∑ ( x i y i ) − ∑ x i ∑ y i ( n − 1 ) ∑ ( x i 2 ) − ( ∑ x i ) 2 k=\frac{(n-1)\sum(x_iy_i)-\sum x_i\sum y_i}{(n-1)\sum(x_i^2)-(\sum x_i)^2} k=(n−1)∑(xi2)−(∑xi)2(n−1)∑(xiyi)−∑xi∑yi b = ∑ ( x i 2 ) ∑ y i − ∑ x i ∑ ( x i y i ) ( n − 1 ) ∑ ( x i 2 ) − ( ∑ x i ) 2 b=\frac{\sum(x_i^2)\sum y_i-\sum x_i\sum(x_iy_i)}{(n-1)\sum(x_i^2)-(\sum x_i)^2} b=(n−1)∑(xi2)−(∑xi)2∑(xi2)∑yi−∑xi∑(xiyi)该公式也是通过最小二乘法推导得到的。