GM(1,1)灰色预测详解+MATLAB预测下一年噪声值

本文为北海的数模课程学习笔记，课程出自微信公众号：数学建模BOOM。

求赞！求收藏！求关注！

目录

 GM(1,1)灰色预测模型

名词解释

灰色系统例子——GDP

如何进行灰色预测

适用赛题

典型例题与原理讲解

原始数据的级比检验

例题与累加“制造规律”

 问题转化为求微分方程中的参数​编辑

表达式处理

 最小二乘法求解

 模型求解

模型检验

思路总结

代码求解

 GM(1,1)灰色预测模型

GM(1,1)灰色预测是灰色系统理论中的一种线性一阶微分方程模型，用于对具有不完全信息的数据序列进行趋势预测，特别适用于数据量小且变化规律不明显的情况。

名词解释

• 黑色系统就是黑盒。
• 白色系统就是白盒。
• 灰色系统就是介于黑与白直接。
• (1,1)是只含有一个变量的一阶微分方程模型

灰色系统例子——GDP

部分已知：往年数据已知。

部分未知：下一年数据未知。但是可以进行预测。

样本数据少：不确定性高。

如何进行灰色预测

• 根据原始数据，通过 累加等 方式削弱随机性、 获得有规律 的新序列

• 建立相应的 微分方程 模型，得到离散点处的解

• 再通过累减求得的原始数据的估计值，从而对原始数据 预测

适用赛题

❑ 数列预测

• 特点： 定时求量 ，已知 xx 年到 xx 年的数据，请预测下一年的数值

• 常见 GDP 、人口数量、耕地面积、粮食产量等问题

• 针对的问题往往短期波动小、可预测，但长期可能变化大、难以准确预测

❑ 灾变预测——可转换为数列预测，都是要解微分方程

• 特点： 定量求时 ，已知 xx 年到 xx 年的数据和某灾变的阈值，预测下一次灾变发生的时间

• 常见洪涝灾害、虫灾等问题

• 模型中需要把超出阈值的数据（异常数据）对应的时间组成新序列

❑ 拓扑预测——可转换为灾变预测

• 特点： 对数据波形进行预测 ，求的是多个模型构成的模型群， 等于求解多个灾变预测

• 与灾变预测类似，不过有较详细的分级，例如虫灾“轻微”“中度”“重度”——相当于划了不同的阈值线，是

❑ 注意事项

• 需要的数据量少 ，而且数据量太多了没意义，例如用近 100 年的 GDP 去预测下一年毫无意义

• 只能短期预测 ，究竟多短没有严格限制

典型例题与原理讲解

原始数据的级比检验

在建立模型的最开始做的。保证这个问题适合用灰色预测，避免白忙活。

例题与累加“制造规律”

 问题转化为求微分方程中的参数

表达式处理

 最小二乘法求解

 记住公式。

 模型求解

模型检验

 

• 如果指的是k取遍所有都小于0.2。
•  中的a是前面通过最小二乘法求得的a。
• 相对误差检验是检验预测效果的。
• 级比偏差检验是对a进行检验。

思路总结

• 核心就是最小二乘法求参数，再求微分方程

• 在最开始要对原始数据进行检验，在最后要对结果进行检验

代码求解

clc,clear
x0=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]'; % 注意这里为列向量
n=length(x0);
lamda=x0(1:n-1)./x0(2:n);  % 计算级比
range=minmax(lamda');      % 计算级比的范围
lamrange_min=exp(-2/(n+1));     % 允许的范围下界
lamrange_max=exp(2/(n+1));      % 允许的范围上界
% 原始数据的级比检验
if (range(1)>lamrange_min && range(2)<lamrange_max)
    fprintf("原始数据通过级比检验")
else
    fprintf("原始数据未能通过级比检验！！！")
    return
end

x1=cumsum(x0)  %累加序列
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
% 注意此处是反斜杠，即左除
u=B\Y       % 拟合参数u(1)=a,u(2)=b，左除，相当于B的逆矩阵乘Y

syms x(t)   %定义符号变量

%dsolve求微分方程的符号解，即累加序列的函数表达式，而不是数值
x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1));     % diff是一阶导或差分
xt=vpa(x,6) %以保留6位小数的格式显示微分方程的解,vpa设置精度

% 函数subs(x,t,[0:n-1])意味着把t=0到t=n-1依次代入x的表达式求出值
nihe1=subs(x,t,[0:n-1]); %求已知数据对应的拟合值
nihe1=double(nihe1);     %符号数转换成数值类型，否则无法作差分运算
nihe=[x0(1),diff(nihe1)] %差分运算，得到原始序列拟合值

yuce_result=double(subs(x,t,7))-double(subs(x,t,6))  % 原始序列的预测值，即第8年噪声

delta=abs((x0'-nihe )./x0')  %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'  %计算级比偏差值