【备战美赛】预测模型

最新推荐文章于 2024-05-30 16:17:47 发布

Fellyhosn

最新推荐文章于 2024-05-30 16:17:47 发布

阅读量1.1k

点赞数 1

分类专栏：建模文章标签：算法 matlab

本文链接：https://blog.csdn.net/Fellyhosn/article/details/122541172

版权

建模专栏收录该内容

6 篇文章 1 订阅

订阅专栏

预测模型指根据已有的数据对未来进行精准地预测，它包括灰度预测、插值与拟合、时间序列预测等方法。我们在中学数学其实也接触过类似的问题，当时用的最多的是用最小二乘法拟合一个一次函数并给出预测值。实际问题往往更加复杂，需要用到的拟合工具也更多。

灰度预测（主要GM11）

如果信息全部已知，那么这个系统叫做白色系统；如果一个系统信息全部未知，那么这个系统叫做黑色系统。如果一个系统信息部分已知部分未知，那么这个系统就叫做灰色系统。（取名非常形象）

灰度预测系统的特点是可以通过少量的数据进行预测，并且预测较为准确，简单常用。GM11指的是一个变量+一阶微分的系统，相对来讲是最简单的一个灰度预测（实际上我在网上还没有看到其他的更复杂维度的灰度预测模型）。所以其缺点也很明显，根据少量的数据进行预测，准确度不受到认可，在美赛中不是很常用；并且在精度上有时候不如用多项式拟合。

算法过程

①一列初始的数据 $X^{(0)}=(X^{(0)}(1),\cdots,X^{(0)}(n))$

②进行一次累加，得到前缀和序列 $X^{(1)}=(X^{(1)}(1),\cdots,X^{(1)}(n))$ ,其中 $X^{(1)}(k)=\displaystyle\sum_{i=1}^kX^{(0)}(i)$

③生成均值序列 $Z^{(1)}=(Z^{(1)}(2),\cdots,Z^{(1)}(n))$ ,其中 $Z^{(1)}(k)=\alpha X^{(1)}(k)+(1-\alpha )X^{(1)}(k-1)$ ，α的取值在(0,1)，一般取0.5。此过程也叫光滑处理。

④得到灰微分方程的基本形式 $X^{(0)}(k)+aZ^{(1)}(k)=b$ ，其中-a为发展系数，b为灰作用量。

以及得到白化微分方程 $\frac{dX^{(1)}}{dt}+aX^{(1)}(t)=b$ .

⑤将灰微分方程移项，得到 $-aZ^{(1)}(k)+b=X^{(0)}(k)$ ，于是得到

$\left$ \begin{matrix} -Z^{(1)}(2) & 1 \\ -Z^{(1)}(3) & 1 \\ \vdots & \vdots \\ -Z^{(1)}(n) & 1 \\ \end{matrix} \right$ \left$ \begin{matrix} a \\ b \\ \end{matrix} \right$ = \left$ \begin{matrix} X^{(0)}(2) \\ X^{(0)}(3) \\ \vdots \\ X^{(0)}(n) \\ \end{matrix} \right$$

即为 $X\beta=Y$ 的形式。通过最小二乘法，求出 $\widehat \beta =(X^TX)^{-1}X^TY$ ，因此求出 $\widehat a,\widehat b$ 。

⑥将求出的 $\widehat a,\widehat b$ 代入白化微分方程，并用常微分方程里的常数变易法得到：

$X^{(1)}(k)=(X^{(0)}(1)-\frac{b}{a}) e^{-\alpha(k-1)}+\frac{b}{a}$

再通过 $X^{(1)}(k)$ 的表达式求出 $X^{(0)}(k)$ 的表达式，并进行预测，即

$X^{(0)}(k)=X^{(1)}(k)-X^{(1)}(k-1)=(X^{(0)}(1)-\frac{b}{a}) e^{-\alpha(k-1)}(1-e^a))$

从表达式可以看出，这是一个一阶指数型的表达式。

例题以及代码实现

例：假定给定了2015~2021年重庆市教育投入经费如下表，运用GM(1,1)模型来预测2022年重庆市教育投入总经费。

2015~2021年重庆市教育经费总投入（单位：万元）
2015	2016	2017	2018	2019	2020	2021
6565622	6979973	7971003	8863208	9483526	10216274	11442967

matlab代码如下：

clc
clear all
close all
x0=[6565622,6979973,7971003,8863208,9483526,10216274,11442967]
%原始数据
n=length(x0);
x1=zeros(1,n);
x1(1)=x0(1);
for j=2:n
    x1(j)=x1(j-1)+x0(j); %求前缀和
end
alpha=0.5;%光滑处理的参数
z1=zeros(1,n)
z1(1)=0;
for j=2:n
    z1(j)=alpha*x1(j)+(1-alpha)*x1(j-1);
end
Y=zeros(n-1,1);X=zeros(n-1,2)
for i=2:n
    Y(i-1)=x0(i);
    X(i-1)=-z1(i);
    X(i-1,2)=1;
end
%生成对应的矩阵
B=inv(X'*X)*X'*Y;
a=B(1);
b=B(2);
x=0:0.1:7;
y3=(x0(1)-b/a).*exp(-a*x).*(1-exp(a));
plot(x,y3,'*');
hold on
n0=0:6;
plot(n0,x0,'ro')
xlabel('年份');ylabel('经费  单位:万元')

最终得到结果如图

单独提取出y3数组，得出预计到2022年重庆市教育投入总经费约为1249亿元。

插值与拟合

插值与拟合是预测模型中常用的方法，两者相似但有细微的区别。两者都是针对给定的一组数据得出一条近似的曲线（也可以是多维），插值要求这条曲线经过每一个给定的数据点，而拟合只要求大致符合，不必经过所有数据点。matlab中有特定的函数来完成数据的插值与拟合。

最近邻插值（零阶插值）

插值结果是一个间断的零阶函数，感觉应该不太会用到。

线性插值

算法思路很简单，把相邻的两个点连接起来，最终构成有多个一次函数段连起来的连续函数。显然精确度不会很高。

分段三次Hermite插值（cubic）

运用了厄密特插值法，相当于是线性插值的一个拓展。思路还是将相邻的两个点作为一对，两点间连接一条三次函数曲线，这样就比线行插值的精确度更高。其几何意义是插值后的曲线与原函数在插值点有公共切线。

样条插值（spline）

最终会得到一条二阶导连续的函数，但算法更为复杂。

插值的matlab代码实现

x=0:4:20;
y=[37 51 45 74 83 88];
xx=0:1:20;
y1=interp1(x,y,xx,'nearest');%最近邻插值
subplot(2,2,1);
plot(x,y,'o',xx,y1)
title('最近邻插值')
hold on;
y2=interp1(x,y,xx);%线性插值
subplot(2,2,2);
plot(x,y,'o',xx,y2)
title('线性插值')
hold on;
y3=interp1(x,y,xx,'pchip');%分段三次Hermite插值
subplot(2,2,3);
plot(x,y,'o',xx,y3)
title('分段三次Hermite插值')
hold on;
y4=interp1(x,y,xx,'spline');%样条插值
subplot(2,2,4);
plot(x,y,'o',xx,y4)
title('样条插值')

得出插值的图像如下：

可以粗略地看出样条插值和分段三次Hermite插值的效果更好。matlab同样能实现多维的插值，代码如下：

x=0:4:16;y=0:4:16;
z=[62 73 80 85 87;76 88 97 72 105;88 108 63 125 128;98 118 132 145 142;106 123 139 120 150];
[X,Y]=meshgrid(0:16,0:16);
Z1=interp2(x,y,z,X,Y,'nearest');
Z2=interp2(x,y,z,X,Y);
Z3=interp2(x,y,z,X,Y,'cubic');
Z4=interp2(x,y,z,X,Y,'spline');
subplot(2,2,1),surfc(X,Y,Z1),title('nearest')
subplot(2,2,2),surfc(X,Y,Z2),title('linear')
subplot(2,2,3),surfc(X,Y,Z3),title('cubic')
subplot(2,2,4),surfc(X,Y,Z4),title('spline')

效果图如下：

除了以上的插值方法，还有散乱点插值、范克利金插值等。在matlab中可以用interp3和interpn进行3维和多维的插值。

最小二乘法拟合

拟合算法也是通过已有的数据对未知的数据进行预测。在matlab中一般有多项式拟合和cftool工具拟合两种办法。

clc
clear all
x=[1 5 8 15 26 29 35 44 70];
y=[3 4 9 23 35 49 32 58 99];
p=polyfit(x,y,6);
xx=0:0.1:80;
yy=polyval(p,xx);
plot(xx,yy,x,y,'o');
title('多项式拟合')

效果图如下：