模型估计-Charlie小站


最早对于模型估计的认识，来源于底下这组令我深恶痛绝的线性回归公式
$$\left\{
\begin{eqnarray}
\hat{b} &=& \frac{\sum_{i=1}^n(x_i-\bar{x})(y_I-\bar{y})}{\sum_{i=1}^n(x_i-\bar{x}^2)} = \frac{\sum_{i=1}^nx_iy_i-n\bar{x}\bar{y}}{\sum_{i=1}^nx_i^2-n\bar{x}^2} \\
\hat{a} &=& \bar{y}-\hat{b}\bar{x} \\
\end{eqnarray}
\right.$$
不是我记不住这个公式，而是在当时（现在也一样）落后的考试制度下，必须用计算器把一个个的数算出来。一般的考试时都会给5、6组数据，光摁计算器就要少说花两分钟。

另一个记忆深刻的是高中时期伏安法测电阻。这个的印象好多了，只要画个图，再算个斜率就得到电阻，整个求解过程满是轻松愉悦的气氛。

现在，也经常遇到的模型估计问题。比如，经常会需要做曲线拟合。虽然这类问题已经有无数现成的工具可以一步解决，但有些时候还是要自己动手，丰衣足食。下面开始正题，简单地介绍一下如何进行模型的估计。

模型估计问题看起来不难，但随便找本书都能写上百十来页。这篇文章就挑些简单的作简要的描述，难的我也不会。本文主要涉及两个方程，一个是$$\mathbf{A}\mathbf{x}=\mathbf{0}$$另一个是
$$\mathbf{A}\mathbf{x}=\mathbf{b}$$
或者换一个写法，就是
$$\arg\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x}\|_2$$或$$\arg\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x}-\mathbf{b}\|_2$$
其中，$\mathbf{A}$是矩阵，$\mathbf{x}$，$\mathbf{b}$，$\mathbf{0}$是维数与$\mathbf{A}$相适应的列向量，这里设$\mathbf{x}$是一个$n$维向量。

SVD分解

第一个方程一般用奇异值分解（SVD）求解。简单地说，就是可以把$\mathbf{A}$分解成
$$\mathbf{A} = \mathbf{U}\mathbf{D}\mathbf{V}^T$$
$\mathbf{U}$和$\mathbf{V}$都是正交矩阵（严格来说，是酉矩阵），$\mathbf{D}$为对角矩阵，而且对角线上的元素非负，从左上角到右下角递减。至于SVD为什么存在，移步维基百科。$\mathbf{V}$的最后一列就是使得$\|\mathbf{A}\mathbf{x}\|_2$取得最小值的单位向量。利用SVD的表达式，要验证这一点并不难。因为$\mathbf{V}$是正交矩阵，所以它的所有列$\{\mathbf{v}_k\}$构成$\mathbf{x}$所在空间的一组标准正交基，于是$\mathbf{x} = \sum_{k}a_k\mathbf{v}_k$。设$\mathbf{D}=\mathrm{diag}\{\lambda_k\}$，利用SVD表达式计算$\mathbf{Ax}$的范数的平方
$$\begin{eqnarray}
\|\mathbf{Ax}\|_2^2 &=& \|\mathbf{UDV}^T\sum_ka_k\mathbf{v}_k\|_2^2 \\
&=& \|\mathbf{DV}^T\sum_ka_k\mathbf{v}_k\|_2^2 \quad\quad （\mathbf{U}\text{具有保范性}）\\
&=& \|\mathbf{D}\sum_ka_k\mathbf{V}^T\mathbf{v}_k\|_2^2 \quad\quad （\text{利用}\mathbf{V}的正交性）\\
&=& \sum_k|\lambda_ka_k|_2^2 \quad\quad （\text{利用}\mathbf{V}的正交性和\mathbf{D}的对角性）\\
&\geq& \sum_k|\lambda_na_k|_2^2 \quad\quad （\lambda_{k+1}\geq\lambda_k）\\
&=& \lambda_n^2 \quad\quad （\mathbf{x}为单位向量）
\end{eqnarray}$$
上式说明在任何情况下，$\mathbf{Ax}$的范数都不可能小于$\lambda_n$。而当$a_n = 1$时，$\|\mathbf{Ax}\|_2=\lambda_n$。

最小二乘

一般称$$\arg\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x}-\mathbf{b}\|_2$$的解为最小二乘解，其结果为$$\mathbf{x}^* = (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{b}$$这个解可以直接通过对代价函数（即范数方程）求导，令其导数为零求得。另一个方案是利用几何结论。由几何关系，最优解应该满足$$\mathbf{A}^T(\mathbf{Ax}-\mathbf{b})=\mathbf{0}$$。这是一个很多人都知道的结论，但是有很多人反应不过来。

以三维空间为例，这个问题用人话描述就是已知空间中的一个点A，现要在平面上找到一个与A最近的点，怎么找？就像上图那样，过A点做平面的垂线，交点B就是满足条件的点。所谓与平面垂直，是指AB与平面上的所有向量都垂直，即积为零。由于这个平面是一个向量空间，所以只要AB与空间的一组基底正交即可。

这里的$\mathbf{Ax}$是$\mathbf{A}$的列向量所张成的空间中的一个向量。$\mathbf{Ax}-\mathbf{b}$就相当于三维空间中的AB。AB与空间中的所有向量正交，只需$\mathbf{Ax}-\mathbf{b}$与$\mathbf{A}$的每一列都正交，因此$$\mathbf{A}^T(\mathbf{Ax}-\mathbf{b})=\mathbf{0}$$从而可以得到最小二乘解为$$\mathbf{x}^* = (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{b}$$

参数估计举例

接下来以最为简单的直线拟合来说明如何通过SVD分解或者是直接求最小二乘解进行模型的估计。我通过MATLAB生成了一组数据

x = [1:5];
y = 3*x + 4 + randn(1,5);

构造的数据如下表所示

$x$	1	2	3	4	5
$y$	6.7950	9.87592	14.4897	17.4090	20.4172

现在要通过模型构造出方程$\mathbf{A}\mathbf{x}=\mathbf{0}$或者是$\mathbf{A}\mathbf{x}=\mathbf{b}$。设直线方程为$$kx+d-y=0$$，写成矩阵形式为$$[x~~1~~-y][k~~d~~1]^T = 0$$或者是$$[x~~1][k~~d]^T = y$$
首先，通过SVD分解求解模型，利用MATLAB求解

n = length(x);
A = zeros(n,3);
for i=1:length(x)
    A(i,:) = [x(i) 1 -y(i)];
end
[~,~,v] = svd(A);
v_end = v(:,end);
v_end = v_end/v_end(3);
k = v_end(1)
d = v_end(2)

再用最小二乘模型，利用MATLAB求解

n = length(x);
A = zeros(n,2);
b = zeros(n,1);
for i=1:length(x)
    A(i,:) = [x(i) 1];
    b(i,:) = [y(i)];
end
result = (A'*A)\A'*b;
k = result(1)
d = result(2)

使用两种方法估计出的模型为