谈谈序列的一种拟合方法

问题提出:

  有递推式: $$y_{t} = 0.2 + 0.3y_{t-7}+e\,\,,\,\, e\sim N(0,1),t为正整数$$
  有初始化序列:$$[y_{1},y_{2},\cdots,y_{7} ]=[10,20,30,40,50,60,70]$$
  根据目标函数生成了100个数据 $$y_{8},y_{9},\cdots,y_{107}$$

  如图:

  求:假设我们不知到目标函数的具体参数,但我们知道模型是\(y_{t} = a + b\cdot y_{t-7}\),并且有数据\(y_{1},y_{2},\cdots,y_{107}\)(即上图)
  请解出\(a,b\)的值。

分析

  题目的意思是已知函数模型以及数据集,求解模型里的参数。
  首先考虑用怎么度量解(a,b)的好坏,这里选择平滑好用的度量方法—-欧氏距离。以下给出具体的代价函数:
$$ cost\,\,\, function : J(a,b) = \frac{1}{100}\cdot \sum_{i=8}^{107}((a+b\cdot y_{i-7}) - y_{i})^{2}$$
  显然,模型是一个7阶的差分方程,既然是差分方程的话,是否可以考虑用\(Z\)变换去做?可否在二维空间里用群粒子算法或者遗传算法求出解(a,b)?其实,对于这总常系数差分序列的拟合,有一种简单、巧妙的处理方法—-转化为线性回归问题(least-squares programming 问题)。

破题


  由序列递推事可得:
$$ y_{8} = a + b\cdot y_{1} \\
y_{9} = a + b\cdot y_{2} \\
y_{10} = a + b\cdot y_{3} \\
\vdots \\
y_{107} = a + b\cdot y_{100}
$$
  显然这里有一个参数矩阵。
  把每一笔数据\(y_{i}\)用2个维度的向量表示,而不是一个数值。
$$ X = \begin{bmatrix}
1 & y_{1} \\
1 & y_{2} \\
1 & y_{3} \\
\vdots \\
1 & y_{100}
\end{bmatrix} \,\,,\,X \in \mathbb{R}^{100 \times 2}$$
注意,第一笔数据特征是\([1,y_{1}]\),标签值是$y_{8}$

$$Y_p = XW = \begin{bmatrix}
1 &y_{1} \\
1 & y_{2} \\
1 & y_{3} \\
\vdots \\
1 & y_{100}
\end{bmatrix} \cdot
\begin{bmatrix}
a \\
b \\
\end{bmatrix}
$$

$$ Y = \begin{bmatrix}
y_{8} \\
y_{9} \\
\vdots \\
y_{107}
\end{bmatrix}
$$
  到这里,问题已经转化为如下的least-squares programming 问题了。
$$ \min_{W}\,\,\, J(w_{0},w_{1})=J(a,b) = \frac{1}{100}\cdot \sum_{i=8}^{107}((a+b\cdot y_{i-7}) - y_{i})^{2} $$
  有解析解(后面有具体推导):$$W = X^{\dagger} Y = (X^{T}X)^{-1}X^{T}Y$$
  由此直接算出\(W\),即\( a,b\)的值。
  拟合出来的效果如下所示:

1
2
3
4
5
6
W =
0.116573700765463
0.295857735126161
cost =
1.042281571564902

更一般的做法

  假设现在只知道目标函数是7阶的常系数差分方程,而不知道具体的函数模型,如何求解?
  方法同上,只是需要做一点点修改。
  把每一笔数据\(y_{i}\)用8个维度的向量表示成\(X_{qi}\)。
$$
X_{q(i-7)} = [ 1,y_{i-7},y_{i-6},y_{i-5},y_{i-4},y_{i-3},y_{i-2},y_{i-1} ] \\
i=8,9,10,\cdots,m+7
$$
  即拟合模型是:
$$
y_{t}= a+ b_{1}y_{t-7} + b_{2}y_{t-6}+ b_{3}y_{t-5}+ b_{4}y_{t-4}+ b_{5}y_{t-3}+ b_{6}y_{t-2}+ b_{7}y_{t-1}
$$
没错,相信已经有读者知道这是自回归。
  同理求出 $$W_{q} = X_{q}^{\dagger} Y = (X_{q}^{T}X_{q})^{-1}X_{q}^{T}Y\,\,\,,\,W_{q} \in \mathbb{R^8}$$
  拟合出来的效果如下图所示:

1
2
3
4
5
6
7
8
9
10
11
12
Wq =
0.254232073618409
0.289206190751554
-0.000633472123790
0.000018960275785
0.008956957040519
-0.019245627029181
0.002814318547785
-0.003671158221363
cost =
1.068767118660750


  相信细心的读者已经注意到了\(W_q\)向量的前面两个数值。
  这种序列拟合的方法到这里已经介绍完毕了,下面是一些进一步的分析,读者去留随意。Matlab代码下载点击这里

噪声,噪声,不死的噪声

  我们知道,在采集数据的时候,不死的噪声总是无处不在,它们获得或多或少都会给模型的最终表现带来消极的影响。我们把初始化序列变为原来的0.1倍(试图增加噪声的影响力),
  拟合的效果如图:


  我们可以看到,单从数据点上看的话,应该说情况太糟糕了,可又发现,即使这样拟合出来的曲线却没有那么糟糕,为什么呢?

  这里的噪声\(e\sim N(0,1)\),是正态分布,均值为0。所以,也许(笔者暂时未能给出论证)对于线性回归模型(linear regression)来说该噪声是比较友好的,不会太捣蛋。

  最后有一个残忍的事实,尽管这个时候噪声看起来不捣蛋,但此时拟合出来的曲线似乎没多大意义了(t越大,预测出来的准确度严重受到噪声的影响。),从这个角度看,噪声还是捣蛋了%>_<%。


小结:

  1. 这里如何消减噪声的带来的消极影响暂时还没找到很好办法。
  2. 对于序列的处理方法有很多,比如移动平均模型(MA)、自回归移动平均模型(ARMA)以及差分自回归移动平均模型(ARIMA)
  3. 提出一个新问题:把问题进一步一般化,假设我们拿到数据集并且只知道是常系数差分递推序列(不知道多少阶),如何求解?

线性回归(Linear regression)的解析解推导

  台湾大学林老师的线上课程《机器学习基石》-第九讲有很明了的介绍,这里直接引用。
  首先,我们表达式写成矩阵的形式。矩阵\(X\)的每一行是一笔数据,每一列是一个特征值。列向量\(W\)是权重参数向量。

  得到

  我们希望能找到一个\(W\)使得这个cost function 在各个方向上的梯度(偏导)都为0

  为了方便求导,我们代换一下

  求导

  得到

  问题简化成如下

  推出结论

注意,这里的\(X\)是可逆的(即满秩)。当\(X\)是奇异矩阵(非满秩)的情况下就不太一样,这里不做讨论