五、加权的最小二乘法
上一篇文章介绍了最小二乘估计的基本原理和估计精度，并且举了例子来说明这个最小二乘估计是如何应用的。在上一篇最后，我提了一个概念，就是最小二乘估计精度不高。为什么这么说呢？还是用一个例子来说明吧。
跟上一个例子类似，也是通过一台仪器对一个数据进行观测，可这次不同了，我用了两台仪器，而且两台仪器的精度也不同。这次简单些，每台仪器只观测一次，两台仪器获得了两个观测量：
Z
1
Z_1
Z1、
Z
2
Z_2
Z2。而这两台仪器的测量方差分别是r和4r。
没啥说的。定义各个向量和矩阵就好了。
Z
=
(
Z
1
Z
2
)
Z=\begin{pmatrix}Z_1\\Z_2\end{pmatrix}
Z=(Z1Z2)，
H
=
(
1
1
)
H=\begin{pmatrix}1\\1\end{pmatrix}
H=(11)，
R
=
(
r
0
0
4
r
)
R=\begin{pmatrix}r&0\\0&4r\end{pmatrix}
R=(r004r)。
代入公式，
X
^
=
(
H
T
H
)
−
1
H
T
Z
\hat X=(H^ TH)^{-1}H^TZ
X^=(HTH)−1HTZ得到：
X
^
=
1
2
(
Z
1
+
Z
2
)
\hat X=\frac{1}{2}(Z_1+Z_2)
X^=21(Z1+Z2)
早就猜到了，估计结果是观测的算术平均值，没毛病。可问题出在估计方差上。再把上面的已知代入估计方差的计算公式
E
[
X
~
X
~
T
]
=
(
H
T
H
)
−
1
H
T
R
H
(
H
T
H
)
−
1
E[\widetilde X\widetilde X^T]=(H^ TH)^{-1}H^TRH(H^ TH)^{-1}
E[X
X
T]=(HTH)−1HTRH(HTH)−1就会得到：
E
[
X
~
2
]
=
5
4
r
E[\widetilde X^2]=\frac{5}{4}r
E[X
2]=45r
发现问题了吗，用两台仪器对数据进行两次测量的精度还不如只用第一台仪器测量一次的精度。
究其原因，上一篇文章提过一嘴，就是一般的最小二乘估计是不分优劣的使用量测值的，从而导致估计精度不高。那有没有办法在估计的过程考虑一下量测值的质量呢？有的，就是加个权重。
还记得最小二乘估计的最优指标吧，让
(
Z
−
H
X
^
)
T
(
Z
−
H
X
^
)
(Z-H\hat X)^ T(Z-H\hat X)
(Z−HX^)T(Z−HX^)
达到最小。我们把这个最优指标加个权重系数W，变成求让
(
Z
−
H
X
^
)
T
W
(
Z
−
H
X
^
)
(Z-H\hat X)^ TW(Z-H\hat X)
(Z−HX^)TW(Z−HX^)
达到最小时的
X
^
\hat X
X^值。这样解出的最小二乘估计就会变成下面这个样子：
X
^
=
(
H
T
W
H
)
−
1
H
T
W
Z
\hat X=(H^ TWH)^{-1}H^TWZ
X^=(HTWH)−1HTWZ
可以看出，当W是单位矩阵时，上面这个跟一般的最小二乘估计是一样的。精度方面的话，加权后的估计方差：
E
[
X
~
X
~
T
]
=
(
H
T
W
H
)
−
1
H
T
W
R
W
H
(
H
T
W
H
)
−
1
E[\widetilde X\widetilde X^T]=(H^ TWH)^{-1}H^TWRWH(H^ TWH)^{-1}
E[X
X
T]=(HTWH)−1HTWRWH(HTWH)−1
好了，下面的问题就是如何选取那个W方阵了。可以证明，当
W
=
R
−
1
W=R^{-1}
W=R−1时，估计的方差最小。证明过程并不复杂，这里就不写了，浪费精力。于是就有了下面两个式子：
X
^
=
(
H
T
R
−
1
H
)
−
1
H
T
R
−
1
Z
\hat X=(H^ TR^{-1}H)^{-1}H^TR^{-1}Z
X^=(HTR−1H)−1HTR−1Z
E
[
X
~
X
~
T
]
=
(
H
T
R
−
1
H
)
−
1
E[\widetilde X\widetilde X^T]=(H^ TR^{-1}H)^{-1}
E[X
X
T]=(HTR−1H)−1
这个改进后的加权最小二乘估计还有一个名字，叫做马尔可夫估计。我之前不是吐槽说一般最小二乘估计的公式里没有体现量测数据的精度吗，这回公式里体现出来了。来吧看看估计的效果吧，还是上面那个例子：
R
=
(
r
0
0
4
r
)
R=\begin{pmatrix}r&0\\0&4r\end{pmatrix}
R=(r004r)所以
R
−
1
=
(
1
r
0
0
1
4
r
)
R^{-1}=\begin{pmatrix}\frac{1}{r}&0\\0&\frac{1}{4r}\end{pmatrix}
R−1=(r1004r1)
代入计算，就可以得到：
X
^
=
4
5
Z
1
+
1
5
Z
2
\hat X=\frac{4}{5}Z_1+\frac{1}{5}Z_2
X^=54Z1+51Z2
原来如此，精度高的仪器的测量值在估计结果中的比重大一些，精度低的仪器测量值比重小一些，合理！看一下估计精度：
E
[
X
~
2
]
=
4
5
r
E[\widetilde X^2]=\frac{4}{5}r
E[X
2]=54r
估值精度果然提高了，而且比任意一次测量的精度都高。
六、递推最小二乘估计（RLS）
不管是一般的，还是加权的最小二乘估计，有个特点，就是要想运用，就必须获得全部的量测值。如果量测数量非常巨大，比如做曲线拟合时，经常是对几千几万次量测数据进行处理，耗存储空间的事咱可以暂且不说，光是庞大的矩阵计算，其计算量可不是闹着玩的。
于是就有了递推最小二乘估计。就是把所有的量测值，拆成k次量测。其中第i次量测的方程为：
Z
i
=
H
i
X
+
V
i
Z_i=H_iX+V_i
Zi=HiX+Vi
所谓递推，是指已经知道了前k次量测后的最小二乘估计结果
X
^
k
\hat X_k
X^k，当第k+1次量测值
Z
k
+
1
Z_{k+1}
Zk+1出现时，求
X
^
k
+
1
\hat X_{k+1}
X^k+1。推导过程有点复杂，有兴趣的同学可以自己试试。这里还是直接给出结论吧。
先定义个
P
k
P_k
Pk。
P
k
=
(
H
‾
k
T
W
‾
k
H
‾
k
)
−
1
P_k=(\overline H^T_k\overline W_k\overline H_k)^{-1}
Pk=(HkTWkHk)−1
带横线的量表示前k次量测的参数集合，我们把加权最小二乘估计公式里的那个求逆的部分单独用
P
k
P_k
Pk表示了。可以看出，如果是马尔可夫估计，
P
k
P_k
Pk其实就是前k次量测后的估值结果的方差。
那么递推算法可以由下面两个公式实现：
P
k
+
1
=
P
k
−
P
k
H
k
+
1
T
(
W
k
+
1
−
1
+
H
k
+
1
P
k
H
k
+
1
T
)
−
1
H
k
+
1
P
k
P_{k+1}=P_k-P_kH^T_{k+1}(W_{k+1}^{-1}+H_{k+1}P_kH^T_{k+1})^{-1}H_{k+1}P_k
Pk+1=Pk−PkHk+1T(Wk+1−1+Hk+1PkHk+1T)−1Hk+1Pk
X
^
k
+
1
=
X
^
k
+
P
k
+
1
H
k
+
1
T
W
k
+
1
(
Z
k
+
1
−
H
k
+
1
X
^
k
)
\hat X_{k+1}=\hat X_k+P_{k+1}H^T_{k+1}W_{k+1}(Z_{k+1}-H_{k+1}\hat X_k)
X^k+1=X^k+Pk+1Hk+1TWk+1(Zk+1−Hk+1X^k)
变量上面不带横线就表示当前这次量测对应的量。这是个通用公式，不想加权，就让
W
k
W_k
Wk是单位矩阵，用马尔可夫估计，就让
W
k
=
R
k
−
1
W_k=R_k^{-1}
Wk=Rk−1。
上面那个估值公式可以描述为：前k+1次量测的估值结果是将前k次量测的估值结果根据第k+1次量测值进行修正而获得。要想用这个递推公式，需要初值，这好办，只要将最初的几次量测利用一般或者加权的最小二乘估计就可以，然后就可以一步一步的递推了。
公式抽象，还是例子更直观一些，来吧，上例子。
又一个例子
对一个变量进行直接观测，已经观测了k次了，并且已经获取了前k次的一般最小二乘估计结果 X ^ k \hat X_k X^k，求当新的量测 Z k + 1 Z_{k+1} Zk+1获得后， X ^ k + 1 \hat X_{k+1} X^k+1是多少。
一般最小二乘估计，那就不用加权了。所以
P
k
=
(
H
‾
k
T
H
‾
k
)
−
1
=
(
(
1
1
⋯
1
)
(
1
1
⋮
1
)
)
−
1
=
1
k
P_k=(\overline H^T_k\overline H_k)^{-1}=(\begin{pmatrix}1&1&\cdots&1\end{pmatrix}\begin{pmatrix}1\\1\\\vdots\\1\end{pmatrix})^{-1}=\frac{1}{k}
Pk=(HkTHk)−1=((11⋯1)⎝⎜⎜⎜⎛11⋮1⎠⎟⎟⎟⎞)−1=k1
代入那两个递推公式，就可以得到：
P
k
+
1
=
1
k
+
1
P_{k+1}=\frac{1}{k+1}
Pk+1=k+11
X
^
k
+
1
=
X
^
k
+
1
k
+
1
(
Z
k
+
1
−
X
^
k
)
\hat X_{k+1}=\hat X_k+\frac{1}{k+1}(Z_{k+1}-\hat X_k)
X^k+1=X^k+k+11(Zk+1−X^k)
这个结论非常好用，以后我们在求平均值的时候，如果已经有了前k个数的平均值，再来一个数，就不需要重新累加k+1个数再求平均，只需要在前k个数平均值上加上第k+1个数与前k个数平均值之差的k+1分之1就好了。
七、最小二乘法总结
最小二乘估计的最大优点就是算法简单，尤其是一般最小二乘估计，不必知道量测误差的统计信息。但这也是它的一种局限。
另外，最小二乘算法只能估计确定的常值向量，不能估计随时间变化的向量。
最后要说明，最小二乘估计的最优指标是量测的方差最小，这个指标其实并不能保证估计的误差最小，所以从这方面来讲，它的估计精度也是不高的。那有没有一种最优指标直接定义估计结果的某些性能达到最优呢？下一篇再讲吧。
扫一扫
1898
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
