在信息处理问题中，观测到的有用信号总受到观测噪声的干扰。从带有观测噪声的数据中得到所需要的各种参量的估计值，这就是估计问题。为了衡量估计的质量，必须要有一个估计准则。估计准则是被估计参数的损失函数或目标函数，任何一种最优估计都是满足损失函数最小或目标函数最大的估计，估计准则不同，得到的估计也不同。常用的估计准则有：最小二乘、最小方差、极大似然、极大验后和贝叶斯风险最小。

最小二乘是以拟合误差为自变量定义的损失函数，最小二乘估计是对损失函数极小化推导而得到的，它适用于对参数的统计规律未知的情况；最小方差估计以参数估计误差的二阶矩为损失函数，它需要已知观测值和参数有关的矩；极大似然估计和极大验后估计以参数的概率密度为目标函数，并对目标函数极大为条件导出，因此它需要更多的先验统计信息。

通常在估计前，要确立观测值与参数的数学关系，以建立数学模型。如果实施的是线性估计，还需要对其中的函数模型进行线性化。为此，本章首先介绍如何建立参数估计的数学模型和模型的线性化过程，然后推导最小二乘估计、最小方差估计、极大似然估计和极大验后估计，最后介绍贝叶斯估计，并分析贝叶斯估计与其他估计方法之间的联系和各自的特性。

与第 3 章开始介绍的“状态”估计方法不同，本章介绍的估计方法不考虑被估计对象的变化或运动规律，因此在数学模型中没有描述被估计对象变化规律的状态方程，对当前时刻来说，被估计量是“静态”的，所以也被称为“参数估计”。

参数估计问题的数学模型

数学模型就是用数学的语言，如用变量、方程和不等式等来描述研究对象的特征及其各个变量内在联系，建立数学模型是实现估计的第一步。通常情况下，数学模型中的函数模型是非线性的，非线性函数模型的估计需要采用复杂的优化算法来求解，这给解决现实问题带来困难。在应用中，通常将非线性函数转化为线性函数，转化后的线性函数虽然是原来非线性函数模型的近似，给估计带来一定程度上的损失，但极大地简化了计算，使估计更容易实现。

数学模型的建立

在现实中，我们感兴趣的对象大多是不可以直接量测的，如 GNSS 卫星导航量测的是卫星到用户的距离，而不是用户的位置，因此，首先要建立观测值与用户位置的必然关系。我们用观测方程来描述观测值与被估计参数之间的必然关系，这个必然关系也称为函数模型。由于观测值是受到随机误差（噪声）干扰的随机变量，所以数学模型既要考虑观测值与待估参数的确定性关系，也要考虑观测值和待估参数的不确定性。在理想的情况下，用概率分布来描述随机误差；在分布未知的情况下，用随机变量的特征值，如方差来度量模型的不确定性，这些对观测噪声随机特性的描述称为随机模型。

在建立数学模型时，我们力求能够真实、系统和完整地反映现实问题，但在建模过程中，模型不易过于复杂而难以计算。所以在确保模型一定准确性的条件下，可以忽略那些非本质的、对客观真实程度影响不大的部分，从而使数学模型更加简明实用。

1. 观测方程和随机模型

以下面的例子来说明如何建立观测方程和随机模型。

例 2.1如图 2.1 所示，某一质点沿着直线做匀速运动，其轨迹为图中的实体直线。质点的纵坐标与质点运动的速度 \(\beta\) 和在初始时刻（\(t_0=0\)）的位置 \(\alpha\) 可以描述为 \[\widetilde{Z}=\alpha+t\beta \tag{2.1.1}\] 为了估计 \(\beta\) 和 \(\alpha\)，现在不同时刻 \(t_1,\ t_2,\ \cdots,\ t_6\) 不等精度地观测了质点的纵坐标，观测值为 \(Z_1,\ Z_2,\ \cdots,\ Z_6\)（图中的圆点），且各观测值间随机独立。各观测值和观测中误差见表 2.1，请建立观测值与待估计参数 \(\beta\) 和 \(\alpha\) 的数学关系。

若设参数为 \(\bm{X}=\left[\begin{array}{ll}\alpha & \beta\end{array}\right]^{\mathrm{T}}\)，在时刻 \(t_1,\ t_2,\ \cdots,\ t_{\ell}\) 质点的实际坐标为 \(\widetilde{\bm{Z}}=\left[\begin{array}{llll}\widetilde{Z}_1, & \widetilde{Z}_2, & \cdots & \widetilde{Z}_{\ell}\end{array}\right]^{\mathrm{T}}\)，其中 \(\ell\) 是观测值的个数，这里 \(\ell=6\)。\(\widetilde{Z}_i\ (i=1,\ 2,\ \cdots,\ \ell)\) 与参数的关系表达为： \[\begin{cases} \widetilde{Z}_1=\alpha+t_1\beta\\ \cdots\\ \widetilde{Z}_i=\alpha+t_i\beta\\ \cdots\\ \widetilde{Z}_{\ell}=\alpha+t_{\ell}\beta \end{cases} \tag{2.1.2}\]

令 \[\bm{H}=\begin{bmatrix}1 & t_1\\ 1 & t_2\\ \vdots & \vdots\\ 1 & t_{\ell}\end{bmatrix} \tag{2.1.3}\] 那么 \[\widetilde{\bm{Z}}=\bm{H}\bm{X} \tag{2.1.4}\] 矩阵 \(\bm{H}\) 反映了 \(\widetilde{\bm{Z}}\) 与参数之间的关系，也称为设计矩阵。参数 \(\bm{X}\) 与 \(\widetilde{\bm{Z}}\) 有如上式的确定的函数关系，在这样确定的函数关系下，只需要知道在两个不同时刻 \(t_1\) 和 \(t_2\) 的观测 \(Z_1\) 和 \(Z_2\) 后，就可以解出参数 \(\alpha\) 和 \(\beta\)。所以，对于此问题的必要观测值数为 \(n=2\)。由于在对 \(\widetilde{Z}_i\) 进行观测时，不可避免地受到观测误差 \(\Delta_i\) 的干扰，所以观测值 \(Z_i\) 为 \[Z_i=\widetilde{Z}_i+\Delta_i \tag{2.1.5}\] 那么，观测值 \(Z_i\) 与参数的关系为 \[Z_i=\alpha+t_i\beta+\Delta_i\quad (i=1,\ 2,\ \cdots,\ \ell) \tag{2.1.6}\] 上式中，观测值 \(Z_i\) 是待估计参数的函数，这样的方程称为观测方程。设 \[\bm{Z}=\begin{bmatrix}Z_1\\ Z_2\\ \vdots\\ Z_{\ell}\end{bmatrix}\ ,\quad \bm{\Delta}=\begin{bmatrix}\Delta_1\\ \Delta_2\\ \vdots\\ \Delta_{\ell}\end{bmatrix} \tag{2.1.7}\] 观测方程可以表示为： \[\underset{\ell\times 1}{\bm{Z}}=\underset{\ell\times n}{\bm{H}}\ \underset{n\times 1}{\bm{X}}+\underset{\ell\times 1}{\bm{\Delta}} \tag{2.1.8}\] 如果将 \(\bm{X}\) 视为有用信号，\(\bm{\Delta}\) 就是对信号的干扰部分，估计问题就是如何将 \(\bm{Z}\) 中的有用信号部分提取出来，从而求得待估参数 \(\bm{X}\)。

在列观测方程时，通常要求 \(\bm{H}\) 为列满秩矩阵，即 \[\mathrm{rank}(\bm{H})=n \tag{2.1.9}\]

观测误差 \(\bm{\Delta}\) 的随机特性可由统计特征值给出。通常假设随机误差 \(\bm{\Delta}\) 的期望为零，即 \[E(\bm{\Delta})=\underset{\ell\times 1}{\bm{0}}=\begin{bmatrix}0\\ 0\\ \vdots\\ 0\end{bmatrix} \tag{2.1.10}\] 方差为 \[\mathrm{Var}(\bm{\Delta})=\underset{\ell\times\ell}{\bm{D}}=\begin{bmatrix} \sigma_{z_1}^2 & \sigma_{z_1z_2} & \cdots & \sigma_{z_1z_{\ell}}\\ & \sigma_{z_2}^2 & \cdots & \sigma_{z_2z_{\ell}}\\ \text{symmetric} & & & \vdots\\ & & & \sigma_{z_{\ell}}^2 \end{bmatrix} \tag{2.1.11}\]

在式 (2.1.8) 中，如果不考虑 \(\bm{X}\) 的随机特性或者先验随机特性未知，即认为 \(\bm{X}\) 为非随机量，那么 \(\bm{H}\bm{X}\) 为确定的非随机部分，观测值 \(\bm{Z}\) 的随机特性就由观测误差决定，所以 \[\begin{aligned} E(\bm{Z})&=\bm{H}\bm{X}\\ \mathrm{Var}(\bm{Z})&=\bm{D} \end{aligned} \tag{2.1.12}\]

观测方程 (2.1.8) 是根据物理现实或者几何条件建立起的观测值与待估参数之间确定的函数关系，称为函数模型；式 (2.1.12) 给出了观测值的期望、观测值的精度和误差之间的相关性，它描述的是观测值的随机特性，称为随机模型。式 (2.1.8) 和式 (2.1.12) 一起给出了观测值与参数的关系，称为估计问题的数学模型，这样的数学模型也称为高斯-马尔可夫模型。

观察式 (2.1.8)，未知量有参数 \(\bm{X}\) 和观测误差 \(\bm{\Delta}\)，共有 \((\ell+n)\) 个。若将式 (2.1.8) 表示为线性方程组 \[\begin{bmatrix}\underset{\ell\times n}{\bm{H}} & \underset{\ell\times\ell}{\bm{I}}\end{bmatrix} \begin{bmatrix}\underset{n\times 1}{\bm{X}}\\ \underset{\ell\times 1}{\bm{\Delta}}\end{bmatrix}=\underset{\ell\times 1}{\bm{Z}} \tag{2.1.13}\] 系数矩阵为行满秩矩阵：\(\mathrm{rank}\begin{bmatrix}\bm{H} & \bm{I}\\ \scriptstyle\ell\times n & \scriptstyle\ell\times\ell\end{bmatrix}=\ell\)。由于系数矩阵增广矩阵的秩也为 \(\ell\)，未知量的个数多于方程个数，所以式 (2.1.13) 有无穷多组解。如何在这无穷多组解中选取“最优”的解就是最优估计所要解决的问题。

根据表 2.1 中的观测值，可以依次给出例 2.1 的观测方程和随机模型。

观测方程为 \[\begin{cases} 4.16=\alpha+1\beta+\Delta_1\\ 4.52=\alpha+2\beta+\Delta_2\\ \cdots\\ 9.26=\alpha+6\beta+\Delta_6 \end{cases} \tag{2.1.14}\]

随机模型为： \[E\begin{bmatrix}\Delta_1\\ \Delta_2\\ \vdots\\ \Delta_{\ell}\end{bmatrix}=\begin{bmatrix}0\\ 0\\ \vdots\\ 0\end{bmatrix}\ ,\quad \bm{D}=\begin{bmatrix}0.5^2 & & & & &\\ & 0.4^2 & & & &\\ & & 0.4^2 & & &\\ & & & 0.5^2 & &\\ & & & & 0.5^2 &\\ & & & & & 0.5^2\end{bmatrix}\ \mathrm{m}^2 \tag{2.1.15}\] 由于观测值相互随机独立，所以方差 \(\bm{D}\) 矩阵即为对角矩阵。

例 2.2在 GPS 定位中，设 GPS 信号发送时刻可见卫星的坐标为（WGS84）\((X^{s_i},Y^{s_i},Z^{s_i})\)（\(i=1,2,\cdots,\ell\)）。为了得到 GPS 接收机在接收信号时的位置 \((X_r,\ Y_r,\ Z_r)\)，观测了接收机与每颗卫星的距离 \(\bm{Z}=[\rho_1\ \rho_2\ \cdots\ \rho_{\ell}]^{\mathrm{T}}\)（假设观测量已经根据经验模型进行了卫星钟差和传播路径中的系统误差改正）。请建立观测值与待估计参数 \((X_r,\ Y_r,\ Z_r)\) 的数学关系。

由几何知识可知，距离交汇可以得到 GPS 接收机的位置。现在取其中的观测值 \(\rho_i\) 来建立与 \((X_r,\ Y_r,\ Z_r)\) 的函数关系： \[\rho_i=\sqrt{\left(X^{s_i}-X_r\right)^2+\left(Y^{s_i}-Y_r\right)^2+\left(Z^{s_i}-Z_r\right)^2}+\Delta_i\ ,\quad i=1,\ 2,\ \cdots,\ \ell \tag{2.1.16}\] 式中，\(\Delta_i\) 为观测值 \(\rho_i\) 的随机误差。在 GPS 观测中，要求接收机钟与卫星钟同步，但实际上接收机钟的稳定性较差，无法做到导航系统时间同步，所以在建模时需要对这部分系统误差进行补偿。由于接收机钟差造成的测距误差对所有卫星观测值是一样的，所以这里用参数 \(\tau\) 来吸收接收机钟差造成的测距误差：\(\tau=c\cdot\Delta t\)，\(c\) 为信号在真空中传播的速度，\(\Delta t\) 为接收机钟与系统时间不同步的误差，单位为秒。因此，观测方程为 \[\rho_i=\sqrt{\left(X^{s_i}-X_r\right)^2+\left(Y^{s_i}-Y_r\right)^2+\left(Z^{s_i}-Z_r\right)^2}+\tau+\Delta_i\ ,\quad i=1,\ 2,\ \cdots,\ \ell \tag{2.1.17}\] 需要估计的参数为 \(\bm{X}=\left[\begin{array}{llll}X_r & Y_r & Z_r & \tau\end{array}\right]^{\mathrm{T}}\)，这四个参数互不相关。若令 \[f_i(\bm{X})=\sqrt{\left(X^{s_i}-X_r\right)^2+\left(Y^{s_i}-Y_r\right)^2+\left(Z^{s_i}-Z_r\right)^2}+\tau \tag{2.1.18}\] 式 (2.1.17) 为： \[\rho_i=f_i(\bm{X})+\Delta_i \tag{2.1.19}\] 令 \[\bm{Z}=\begin{bmatrix}\rho_1\\ \rho_2\\ \vdots\\ \rho_{\ell}\end{bmatrix}\ ,\quad \bm{\Delta}=\begin{bmatrix}\Delta_1\\ \Delta_2\\ \vdots\\ \Delta_{\ell}\end{bmatrix}\ ,\quad \bm{F}(\bm{X})=\begin{bmatrix}f_1(\bm{X})\\ f_2(\bm{X})\\ \vdots\\ f_{\ell}(\bm{X})\end{bmatrix} \tag{2.1.20}\] 那么，观测方程为 \[\bm{Z}=\bm{F}(\bm{X})+\bm{\Delta} \tag{2.1.21}\] 观测值的随机特性为： \[\begin{aligned} E(\bm{Z})&=\bm{F}(\bm{X})\\ \mathrm{Var}(\bm{Z})&=\bm{D} \end{aligned} \tag{2.1.22}\] \(\bm{D}\) 为观测误差 \(\bm{\Delta}\) 方差矩阵，由于 \(\bm{F}(\bm{X})\) 为非随机量，所以 \(\bm{D}\) 也是观测值的方差矩阵。如果观测值相互随机独立，\(\bm{D}\) 为对角矩阵。在 GPS 观测中，观测值的方差可以通过与卫星的高度角或者信噪比等相关的经验公式得到。式 (2.1.21) 的函数描述了观测量与待估计参数间的关系，随机模型 (2.1.22) 描述了观测值间的相关关系和不确定性，它们一起构成了估计参数的数学模型。

2. 参数的约束条件

在上面的举例中，待估计参数之间没有联系，也就是在没有观测值前，这些参数没有确定的函数关系。但在有些情况下，要求强制参数估计满足某种条件，这时应该将这个条件描述出来，在估计时与观测方程一并考虑。

例 2.3如图 2.3 所示，已知基站 \(A_1\)，\(A_2\) 和 \(A_3\) 坐标，为了得到目标 \(P_1\) 和 \(P_2\) 的平面坐标，在这三个基站上分别观测了基站与 \(P_1\) 和 \(P_2\) 的距离，各观测值随机独立，观测误差均为 \(0.05\,\mathrm{m}\)，基站已知坐标和观测值见表 2.2。此外已经用高精度仪器观测得到 \(P_1\) 和 \(P_2\) 的距离为 \(80.50\,\mathrm{m}\)。设目标 \(P_1\) 和 \(P_2\) 的平面坐标为参数 \(\bm{X}=\left[\begin{array}{llll}X_{P_1} & Y_{P_1} & X_{P_2} & Y_{P_2}\end{array}\right]^{\mathrm{T}}\)，给出观测值与参数的函数关系，并给出参数所应该满足的约束条件。

解：将基站 \(i\) 的已知坐标表示为 \(\left(\begin{array}{ll}X^i & Y^i\end{array}\right)\)；设目标为 \(P_j\)，并设基站 \(i\) 与目标 \(P_j\) 之间的距离观测值为 \(Z_{ij}\)，那么观测值 \(Z_{ij}\) 可用参数和已知基站的坐标表示 \[Z_{ij}=\left(\sqrt{\left(X^i-X_{P_j}\right)^2+\left(Y^i-Y_{P_j}\right)^2}\ \right)+\Delta_{ij}\] 已经用高精度仪器得到了目标 \(P_1\) 和 \(P_2\) 的距离为 \(80.50\,\mathrm{m}\)，那么 \(P_1\) 和 \(P_2\) 两点的坐标参数需要满足条件 \[\sqrt{\left(X_{P_1}-X_{P_2}\right)^2+\left(Y_{P_1}-Y_{P_2}\right)^2}-80.50=0\] 将这个条件表示为更一般的形式 \[\underset{c\times 1}{\bm{\varphi}}(\bm{X})=0 \tag{2.1.23}\] 在上式中 \(c\) 为约束条件方程的个数。式 (2.1.23) 给出了参数应该满足的函数关系，在这个方程中，没有观测值，它表达的是参数之间的函数关系。此时，完整地描述此问题的函数模型为： \[\begin{cases} \bm{Z}=\bm{F}(\bm{X})+\bm{\Delta}\\ \bm{\varphi}(\bm{X})=0 \end{cases} \tag{2.1.24}\] 由于观测值随机独立，中误差为 \(5\,\mathrm{cm}\)，所以随机模型为 \[\bm{D}=\begin{bmatrix}5^2 & & & &\\ & 5^2 & & &\\ & & 5^2 & &\\ & & & 5^2 &\\ & & & & 5^2\end{bmatrix}\mathrm{cm}^2 \tag{2.1.25}\] 式 (2.1.24) 和式 (2.1.25) 一起构成了描述此问题的数学模型。在这样的模型上对参数进行估计，得到的参数估计值一定满足式 (2.1.23) 的约束条件。

函数模型的线性化

1. 函数模型的线性化

在例 2.1 的观测方程式 (2.1.8) 中，参数是观测值的线性函数，但在更一般的情况下，观测量是参数的非线性函数，如例 2.2 的观测方程就是非线性的，例 2.3 中的观测方程和限制条件都是非线性的，接下来对非线性模型进行线性化。

函数模型的线性化通常采用泰勒级数将函数展开，舍弃高阶项后得到。设 \(\bm{X}\) 的近似值为 \(\bm{X}^{*}\)，将 \(f_i(\bm{X})\) 在 \(\bm{X}^{*}\) 处展开： \[\begin{aligned} f_i(\bm{X})=f_i\left(\begin{array}{lllllll}X_1^{*} & \cdots & X_j^{*} & \cdots & X_n^{*}\end{array}\right) &+\left(\frac{\partial f_i}{\partial X_1}\right)_{*}\left(X_1-X_1^{*}\right)+\cdots\\ +\left(\frac{\partial f_i}{\partial X_j}\right)_{*}\left(X_j-X_j^{*}\right) &+\cdots+\left(\frac{\partial f_i}{\partial X_n}\right)_{*}\left(X_n-X_n^{*}\right)+h\left(X_n-X_n^{*}\right) \end{aligned} \tag{2.1.26}\] 其中 \(\left(\dfrac{\partial f_i}{\partial X_j}\right)_{*}\) 为 \(X_j\) 对函数 \(f_i\) 的一阶导数，然后代入 \(\bm{X}^{*}\) 后的值；\(h\left(X_n-X_n^{*}\right)\) 为高阶项。设 \(x_j=X_j-X_j^{*}\)，并忽略高阶项，上式为 \[\begin{aligned} f_i(\bm{X})=\left(\frac{\partial f_i}{\partial X_1}\right)_{*}x_1\cdots+ &\left(\frac{\partial f_i}{\partial X_j}\right)_{*}x_j\cdots+ \left(\frac{\partial f_i}{\partial X_n}\right)_{*}x_n\ +\\ &f_i\left(\begin{array}{lllllll}X_1^{*} & \cdots & X_j^{*} & \cdots & X_n^{*}\end{array}\right)+\Delta_i \end{aligned} \tag{2.1.27}\] 若记： \[\bm{h}_i=\left[\begin{array}{llll}h_{i1} & h_{i2} & \cdots & h_{in}\end{array}\right] =\left[\begin{array}{llll} \left(\dfrac{\partial f_i}{\partial X_1}\right)_{*} & \left(\dfrac{\partial f_i}{\partial X_2}\right)_{*} & \cdots & \left(\dfrac{\partial f_i}{\partial X_n}\right)_{*} \end{array}\right] \tag{2.1.28}\] 和 \[f_i(\bm{X}^{*})=f_i\left(\begin{array}{llll}X_1^{*} & X_2^{*} & \cdots & X_n^{*}\end{array}\right) \tag{2.1.29}\] 以及 \[\bm{x}=\bm{X}-\bm{X}^{*}=\begin{bmatrix}x_1\\ x_2\\ \vdots\\ x_n\end{bmatrix} =\begin{bmatrix}X_1-X_1^{*}\\ X_2-X_2^{*}\\ \vdots\\ X_n-X_n^{*}\end{bmatrix} \tag{2.1.30}\] 那么，观测方程为： \[Z_i=\bm{h}_i\bm{x}+f_i(\bm{X}^{*})+\Delta_i \tag{2.1.31}\] 从式 (2.1.29) 可以看出，\(f_i(\bm{X}^{*})\) 是由参数的初始值 \(\bm{X}^{*}\) 计算得到，也可以看作近似观测值 \[Z_i^{*}=f_i(\bm{X}^{*})=f_i\left(\begin{array}{llll}X_1^{*} & X_2^{*} & \cdots & X_n^{*}\end{array}\right) \tag{2.1.32}\] 将其从观测值 \(Z_i\) 中减去，得到： \[z_i=Z_i-f_i(\bm{X}^{*})=\bm{h}_i\bm{x}+\Delta_i \tag{2.1.33}\] 令：

设 \[\bm{z}=\begin{bmatrix}z_1\\ z_2\\ \vdots\\ z_{\ell}\end{bmatrix}\ ,\quad \bm{H}=\begin{bmatrix}\bm{h}_1\\ \bm{h}_2\\ \vdots\\ \bm{h}_{\ell}\end{bmatrix}\ ,\quad \bm{F}(\bm{X}^{*})=\begin{bmatrix}f_1(\bm{X}^{*})\\ f_2(\bm{X}^{*})\\ \vdots\\ f_{\ell}(\bm{X}^{*})\end{bmatrix} \tag{2.1.34}\] 最后 \[\bm{z}=\bm{Z}-\bm{F}(\bm{X}^{*})=\bm{H}\bm{x}+\bm{\Delta} \tag{2.1.35}\] 观察式 (2.1.35)，它与线性观测方程 (2.1.8) 完全一样，这样就将非线性观测方程转化为了线性方程。

现在按照以上线性化方法将式 (2.1.17) 在 \(\bm{X}^{*}=\left[\begin{array}{llll}X^{*} & Y^{*} & Z^{*} & \tau^{*}\end{array}\right]\) 处线性化，得到 \[\Delta\rho_i=\rho_i-\rho_i^{*}=\frac{-\Delta X_i^{*}}{S_i^{*}}x+\frac{-\Delta Y_i^{*}}{S_i^{*}}y+\frac{-\Delta Z_i^{*}}{S_i^{*}}z+\Delta\tau+\Delta_i \tag{2.1.36}\] 其中 \[\begin{cases} S_i^{*}=\sqrt{\left(X^{s_i}-X^{*}\right)^2+\left(Y^{s_i}-Y^{*}\right)^2+\left(Z^{s_i}-Z^{*}\right)^2}\\ \rho_i^{*}=\sqrt{\left(X^{s_i}-X^{*}\right)^2+\left(Y^{s_i}-Y^{*}\right)^2+\left(Z^{s_i}-Z^{*}\right)^2}+\tau^{*}\\ \Delta X_i^{*}=X^{s_i}-X^{*},\ \Delta Y_i^{*}=Y^{s_i}-Y^{*},\ \Delta Z_i^{*}=Z^{s_i}-Z^{*}\\ x=X-X^{*},\ y=Y-Y^{*},\ z=Z-Z^{*},\ \Delta\tau=\tau-\tau^{*} \end{cases} \tag{2.1.37}\] 例 2.2 的观测方程写成矩阵形式为： \[\begin{bmatrix}\Delta\rho_1\\ \Delta\rho_2\\ \vdots\\ \Delta\rho_{\ell}\end{bmatrix} =\begin{bmatrix}\rho_1-\rho_1^{*}\\ \rho_2-\rho_2^{*}\\ \vdots\\ \rho_{\ell}-\rho_{\ell}^{*}\end{bmatrix} =\begin{bmatrix} -\dfrac{\Delta X_1^{*}}{S_1^{*}} & -\dfrac{\Delta Y_1^{*}}{S_1^{*}} & -\dfrac{\Delta Z_1^{*}}{S_1^{*}} & 1\\[10pt] -\dfrac{\Delta X_2^{*}}{S_2^{*}} & -\dfrac{\Delta Y_2^{*}}{S_2^{*}} & -\dfrac{\Delta Z_2^{*}}{S_2^{*}} & 1\\[10pt] \vdots & \vdots & \vdots & \vdots\\[4pt] -\dfrac{\Delta X_{\ell}^{*}}{S_{\ell}^{*}} & -\dfrac{\Delta Y_{\ell}^{*}}{S_{\ell}^{*}} & -\dfrac{\Delta Z_{\ell}^{*}}{S_{\ell}^{*}} & 1 \end{bmatrix} \begin{bmatrix}x\\ y\\ z\\ \Delta\tau\end{bmatrix} +\begin{bmatrix}\Delta_1\\ \Delta_2\\ \vdots\\ \Delta_{\ell}\end{bmatrix} \tag{2.1.38}\] 若函数模型除了观测方程外，还有如式 (2.1.23) 的约束条件，且约束条件为非线性函数，那么约束条件也要与观测方程一并进行线性化。线性化方法与观测方程的线性化一样，在近似值 \(\bm{X}^{*}\) 处用泰勒级数展开，舍去二阶和二阶以上的高阶项，有 \[\underset{c\times n}{\bm{C}}\ \underset{n\times 1}{\bm{x}}+\underset{c\times 1}{\bm{\varphi}(\bm{X}^{*})}=0 \tag{2.1.39}\] 其中， \[\bm{C}=\begin{bmatrix}\bm{C}_1\\ \bm{C}_1\\ \vdots\\ \bm{C}_c\end{bmatrix} =\begin{bmatrix} \left(\dfrac{\partial\varphi_1(\bm{X})}{\partial X_1}\right)_{*} & \left(\dfrac{\partial\varphi_1(\bm{X})}{\partial X_2}\right)_{*} & \cdots & \left(\dfrac{\partial\varphi_1(\bm{X})}{\partial X_n}\right)_{*}\\[10pt] \left(\dfrac{\partial\varphi_2(\bm{X})}{\partial X_1}\right)_{*} & \left(\dfrac{\partial\varphi_2(\bm{X})}{\partial X_2}\right)_{*} & \cdots & \left(\dfrac{\partial\varphi_2(\bm{X})}{\partial X_n}\right)_{*}\\[10pt] \vdots & \vdots & \ddots & \vdots\\[4pt] \left(\dfrac{\partial\varphi_c(\bm{X})}{\partial X_1}\right)_{*} & \left(\dfrac{\partial\varphi_c(\bm{X})}{\partial X_2}\right)_{*} & \cdots & \left(\dfrac{\partial\varphi_c(\bm{X})}{\partial X_n}\right)_{*} \end{bmatrix} \tag{2.1.40}\] \[\bm{\varphi}(\bm{X}^{*})=\begin{bmatrix}\varphi_1(\bm{X}^{*})\\ \varphi_2(\bm{X}^{*})\\ \vdots\\ \varphi_c(\bm{X}^{*})\end{bmatrix} \tag{2.1.41}\]

例 2.4将例 2.3 中的观测方程和约束条件线性化。

解：此问题中的未知参数为 \(\bm{X}=\left[\begin{array}{llll}X_{P_1} & Y_{P_1} & X_{P_2} & Y_{P_2}\end{array}\right]^{\mathrm{T}}\)，观测方程为 \(Z_{ij}=\left(\sqrt{\left(X^i-X_{P_j}\right)^2+\left(Y^i-Y_{P_j}\right)^2}\right)+\Delta_{ij}\)

设待估计参数的近似值为 \(\bm{X}=\left[\begin{array}{llll}X_{P_1}^{*} & Y_{P_1}^{*} & X_{P_2}^{*} & Y_{P_2}^{*}\end{array}\right]^{\mathrm{T}}\)，现以基站 \(A_1\) 与目标 \(P_1\) 的边长观测值 \(Z_{A_1P_1}\) 为例进行线性化。用泰勒公式将上式展开并舍去二阶和高阶项： \[Z_{A_1P_1}=Z_{A_1P_1}^{*}+\begin{bmatrix}-\dfrac{X^{A_1}-X_{P_1}^{*}}{Z_{A_1P_1}^{*}} & -\dfrac{Y^{A_1}-Y_{P_1}^{*}}{Z_{A_1P_1}^{*}} & 0 & 0\end{bmatrix} \begin{bmatrix}x_{P_1}\\ y_{P_1}\\ x_{P_2}\\ y_{P_2}\end{bmatrix}+\Delta_{A_1P_1} \tag{2.1.42}\] 其中， \[Z_{A_1P_1}^{*}=\left(\sqrt{\left(X^{A_1}-X_{P_1}^{*}\right)^2+\left(Y^{A_1}-Y_{P_1}^{*}\right)^2}\right)\ ,\quad \bm{x}=\begin{bmatrix}x_{P_1}\\ y_{P_1}\\ x_{P_2}\\ y_{P_2}\end{bmatrix} =\begin{bmatrix}X_{P_1}-X_{P_1}^{*}\\ Y_{P_1}-Y_{P_1}^{*}\\ X_{P_2}-X_{P_2}^{*}\\ Y_{P_2}-Y_{P_2}^{*}\end{bmatrix} \tag{2.1.43}\] 设 \(z_{A_1P_1}=Z_{A_1P_1}-Z_{A_1P_1}^{*}\) 那么，观测方程为 \[z_{A_1P_1}=\begin{bmatrix}-\dfrac{X^{A_1}-X_{P_1}^{*}}{Z_{A_1P_1}^{*}} & -\dfrac{Y^{A_1}-Y_{P_1}^{*}}{Z_{A_1P_1}^{*}} & 0 & 0\end{bmatrix} \begin{bmatrix}x_{P_1}\\ y_{P_1}\\ x_{P_2}\\ y_{P_2}\end{bmatrix}+\Delta_{A_1P_1} \tag{2.1.44}\] 例 2.3 中约束条件为 \[\varphi(\bm{X})=\sqrt{\left(X_{P_1}-X_{P_2}\right)^2+\left(Y_{P_1}-Y_{P_2}\right)^2}-80.50=0\] 线性化后的约束条件为 \[\varphi(\bm{X})=\varphi(\bm{X}^{*})+\left(\frac{\partial\varphi}{\partial\bm{X}}\right)_{*}\bm{x}=0 \tag{2.1.45}\] 其中 \[\varphi(\bm{X}^{*})=\sqrt{\left(X_{P_1}^{*}-X_{P_2}^{*}\right)^2+\left(Y_{P_1}^{*}-Y_{P_2}^{*}\right)^2}-80.50\] \[\left(\frac{\partial\varphi}{\partial\bm{X}}\right)_{*}\bm{x} =\begin{bmatrix}\dfrac{X_{P_1}^{*}-X_{P_2}^{*}}{S_{P_1P_2}^{*}} & \dfrac{Y_{P_1}^{*}-Y_{P_2}^{*}}{S_{P_1P_2}^{*}} & -\dfrac{X_{P_1}^{*}-X_{P_2}^{*}}{S_{P_1P_2}^{*}} & -\dfrac{Y_{P_1}^{*}-Y_{P_2}^{*}}{S_{P_1P_2}^{*}}\end{bmatrix} \begin{bmatrix}x_{P_1}\\ y_{P_1}\\ x_{P_2}\\ y_{P_2}\end{bmatrix}\] \[S_{P_1P_2}^{*}=\sqrt{\left(X_{P_1}^{*}-X_{P_2}^{*}\right)^2+\left(Y_{P_1}^{*}-Y_{P_2}^{*}\right)^2} \tag{2.1.46}\]

2. 线性化带来的模型误差

从以上的线性化过程也可以看出，在舍弃高阶项时也给函数模型也带来了误差，即线性化带来的误差，线性化带来的误差如图 2.4 所示。图中的曲线为一维参数情况下的非线性函数 \(f(X)\)，点 \(X^{*}\) 处为函数值 \(f(X^{*})\)，虚线是函数在 \(X^{*}\) 处的切线。可以看出，线性化后舍去高阶项后的取值在 \(b\) 点处，即 \[f_b(X)=f(X^{*})+\left(\frac{\mathrm{d}f}{\mathrm{d}X}\right)_{*}(X-X^{*}) \tag{2.1.47}\] 而未线性化的函数值 \(f(X)=f(X^{*}+\Delta X)\) 在 \(c\) 处取值。线性化前后的差异为 \[bc=f(X)-\left[\,f(X^{*})+\left(\frac{\mathrm{d}f}{\mathrm{d}X}\right)_{*}(X-X^{*})\,\right] \tag{2.1.48}\] \(bc\) 即为泰勒级数中舍去的高阶项。函数的非线性化程度越高，\(bc\) 越大，即线性化带来的模型误差就越大。此外，\(X^{*}\) 与 \(X\) 的差异越大，\(bc\) 也就越大。所以在实际应用时，一般根据经验或者预测方法取得参数的初始值 \(X^{*}\)，然后在估计时采用迭代方法使 \(X^{*}\) 尽可能地接近 \(X\) 来减小线性化带来的误差。在后面的最小二乘估计中，将介绍如何通过迭代计算来减小线性化带来的模型误差。

最小二乘估计

最小二乘估计方法是由德国数学家高斯（C. F. Gauss，1777—1855）提出的，它的准则是使得残差的加权平方和最小。本节介绍的是基于 2.1 节中介绍的高斯-马尔可夫模型推导的最小二乘估计，它简单直观，易于编程实现，是目前应用最广泛的估计方法之一。

在现实应用中，最小二乘估计有不同的实现方式：仅利用当前时刻观测值进行的最小二乘估计，也被称为“snapshot”最小二乘估计；利用所有观测值，将所有观测值“堆放”在一起，被称为“batch”最小二乘估计；利用当前观测值对参数不断地进行更新，被称为“递推”的最小二乘估计。无论哪一种实现方式，其实质都是基于高斯-马尔可夫模型并且使残差的加权平方和最小的估计。

最小二乘估计

1. 最小二乘估计

上一节得到的高斯-马尔可夫模型为： \[\underset{\ell\times 1}{\bm{Z}}=\underset{\ell\times n}{\bm{H}}\ \underset{n\times 1}{\bm{X}}+\underset{\ell\times 1}{\bm{\Delta}} \tag{2.2.1}\] \[\begin{aligned} E(\bm{Z})&=\bm{H}\bm{X}\\ \mathrm{Var}(\bm{Z})&=\bm{D} \end{aligned} \tag{2.2.2}\] 现假设通过某种方法估计得到参数 \(\bm{X}\) 的估计 \(\hat{\bm{X}}\)，代入观测方程后可得到估计的观测值 \[\underset{\ell\times 1}{\hat{\bm{Z}}}=\underset{\ell\times n}{\bm{H}}\ \underset{n\times 1}{\hat{\bm{X}}} \tag{2.2.3}\] 设估计观测值 \(\hat{\bm{Z}}\) 与观测值 \(\bm{Z}\) 的差异为： \[\bm{v}=\hat{\bm{Z}}-\bm{Z}=\bm{H}\hat{\bm{X}}-\bm{Z} \tag{2.2.4}\] 其中 \[\bm{v}=\begin{bmatrix}v_1\\ v_2\\ \vdots\\ v_{\ell}\end{bmatrix} \tag{2.2.5}\] 上式中的 \(v_i\) 表示观测值 \(Z_i\) 的“残差”，\(\bm{v}\) 也称为观测值的残差向量。由 2.1.1 节的分析可知，式 (2.2.4) 有无穷多组解，如果给它某种最优准则，那么满足这种准则下的解就是“最优”解。最小二乘估计的准则是：参数估计使观测值的残差平方和最小，表达为： \[L(\hat{\bm{X}})=\sum_{i=1}^{i=\ell}v_i^2=\min \tag{2.2.6}\] 其中 \(L(\hat{\bm{X}})\) 为目标函数，也称为估计损失函数。用向量的形式表示为： \[L(\hat{\bm{X}})=\bm{v}^{\mathrm{T}}\bm{v}=\min \tag{2.2.7}\]

在方程 (2.2.4) 的无穷多组解中，能够满足上式的解即为最小二乘解 \(\hat{\bm{X}}_{LS}\) \[\hat{\bm{X}}_{LS}=\arg\,\min_{\hat{\bm{X}}}L(\hat{\bm{X}}) \tag{2.2.8}\] 上式中的 \(\arg\,\min\limits_{\hat{\bm{X}}}L(\hat{\bm{X}})\) 表示使 \(L(\hat{\bm{X}})\) 最小值时的变量的取值。准则 (2.2.7) 视所有的观测值 \(Z_1\)，\(Z_2\)，…，\(Z_{\ell}\) 对参数的估计的影响是相同的，但有时观测值的精度是不同的，所以在估计时我们希望精度好的观测值能比精度差的观测值对参数估计产生的影响大，即方差小的观测值对参数估计的影响大，反之亦然。如果观测值随机独立，\(Z_i\) 的方差为 \(\sigma_i^2\)，那么给观测值 \(Z_i\) 赋予的影响因子为 \[w_i=\frac{\sigma_0^2}{\sigma_i^2} \tag{2.2.9}\] 上式中 \(\sigma_0^2\) 为任意正实数，\(w_i\) 与方差 \(\sigma_i^2\) 成反比，称为观测值 \(Z_i\) 的权。从后面的证明也可以看到 \(\sigma_0^2\) 的取值并不影响最小二乘参数估计值。考虑权因子，最小二乘准则为 \[\sum_{i=1}^{i=\ell}v_i^2w_i=\min \tag{2.2.10}\] 上式也称为加权最小二乘准则。若将观测值的权表述为矩阵 \[\bm{W}=\begin{bmatrix}w_1 & & &\\ & w_2 & &\\ & & \ddots &\\ & & & w_{\ell}\end{bmatrix} =\begin{bmatrix}\dfrac{\sigma_0^2}{\sigma_1^2} & & &\\[6pt] & \dfrac{\sigma_0^2}{\sigma_2^2} & &\\[6pt] & & \ddots &\\[6pt] & & & \dfrac{\sigma_0^2}{\sigma_{\ell}^2}\end{bmatrix} \tag{2.2.11}\] 取 \(\sigma_0^2=\sigma_i^2\)，那么观测值 \(Z_i\) 的权为 \(1\)，所以称 \(\sigma_0^2\) 为“单位权方差”。权矩阵 \(\bm{W}\) 给出了观测值间精度的比例关系，当我们无法确定观测值的绝对精度的时候，可以根据经验给出观测值精度的比例关系。例如，知道观测值 \(Z_1\) 的观测中误差为 \(Z_2\) 的两倍，那么在 \(\bm{W}\) 矩阵中，\(w_1=\dfrac{1}{4}w_2\)。当观测值 \(Z_i\) 的中误差非常大，甚至无穷大时，那么其对应的观测值的权 \(w_i\) 可设为零，这意味着 \(Z_i\) 对参数估计的影响为“零”，这与剔除观测值 \(Z_i\) 进行参数估计的效果是一样的。

在更一般的情况下，观测值随机相关，设观测值的方差矩阵为 \(\bm{D}\)，那么观测值向量的权矩阵为 \[\bm{W}=\sigma_0^2\,\bm{D}^{-1} \tag{2.2.12}\] 这时的最小二乘准则为 \[L(\bm{v})=\bm{v}^{\mathrm{T}}\bm{W}\bm{v}=\min \tag{2.2.13}\]

接下来推导满足式 (2.2.13) 的参数估计。设满足式 (2.2.13) 的参数估计为 \(\hat{\bm{X}}_{LS}\)，它满足目标函数： \[L(\hat{\bm{X}}_{LS})=\bm{v}^{\mathrm{T}}\bm{W}\bm{v}=\min \tag{2.2.14}\] 与式 (2.2.7) 比较，式 (2.2.14) 给出了更一般情况下的最小二乘准则。

将误差方程 (2.2.4) 代入目标函数式 (2.2.14) \[\begin{aligned} L(\hat{\bm{X}}_{LS})&=(\bm{H}\hat{\bm{X}}_{LS}-\bm{Z})^{\mathrm{T}}\bm{W}(\bm{H}\hat{\bm{X}}_{LS}-\bm{Z})\\ &=\hat{\bm{X}}_{LS}^{\mathrm{T}}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{X}}_{LS} -\hat{\bm{X}}_{LS}^{\mathrm{T}}\bm{H}^{\mathrm{T}}\bm{W}\bm{Z} -\bm{Z}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{X}}_{LS}+\bm{Z}^{\mathrm{T}}\bm{W}\bm{Z} \end{aligned} \tag{2.2.15}\] \(L(\hat{\bm{X}}_{LS})\) 是 \(\hat{\bm{X}}_{LS}\) 的函数，为了使其最小，由函数极值方法得到 \[\frac{\partial L(\hat{\bm{X}}_{LS})}{\partial\hat{\bm{X}}_{LS}}=2\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{X}}_{LS}-2\bm{H}^{\mathrm{T}}\bm{W}\bm{Z}=0 \tag{2.2.16}\] 即 \[\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{X}}_{LS}-\bm{H}^{\mathrm{T}}\bm{W}\bm{Z}=0 \tag{2.2.17}\] 上式也称为“法方程”。由于 \(\bm{H}\) 为列满秩矩阵，且 \(\bm{W}\) 为满秩方阵，所以 \(\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\) 为满秩矩阵 \[\mathrm{rank}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})=n \tag{2.2.18}\]

根据式 (2.2.16) 解出参数估计 \(\hat{\bm{X}}_{LS}\) \[\hat{\bm{X}}_{LS}=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{Z} \tag{2.2.19}\] 现令 \[\bm{N}=\bm{H}^{\mathrm{T}}\bm{W}\bm{H} \tag{2.2.20}\] 和 \[\bm{Q}_{\hat{\bm{X}}_{LS}}=\bm{N}^{-1} \tag{2.2.21}\] 那么 \[\hat{\bm{X}}_{LS}=\bm{Q}_{\hat{\bm{X}}_{LS}}\bm{H}^{\mathrm{T}}\bm{W}\bm{Z} \tag{2.2.22}\] 也称 \(\bm{Q}_{\hat{\bm{X}}_{LS}}\) 为参数解 \(\hat{\bm{X}}_{LS}\) 的协因数矩阵。\(\hat{\bm{X}}_{LS}\) 是满足残差平方和最小的估计，也称为“最小二乘估计”。

将 \(\bm{W}=\sigma_0^2\,\bm{D}^{-1}\) 代入式 (2.2.19) 得到 \[\hat{\bm{X}}_{LS}=\left(\bm{H}^{\mathrm{T}}(\sigma_0^2\,\bm{D}^{-1})\bm{H}\right)^{-1}\bm{H}^{\mathrm{T}}(\sigma_0^2\,\bm{D}^{-1})\bm{Z} \tag{2.2.23}\] 将上式中的 \(\sigma_0^2\) 约去，得到 \[\hat{\bm{X}}_{LS}=(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{Z} \tag{2.2.24}\] 从上面的推导看出，无论用式 (2.2.19)，还是用式 (2.2.24)，得到的参数估计都是等价的。上面的推导过程也说明了确定权矩阵时 \(\sigma_0^2\) 的数值并不影响 \(\hat{\bm{X}}_{LS}\)，\(\hat{\bm{X}}_{LS}\) 的估计值由设计矩阵 \(\bm{H}\)、观测值 \(\bm{Z}\) 和观测值的精度比例关系决定。当方差矩阵 \(\bm{D}\) 未知时，可用根据经验确定观测值精度的比例关系，用式 (2.2.19) 计算 \(\hat{\bm{X}}_{LS}\)。在已知方差矩阵 \(\bm{D}\) 时，可用式 (2.2.23) 进行估计，这时的单位权中误差 \(\sigma_0^2\) 默认为“1”。

将式 (2.2.19) 代入误差方程式 (2.2.4) 中可计算残差 \[\bm{v}=\left(\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}-\bm{I}\right)\bm{Z} \tag{2.2.25}\] 令 \(\bm{P}_H=\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\)，并将观测方程式 (2.1.1) 代入上式可以得到 \[\bm{v}=-\left(\bm{I}-\bm{P}_H\right)\bm{\Delta} \tag{2.2.26}\] 上式表明矩阵 \(-(\bm{I}-\bm{P}_H)\) 将观测误差映射于残差向量，当观测误差有异常时，如有粗差（错误观测），粗差将会在残差向量上有所体现，所以我们可以通过观察残差向量或者对残差向量进行假设检验来发现观测值中是否有粗差。

可以证明，上式中的矩阵 \((\bm{I}-\bm{P}_H)\) 为幂等矩阵，幂等矩阵的秩等于其矩阵的迹的绝对值，因此 \[\begin{aligned} \mathrm{rank}(\bm{I}-\bm{P}_H)&=\mathrm{abs}\left(\mathrm{tr}\left(\bm{I}-\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\right)\right)\\ &=\mathrm{abs}\left(\mathrm{tr}(\bm{I})-\mathrm{tr}\left((\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\right)\right)\\ &=\mathrm{abs}(\ell-n)\\ &=\ell-n \end{aligned} \tag{2.2.27}\] 上式中的 \(\ell-n\) 也是此估计问题的多余观测数，即自由度。

2. 最小二乘估计的统计特性

最小二乘估计的期望 \[\begin{aligned} E(\hat{\bm{X}}_{LS})&=E\left[\,(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{Z}\,\right]\\ &=E\left[\,(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}(\bm{H}\bm{X}+\bm{\Delta})\,\right]\\ &=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}E(\bm{X})+(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}E(\bm{\Delta})\\ &=\bm{X} \end{aligned} \tag{2.2.28}\] 上式表明最小二乘估计 \(\hat{\bm{X}}_{LS}\) 的期望为 \(\bm{X}\)，即最小二乘估计为无偏估计。

残差的期望为 \[\begin{aligned} E(\bm{v})&=E(\hat{\bm{Z}}-\bm{Z})\\ &=\bm{H}E(\hat{\bm{X}}_{LS})-E(\bm{Z})\\ &=\bm{H}\bm{X}-\bm{H}\bm{X}\\ &=0 \end{aligned} \tag{2.2.29}\] 上式表明残差期望为零，这与观测误差 \(\bm{\Delta}\) 的期望一致。

由于最小二乘估计为无偏估计，所以最小二乘估计的方差也为均方差： \[\mathrm{Var}(\hat{\bm{X}}_{LS})=E\left[\,(\hat{\bm{X}}_{LS}-\bm{X})(\hat{\bm{X}}_{LS}-\bm{X})^{\mathrm{T}}\,\right] \tag{2.2.30}\] 由于 \[\begin{aligned} \hat{\bm{X}}_{LS}-\bm{X}&=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}(\bm{H}\bm{X}+\bm{\Delta})-\bm{X}\\ &=\bm{X}+(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{\Delta}-\bm{X}\\ &=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{\Delta} \end{aligned} \tag{2.2.31}\] 将式 (2.2.31) 代入式 (2.2.30)，得到 \[\mathrm{Var}(\hat{\bm{X}}_{LS})=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}E(\bm{\Delta}\bm{\Delta}^{\mathrm{T}})\bm{W}\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1} \tag{2.2.32}\] 上式中的 \(E(\bm{\Delta}\bm{\Delta}^{\mathrm{T}})\) 即为随机误差的方差 \(\mathrm{Var}(\bm{\Delta})\) \[\mathrm{Var}(\bm{\Delta})=\bm{D} \tag{2.2.33}\] 由于 \[\bm{W}=\sigma_0^2\,\bm{D}^{-1} \tag{2.2.34}\] 得到 \[\begin{aligned} \mathrm{Var}(\hat{\bm{X}}_{LS})&=(\bm{H}^{\mathrm{T}}\sigma_0^2\bm{D}^{-1}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\sigma_0^2\bm{D}^{-1}\bm{D}\sigma_0^2\bm{D}^{-1}\bm{H}(\bm{H}^{\mathrm{T}}\sigma_0^2\bm{D}^{-1}\bm{H})^{-1}\\ &=(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} \end{aligned} \tag{2.2.35}\] 考虑式 (2.2.34)，\(\mathrm{Var}(\hat{\bm{X}}_{LS})\) 也为 \[\mathrm{Var}(\hat{\bm{X}}_{LS})=\sigma_0^2(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1} \tag{2.2.36}\]

由式 (2.2.36) 并根据误差传播定律，可以求得残差的方差为 \[\mathrm{Var}(\bm{v})=\bm{D}-\bm{H}(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1}\bm{H}^{\mathrm{T}} \tag{2.2.37}\] 残差 \(\bm{V}\) 与 \(\hat{\bm{X}}_{LS}\) 的协方差为 \[\begin{aligned} \mathrm{Cov}(\bm{v},\ \hat{\bm{X}}_{LS})&=(\bm{P}_H-\bm{I})\bm{D}\left[\,(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\,\right]^{\mathrm{T}}\\ &=(\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}-\bm{I})\bm{D}\left[\,(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\,\right]^{\mathrm{T}}\\ &=-\bm{D}\bm{W}\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}+\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{D}\bm{W}\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\\ &=-\sigma^2\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}+\sigma^2\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\\ &=0 \end{aligned} \tag{2.2.38}\] 上式表明残差 \(\bm{V}\) 与 \(\hat{\bm{X}}_{LS}\) 不相关。

3. 验后估计单位权方差 \(\hat{\sigma}_0^2\) 和应用

在确定权矩阵时，由于不知道观测值的绝对精度，我们可以设定任意数值的单位权中误差 \(\sigma_0^2\)，这并不影响最小二乘估计。在得到最小二乘估计后，我们可以利用权矩阵和残差对单位权中误差 \(\sigma_0^2\) 进行估计，从而了解观测值的绝对精度。通过观测值估计得到的单位权方差称为验后单位权中误差 \(\hat{\sigma}_0^2\)。下面推导如何利用观测值估计验后单位权中误差 \(\hat{\sigma}_0^2\)。

最小二乘估计的准则的目标函数 \(L(\hat{\bm{X}})\) 的期望为 \[E\left(\,L(\hat{\bm{X}})\,\right)=E(\bm{v}^{\mathrm{T}}\bm{W}\bm{v}) \tag{2.2.39}\] 根据二次型定理（附录 C-4）有 \[\begin{aligned} E\left(\,L(\hat{\bm{X}})\,\right)&=E(\bm{v}^{\mathrm{T}}\bm{W}\bm{v})\\ &=\mathrm{tr}\left(\bm{W}\mathrm{Var}(\bm{v})\right)+E^{\mathrm{T}}(\bm{v})\bm{W}E(\bm{v}) \end{aligned} \tag{2.2.40}\] 将式 (2.2.37) 和式 (2.2.29) 代入上式，得 \[\begin{aligned} E(\bm{v}^{\mathrm{T}}\bm{W}\bm{v})&=\sigma_0^2\,\mathrm{tr}\left(\bm{I}-\bm{W}\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\right)\\ &=\sigma_0^2\left[\,\mathrm{tr}(\bm{I})-\mathrm{tr}\left((\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\right)\,\right]\\ &=\sigma_0^2(\ell-n) \end{aligned} \tag{2.2.41}\] 可以看到 \(\dfrac{\bm{v}^{\mathrm{T}}\bm{W}\bm{v}}{\ell-n}\) 是单位权方差 \(\sigma_0^2\) 的无偏估计，因此有 \[\begin{aligned} \hat{\sigma}_0^2&=\frac{\bm{v}^{\mathrm{T}}\bm{W}\bm{v}}{\ell-n}\\ &=\frac{\sigma_0^2\,\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}}{\ell-n} \end{aligned} \tag{2.2.42}\] \(\hat{\sigma}_0^2\) 由残差计算得到，包含有观测值的信息，当观测值中有粗差，或者数学模型与实际不符时，在 \(\hat{\sigma}_0^2\) 上都可以得到体现。这里不加证明的给出，当数学模型符合式 (2.2.1) 和式 (2.1.2) 的高斯-马尔可夫模型时，\(\dfrac{\hat{\sigma}_0^2(\ell-n)}{\sigma_0^2}\) 服从自由度为 \((\ell-n)\)，中心化参数 \(\lambda=0\) 的卡方分布 \(\chi^2_{(\ell-n,\ \lambda=0)}\) \[t=\frac{\hat{\sigma}_0^2(\ell-n)}{\sigma_0^2}=\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}\sim\chi^2_{(\ell-n,\ \lambda=0)} \tag{2.2.43}\] 所以当得到 \(\hat{\sigma}_0^2\) 后，可对 \(\hat{\sigma}_0^2\) 构成的检验量进行假设检验。原假设为 \[H_0:\ \text{模型符合式 (2.2.1) 和式 (2.2.2) 的假设},t\sim\chi^2_{(\ell-n,\ \lambda=0)} \tag{2.2.44}\] 备选假设为 \[H_{\alpha}:\ \text{模型不符合 (2.2.1) 和式 (2.2.2) 的假设},t\sim\chi^2_{(\ell-n,\ \lambda\neq 0)}\] 给定的显著性水平 \(\alpha\)（误警概率）和自由度 \((\ell-n)\)，可得到分位置 \(T_{1-\alpha}\)（检测限值），如果 \[t<T_{1-\alpha} \tag{2.2.45}\] 则接受原假设，这时可以计算得到验后的估计参数精度 \[\begin{aligned} \mathrm{Var}(\hat{\bm{X}}_{LS})&=\hat{\sigma}_0^2(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\\ &=\frac{\hat{\sigma}_0^2}{\sigma_0^2}(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} \end{aligned} \tag{2.2.46}\] 如果 \(t\geq T_{1-\alpha}\)，接受备选假设，认为观测值向量中有粗差，或者给定的验方差矩阵 \(\bm{D}\) 与实际不符。对模型的假设检验并不能辨别到底是哪种原因造成拒绝原假设，一般在方差准确的情况下检验粗差，或者认为观测值无粗差的情况下来检验方差是否准确。

在对粗差观测进行检验时，如果接受了备选假设，通常将残差最大的观测值剔除，或者通过巴尔达（Baarda）检测来剔除观测值，然后重新进行最小二乘估计，并再次进行假设检验。当具备足够多的多余观测并且只有一个粗差观测值时，上述方法一般可以准确的识别粗差观测。但当存在两个或两个以上的粗差观测时，会错误识别或者无法识别粗差观测。在安全性能要求较高的导航应用中，如 GNSS 为民航提供导航时，如果无法识别和剔除故障卫星（粗差观测），应对用户提出告警，此时的 GNSS 不能提供导航。

4. 最小二乘估计的正交特性

最小二乘估计的正交特性在理论证明中有广泛的应用，了解它的正交特性对最小二乘估计理论也有更好的理解。

由式 (2.2.26) 得到 \[\bm{v}=(\bm{P}_H-\bm{I})\bm{Z} \tag{2.2.47}\] 根据式 (2.2.3) 可得到改正后的观测值为 \[\begin{aligned} \hat{\bm{Z}}&=\bm{H}(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{Z}\\ &=\bm{P}_H\bm{Z} \end{aligned} \tag{2.2.48}\] 由于 \(\bm{P}_H\) 为幂等矩阵，所以有 \[(\bm{P}_H-\bm{I})\bm{P}_H=0 \tag{2.2.49}\] 上式表明矩阵 \((\bm{P}_H-\bm{I})\) 的行与矩阵 \(\bm{P}_H\) 的列的内积为零，这意味着 \((\bm{P}_H-\bm{I})\) 与 \(\bm{P}_H\) 正交。又由于 \[\mathrm{rank}(\bm{I}-\bm{P}_H)=\ell-n \tag{2.2.50}\] 和 \[\mathrm{rank}(\bm{P}_H)=n \tag{2.2.51}\] 这表明由矩阵 \(\bm{P}_H\) 构成的 \(n\) 维向量空间 \(V_H\) 与由 \((\bm{P}_H-\bm{I})\) 构成的 \(\ell-n\) 维向量空间 \(V_H^{\perp}\) 正交，\(V_H^{\perp}\) 为 \(V_H\) 的正交补。\(V_H\) 和 \(V_H^{\perp}\) 构成了 \(\ell\) 维的向量空间 \(V^{\ell}\)，即 \[V^{\ell}=V_H^{\perp}\oplus V_H \tag{2.2.52}\] 由于 \(\hat{\bm{Z}}\) 为 \(\bm{P}_H\) 的线性组合，\(\bm{v}\) 为 \((\bm{P}_H-\bm{I})\) 的线性组合，所以 \(\hat{\bm{Z}}\) 与 \(\bm{v}\) 必然正交，有 \[(\hat{\bm{Z}},\ \bm{v})_W=\hat{\bm{Z}}^{\mathrm{T}}\bm{W}\bm{v}=0 \tag{2.2.53}\] 其中 \(\bm{W}\) 为观测值的权，\((\hat{\bm{Z}},\ \bm{v})_W\) 为 \(\hat{\bm{Z}}\) 与 \(\bm{v}\) 的广义内积。

除了以上的正交关系，还有 \[\bm{H}^{\mathrm{T}}\bm{W}(\bm{P}_H-\bm{I})=0 \tag{2.2.54}\] 上式表明残差向量 \(\bm{v}\) 还与设计矩阵 \(\bm{H}\) 的各列正交。观测值 \(\bm{Z}\)、改正后的观测值 \(\hat{\bm{Z}}\) 和残差向量 \(\bm{v}\) 之间的几何关系可以用图 2.5 表示：矩阵 \(\bm{P}_H\) 将观测值 \(Z\) 投影到空间 \(V_H\) 上得到 \(\hat{\bm{Z}}\)，矩阵 \((\bm{P}_H-\bm{I})\) 将观测值投影到空间 \(V_H^{\perp}\) 得到 \(\bm{v}\)，所以 \(\hat{\bm{Z}}\) 与 \(\bm{v}\) 正交。

附有约束条件的最小二乘估计

1. 参数估计

这里以线性化后的函数模型来推导附有约束条件的最小二乘估计。

观测方程为 \[\bm{z}=\bm{H}\bm{x}+\bm{\Delta} \tag{2.2.55}\] 约束条件为 \[\bm{C}\bm{x}+\bm{\varphi}(\bm{X}^{*})=0 \tag{2.2.56}\] 随机模型为 \[\begin{aligned} E(\bm{\Delta})&=0\\ \mathrm{Var}(\bm{\Delta})&=\bm{D} \end{aligned} \tag{2.2.57}\] 其中 \(\bm{H}\) 为列满秩矩阵；矩阵 \(\bm{C}\) 为行满秩矩阵，其秩为参数的约束条件的个数，表明若有多个约束条件，各个条件之间函数独立，而且 \(c<n\)。有了参数约束条件后，必要观测值数为 \(n-c\)，多余观测数为 \(l-(n-c)\)。

假设 \(\bm{x}\) 的最小二乘估计 \(\hat{\bm{x}}_{LS}\)，那么误差修正后的观测值 \[\hat{\bm{z}}=\bm{H}\hat{\bm{x}}_{LS} \tag{2.2.58}\] 观测值的残差为 \[\bm{v}=\bm{H}\hat{\bm{x}}_{LS}-\bm{z} \tag{2.2.59}\] 同时 \(\hat{\bm{x}}_{LS}\) 还要满足约束条件 \[\bm{C}\hat{\bm{x}}_{LS}+\bm{\varphi}(\bm{X}^{*})=0 \tag{2.2.60}\] 联立式 (2.2.60) 和式 (2.2.59) 得到 \[\begin{bmatrix}\bm{I} & -\bm{H}\\ 0 & \bm{C}\end{bmatrix} \begin{bmatrix}\bm{v}\\ \hat{\bm{x}}_{LS}\end{bmatrix} =\begin{bmatrix}-\bm{z}\\ -\bm{\varphi}(\bm{X}^{*})\end{bmatrix} \tag{2.2.61}\] 在式 (2.2.61) 中，方程的个数为 \(\ell+c\)，被估计量的个数为 \(\ell+n\)，所以方程的个数小于被估价量的个数，且系数矩阵的秩等于其增广矩阵的秩，即 \[\mathrm{rank}\begin{bmatrix}\bm{I} & -\bm{H}\\ 0 & \bm{C}\end{bmatrix} =\mathrm{rank}\left[\begin{array}{cc|c}\bm{I} & -\bm{H} & -\bm{z}\\ 0 & \bm{C} & -\bm{\varphi}(\bm{X}^{*})\end{array}\right] \tag{2.2.62}\] 这表明式 (2.2.62) 是有无穷组解的相容方程。现在要在这无穷多组解中找出能够满足最小二乘准则 \(\bm{v}^{\mathrm{T}}\bm{W}\bm{v}=\min\)，并且满足式 (2.2.60) 的一组解。按照拉格朗日乘数法，目标函数为 \[L(\hat{\bm{x}})=\bm{v}^{\mathrm{T}}\bm{W}\bm{v}+2\bm{K}^{\mathrm{T}}\left(\bm{C}\bm{x}+\bm{\varphi}(\bm{X}^{*})\right) \tag{2.2.63}\] 式中 \(\bm{K}\) 是对应于约束条件方程的联系数向量。现将式 (2.2.59) 代入上式得到 \[L(\hat{\bm{x}})=\hat{\bm{x}}_{LS}^{\mathrm{T}}\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{x}}_{LS} -\hat{\bm{x}}_{LS}^{\mathrm{T}}\bm{H}^{\mathrm{T}}\bm{W}\bm{z}-\bm{z}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{x}}_{LS} +\bm{z}^{\mathrm{T}}\bm{W}\bm{z}+2\bm{K}^{\mathrm{T}}\left(\bm{C}\hat{\bm{x}}_{LS}+\bm{\varphi}(\bm{X}^{*})\right) \tag{2.2.64}\] 为了求 \(L(\hat{\bm{x}})\) 的极小值，将其对 \(\bm{x}\) 求导并令其为零，整理后则有 \[\bm{H}^{\mathrm{T}}\bm{W}\bm{H}\hat{\bm{x}}_{LS}-\bm{H}^{\mathrm{T}}\bm{W}\bm{z}+\bm{C}^{\mathrm{T}}\bm{K}=0 \tag{2.2.65}\] 与约束条件 (2.2.60) 联立有 \[\begin{bmatrix}\bm{H}^{\mathrm{T}}\bm{W}\bm{H} & \bm{C}^{\mathrm{T}}\\ \bm{C} & \bm{0}\end{bmatrix} \begin{bmatrix}\hat{\bm{x}}_{LS}\\ \bm{K}\end{bmatrix} =\begin{bmatrix}\bm{H}^{\mathrm{T}}\bm{W}\bm{z}\\ -\bm{\varphi}(\bm{X}^{*})\end{bmatrix} \tag{2.2.66}\] 上式是在有约束条件下得到的法方程，方程的个数为 \(n+c\)，未知数的个数也为 \(n+c\)，\(\left[\begin{array}{ll}\hat{\bm{x}}_{LS} & \bm{K}\end{array}\right]^{\mathrm{T}}\) 前的系数矩阵满秩，所以可以直接对法方程系数矩阵求逆。

令 \[\bm{N}=\bm{H}^{\mathrm{T}}\bm{W}\bm{H} \tag{2.2.67}\] 解得 \[\begin{bmatrix}\hat{\bm{x}}_{LS}\\ \bm{K}\end{bmatrix} =\begin{bmatrix}\bm{N} & \bm{C}^{\mathrm{T}}\\ \bm{C} & \bm{0}\end{bmatrix}^{-1} \begin{bmatrix}\bm{H}^{\mathrm{T}}\bm{W}\bm{z}\\ -\bm{\varphi}(\bm{X}^{*})\end{bmatrix} \tag{2.2.68}\] 以上是 \(\left[\begin{array}{ll}\hat{\bm{x}}_{LS} & \bm{K}\end{array}\right]^{\mathrm{T}}\) 的整体求解。也可以分步求解：先解出 \(\bm{K}\)，再求解 \(\hat{\bm{x}}_{LS}\)。求解步骤如下：

用矩阵 \(\bm{C}\cdot(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\) 左乘式 (2.2.65) 后减去式 (2.2.60) 消去 \(\hat{\bm{x}}_{LS}\)，可先解出联系数矩阵 \(\bm{K}\) \[\bm{K}=(\bm{C}\bm{N}^{-1}\bm{C}^{\mathrm{T}})^{-1}\left(\bm{C}\bm{N}^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{z}+\bm{\varphi}(\bm{X}^{*})\right) \tag{2.2.69}\] 令 \[\bm{N}_{CC}=\bm{C}\bm{N}^{-1}\bm{C}^{\mathrm{T}} \tag{2.2.70}\] 并将式 (2.2.69) 代入式 (2.2.65)，可解出参数 \(\hat{\bm{x}}_{LS}\) \[\hat{\bm{x}}_{LS}=\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right)\bm{H}^{\mathrm{T}}\bm{W}\bm{z} -\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{\varphi}(\bm{X}^{*}) \tag{2.2.71}\] 这样分步解算的 \(\hat{\bm{x}}_{LS}\) 的估计与式 (2.2.68) 的解是等价的，但求解中矩阵求逆的阶数减少了，所以可以提高计算机的解算效率和数值的准确性。

最后，参数估计为 \[\hat{\bm{X}}_{LS}=\bm{X}^{*}+\hat{\bm{x}}_{LS} \tag{2.2.72}\] 将 \(\hat{\bm{x}}_{LS}\) 代入式 (2.2.59) 可以得到残差 \(\bm{v}\)。

2. 估计的统计特性

将式 (2.2.55) 代入式 (2.2.71)，并取期望 \[\begin{aligned} E(\hat{\bm{x}}_{LS})&=\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right) \bm{H}^{\mathrm{T}}\bm{W}E(\bm{H}\bm{x}+\bm{\Delta})-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{\varphi}(\bm{X}^{*})\\ &=\bm{x}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\left(\bm{C}\bm{x}+\bm{\varphi}(\bm{X}^{*})\right) \end{aligned} \tag{2.2.73}\] 考虑式 (2.2.56)，上式的第二项为零，则 \[E(\hat{\bm{x}}_{LS})=\bm{x} \tag{2.2.74}\] 上式表明，附有约束条件的最小二乘参数估计为无偏估计。

\(\hat{\bm{x}}_{LS}\) 的方差可以根据方差传播定律求得 \[\begin{aligned} \mathrm{Var}(\hat{\bm{x}}_{LS})&=\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right) \bm{H}^{\mathrm{T}}\bm{W}\bm{D}\bm{W}\bm{H}\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right)^{\mathrm{T}}\\ &=\sigma_0^2\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right) \end{aligned} \tag{2.2.75}\] 设 \[\bm{Q}_{\hat{\bm{x}}_{LS}}=\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1} \tag{2.2.76}\] \(\bm{Q}_{\hat{\bm{x}}_{LS}}\) 也是附有约束条件的最小二乘参数估计的协因数矩阵。

与无约束的最小二乘估计一样，可证明 \(\dfrac{\bm{v}^{\mathrm{T}}\bm{W}\bm{v}}{\ell-n+c}\) 是单位权方差的无偏估计 \[\hat{\sigma}_0^2=\frac{\bm{v}^{\mathrm{T}}\bm{W}\bm{v}}{\ell-n+c} \tag{2.2.77}\]

最小二乘估计的迭代计算

在 2.1 节中介绍了模型的线性化，在线性化后，得到观测方程 \[\bm{z}=\bm{Z}-\bm{F}(\bm{X}^{*})=\bm{H}\bm{x}+\bm{\Delta} \tag{2.2.78}\] \[\hat{\bm{x}}_{LS}=(\bm{H}^{\mathrm{T}}\bm{W}\bm{H})^{-1}\bm{H}^{\mathrm{T}}\bm{W}\bm{z} \tag{2.2.79}\] 参数估计为 \[\hat{\bm{X}}_{LS}=\bm{X}^{*}+\hat{\bm{x}}_{LS} \tag{2.2.80}\] 如 2.1 节中介绍的，由于线性化时忽略了模型中的高阶项，给模型带来了误差，模型误差的大小不仅与函数的非线性化程度有关，还与 \(\bm{X}\) 和近似值 \(\bm{X}^{*}\) 的差异有关。一般在线性化时，根据经验或者预测给出 \(\bm{X}^{*}\)，近似值 \(\bm{X}^{*}\) 与 \(\bm{X}\) 的差异越大，模型误差就越大。为了减小近似值 \(\bm{X}^{*}\) 选取带来的模型误差，这里介绍通过迭代计算来减小模型误差的方法。

如图 2.6 所示，设第一次估计前线性化时近似值取值为 \(\bm{X}^{*(1)}\)，得到线性化模型 \[\bm{z}^{(1)}=\bm{Z}-\bm{F}\left(\bm{X}^{*(1)}\right)=\bm{H}^{(1)}\bm{x}+\bm{\Delta} \tag{2.2.81}\]

在这个模型上进行第一次估计，设估计值为 \(\hat{\bm{x}}^{(1)}\)，那么第一次迭代估计的参数为 \[\hat{\bm{X}}^{(1)}=\bm{X}^{*(1)}+\hat{\bm{x}}^{(1)} \tag{2.2.82}\] 为了获得与参数 \(\bm{X}\) 数值更为接近的近似值，可以根据第一次估计的结果重新进行线性化，第二次进行线性化的近似值为： \[\bm{X}^{*(2)}=\hat{\bm{X}}^{(1)}=\bm{X}^{*(1)}+\hat{\bm{x}}^{(1)} \tag{2.2.83}\]

第二次线性化函数模型 \[\bm{z}^{(2)}=\bm{Z}-\bm{F}\left(\bm{X}^{*(2)}\right)=\bm{H}^{(2)}\bm{x}+\bm{\Delta} \tag{2.2.84}\] 在此模型上再进行估计。重复以上的步骤，在第 \(k\) 次线性化后有估计值： \[\hat{\bm{X}}^{(k)}=\bm{X}^{*(k)}+\hat{\bm{x}}^{(k)} \tag{2.2.85}\] 若得到的 \(\hat{\bm{x}}^{(k)}\) 足够小，说明 \(\hat{\bm{X}}^{(k)}\) 已无明显的变化，即可终止迭代。在实际计算时可以将 \(\hat{\bm{x}}^{(k)}\) 中绝对值最大的值（无穷范数）或者将 \(\hat{\bm{x}}^{(k)}\) 的长度（2 范数）与预先设定好的限值进行比较，来判断是否停止迭代计算。\(\hat{\bm{x}}^{(k)}\) 的限值可以根据观测值的精度来确定。最后的参数估计为 \[\hat{\bm{X}}_{LS}=\hat{\bm{X}}^{(k)}=\bm{X}^{(1)*}+\hat{\bm{x}}^{(1)}\cdots+\hat{\bm{x}}^{(k)} \tag{2.2.86}\] 在得到参数估计 \(\hat{\bm{X}}^{(k)}\) 后进行精度评定。\(\hat{\bm{X}}^{(k)}\) 的验前方差为 \[\mathrm{Var}\left(\hat{\bm{X}}_{LS}\right)=\sigma_0^2\left(\bm{H}^{(k)\mathrm{T}}\bm{W}\bm{H}^{(k)}\right)^{-1} \tag{2.2.87}\] 残差为 \[\bm{v}=\bm{H}^{(k)}\bm{x}^{(k)}-\bm{z}^{(k)} \tag{2.2.88}\] 验后单位权方差为 \[\begin{aligned} \hat{\sigma}_0^2&=\bm{v}^{\mathrm{T}}\bm{W}\bm{v}/(\ell-n)\\ &=\sigma_0^2\,\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}/(\ell-n) \end{aligned} \tag{2.2.89}\] \(\hat{\bm{X}}^{(k)}\) 的验后方差为 \[\mathrm{Var}\left(\hat{\bm{X}}_{LS}\right)=\hat{\sigma}_0^2\left(\bm{H}^{(k)\mathrm{T}}\bm{W}\bm{H}^{(k)}\right)^{-1} \tag{2.2.90}\]

这里需要注意的是，在每一次迭代计算的时候都需要更新矩阵 \(\bm{H}\) 和向量 \(\bm{z}\)。如第 \(i\) 次迭代计算时，以 \(\hat{\bm{X}}^{(i-1)}\) 作为近似值进行线性化，得到矩阵 \(\bm{H}^{(i)}\) 和向量 \(\bm{z}^{(i)}\)，然后估计得到 \(\hat{\bm{x}}^{(i)}\)。迭代结束后，最后计算残差、验后单位权方差和精度评定。

如果模型中除了观测方程外，还有约束条件，那么约束条件与观测方程一并线性化，其迭代过程和估计方法与上面相同。这里需要说明的是，迭代计算只是在一定程度上减小了线性化带来的模型误差，但它无法补偿线性化时忽略高阶项带来的模型误差，尤其当函数模型的非线性化程度非常高时，就应该考虑模型的二阶项并用非线性最小二乘来估计。

算例分析

例 2.5用最小二乘估计法估计例 2.1 中的参数，并求估计参数的验后估计精度

解：此问题的观测方程为 \[\underset{6\times 1}{\bm{Z}}=\underset{6\times 2}{\bm{H}}\ \underset{2\times 1}{\bm{X}}+\underset{6\times 1}{\bm{\Delta}}\] 未知参数为 \[\bm{X}=\left[\begin{array}{ll}\alpha & \beta\end{array}\right]^{\mathrm{T}}\] 观测值数 \(\ell\) 为 6，必要观测值数为 2，所以此估计问题的自由度为 \(r=\ell-n=4\)。由式 (2.1.14) 得到 \[\bm{Z}=\begin{bmatrix}4.16\\ 4.52\\ \vdots\\ 9.26\end{bmatrix}\qquad \bm{H}=\begin{bmatrix}1 & 1\\ 1 & 2\\ \vdots & \vdots\\ 1 & 6\end{bmatrix}\] 式 (2.1.15) 给出了随机模型。根据式 (2.2.24) \[\hat{\bm{X}}_{LS}=\left(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H}\right)^{-1}\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{Z}\] 得到 \[\hat{\bm{X}}_{LS}=\begin{bmatrix}\hat{\alpha}\\ \hat{\beta}\end{bmatrix}=\begin{bmatrix}2.23\\ 1.20\end{bmatrix}\mathrm{m} \tag{2.2.91}\] 估计得到的质点轨迹为 \[\hat{Z}=\hat{\alpha}+t\hat{\beta} \tag{2.2.92}\] 如图 2.7 所示，图中实线为质点的实际轨迹；圆点为观测值；虚线即为估计的质点轨迹。

\(\hat{\bm{X}}_{LS}\) 的验前方差矩阵为 \[\mathrm{Var}(\hat{\bm{X}}_{LS})=(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} =\begin{bmatrix}\sigma_{\hat{\alpha}}^2 & \sigma_{\hat{\alpha}}\sigma_{\hat{\beta}}\\ \sigma_{\hat{\alpha}}\sigma_{\hat{\beta}} & \sigma_{\hat{\beta}}^2\end{bmatrix} =\begin{bmatrix}0.1842 & -0.046\\ -0.046 & 0.0133\end{bmatrix}\mathrm{m}^2 \tag{2.2.93}\] 参数 \(\hat{\alpha}\) 和 \(\hat{\beta}\) 的验前中误差分别为 \(\sigma_{\hat{\alpha}}=0.43\,\mathrm{m}\) 和 \(\sigma_{\hat{\beta}}=0.16\,\mathrm{m}\)。在这里，由于假设已知随机模型（式 (2.1.15)），所以验前的单位权方差 \(\sigma_0^2=1\,\mathrm{m}^2\)。

以上得到的是参数估计 \(\left[\begin{array}{ll}\hat{\alpha} & \hat{\beta}\end{array}\right]^{\mathrm{T}}\) 的验前精度，在精度评定的时候默认单位权方差为 \(\sigma_0^2=1\,\mathrm{m}^2\)。下面用观测值来估计验后精度。

将式 (2.2.91) 代入误差方程得到残差 \(\bm{v}=\left[\begin{array}{llllll}-0.73 & 0.09 & 0.77 & 0.23 & -0.97 & 0.13\end{array}\right]^{\mathrm{T}}\mathrm{m}\)，验后单位权方差 \(\hat{\sigma}_0^2\) 为：

\[\begin{aligned} \hat{\sigma}_0^2&=\frac{\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}}{\ell-n}\\ &=9.9722/4\\ &=2.4930\ \mathrm{m}^2 \end{aligned}\] 最后， \[\hat{\sigma}_0=1.58\ \mathrm{m}\] 构造假设检验量 \[t=\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}=9.9722\] 若观测值无粗差，那么 \[t=\frac{\hat{\sigma}_0^2\cdot(\ell-n)}{\sigma_0^2}=\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}\sim\chi^2_{(\ell-n,\ \lambda=0)}\] 给定显著性水平 \(\alpha=0.01\)，自由度为 \(4\)，那么得到卡方分布的分位置 \(T=13.2767\)。可以看到 \(t<T\)，即认为在 \(99\%\) 的置信区间上，此例的观测值中没有显著的粗差。

最后，估计参数的验后估计精度为 \[\mathrm{Var}(\hat{\bm{X}}_{LS})=\frac{\hat{\sigma}_0^2}{\sigma_0^2}(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} =\begin{bmatrix}0.2910 & -0.0729\\ -0.0729 & 0.0210\end{bmatrix}\mathrm{m}^2\]

例 2.6表 2.3 给出了例 2.2 中 GPS 信号发送时的卫星坐标（WGS84）\((X^{s_i}\quad Y^{s_i}\quad Z^{s_i})\)（\(i=1,\ 2,\ \cdots,\ \ell\)）、观测值 \(\bm{Z}=\left[\begin{array}{llll}\rho_1 & \rho_2 & \cdots & \rho_{\ell}\end{array}\right]\) 和观测值中误差 \(\sigma_{\rho_i}\)。假设各观测值间相互独立，用最小二乘法估计接收机的坐标 \(\left(\begin{array}{lll}X_r & Y_r & Z_r\end{array}\right)\)。

GPS 测距定位的原理是距离交汇，但由于观测误差的存在，如图 2.8 所示，各卫星的距离观测值并不相交于用户所在的位置。最小二乘估计的目的是对有误差的观测值进行改正，改正后的观测值相交于一点，并满足残差加权平方和最小。

解：式 (2.1.16) 给出了观测值 \(\rho_i\) 与 \(\left(\begin{array}{lll}X_r & Y_r & Z_r\end{array}\right)\) 的函数关系，观测方程为： \[\rho_i=\left(\sqrt{\left(X^{s_i}-X_r\right)^2+\left(Y^{s_i}-Y_r\right)^2+\left(Z^{s_i}-Z_r\right)^2}+\tau\right)+\Delta_i\ ,\quad i=1,\ 2,\ \cdots,\ \ell\] 被估计参数为 \(\bm{X}=\left[\begin{array}{llll}X_r & Y_r & Z_r & \tau\end{array}\right]^{\mathrm{T}}\)，观测值数为 \(\ell=7\)，必要观测值数为 \(4\)，多余观测数为 \(3\)。显然，估计参数 \(\left(\begin{array}{lll}X_r & Y_r & Z_r\end{array}\right)\) 与观测值 \(\rho_i\) 呈非线性关系，所以需要先进行线性化，线性化后的观测方程为式 (2.1.38)： \[\begin{bmatrix}\Delta\rho_1\\ \Delta\rho_2\\ \vdots\\ \Delta\rho_7\end{bmatrix} =\begin{bmatrix}\rho_1-\rho_1^{*}\\ \rho_2-\rho_2^{*}\\ \vdots\\ \rho_7-\rho_7^{*}\end{bmatrix} =\begin{bmatrix} -\dfrac{\Delta X_1^{*}}{S_1^{*}} & -\dfrac{\Delta Y_1^{*}}{S_1^{*}} & -\dfrac{\Delta Z_1^{*}}{S_1^{*}} & 1\\[10pt] -\dfrac{\Delta X_2^{*}}{S_2^{*}} & -\dfrac{\Delta Y_2^{*}}{S_2^{*}} & -\dfrac{\Delta Z_2^{*}}{S_2^{*}} & 1\\[10pt] \vdots & \vdots & \vdots & \vdots\\[4pt] -\dfrac{\Delta X_7^{*}}{S_7^{*}} & -\dfrac{\Delta Y_7^{*}}{S_7^{*}} & -\dfrac{\Delta Z_7^{*}}{S_7^{*}} & 1 \end{bmatrix} \begin{bmatrix}x\\ y\\ z\\ \Delta\tau\end{bmatrix} +\begin{bmatrix}\Delta_1\\ \Delta_2\\ \vdots\\ \Delta_7\end{bmatrix}\] 现在取参数 \(\bm{X}\) 的近似值 \[\bm{X}^{*(1)}=\left[\begin{array}{llll}0 & 0 & 0 & 0\end{array}\right]^{\mathrm{T}}\] 经第一次线性化后，观测方程为 \[\underset{7\times 1}{\bm{z}^{(1)}}=\underset{7\times 4}{\bm{H}^{(1)}}\ \underset{4\times 1}{\bm{x}}+\underset{7\times 1}{\bm{\Delta}}\] 观测向量 \(\bm{z}^{(1)}\) 为 \[\bm{z}^{(1)}=\left[\begin{array}{lllllll}\Delta\rho_1 & \Delta\rho_2 & \cdots & \Delta\rho_7\end{array}\right]^{\mathrm{T}}\] 设计矩阵 \(\bm{H}^{(1)}\) 及第一次迭代估计得到的改正数见表 2.4“第一次迭代结果”。

现在取参数 \(\bm{X}=\left[\begin{array}{llll}X_r & Y_r & Z_r & \tau\end{array}\right]^{\mathrm{T}}\) 的近似值为 \(\bm{X}^{*(1)}=\left[\begin{array}{llll}0 & 0 & 0 & 0\end{array}\right]^{\mathrm{T}}\)，将观测值、卫星坐标和近似值代入式 (2.1.38)，得到第一次线性化后 \(\left[\begin{array}{llll}x & y & z & \Delta\tau\end{array}\right]^{\mathrm{T}}\) 的系数矩阵 \(\bm{H}\) 和观测值 \(\Delta\bm{\rho}\)。第一次线性化系数矩阵 \(\bm{H}\)、观测值 \(\Delta\bm{\rho}\) 和估计结果见表 2.4。

根据表 2.3 给出的中误差构成对角方差矩阵 \(\bm{D}\)： \[\bm{D}=\mathrm{diag}\left(\left[\begin{array}{lllllll}1.139 & 2.732 & 1.651 & 0.978 & 1.364 & 1.620 & 3.893\end{array}\right]\right)\mathrm{m}^2\] 设单位权中方差为 \(\sigma_0^2=1\,\mathrm{m}^2\)，按照图 2.6 的迭代流程，当 \(\max\left|\hat{\bm{x}}^{(k)}\right|<1\,\mathrm{m}\) 时停止迭代。程序共进行了 5 次迭代，第 5 次迭代的结果见表 2.4。

终止迭代后，输出最后一次迭代的结果，并进行精度评定。参数的验前方差为

参数的验前方差为 \[\mathrm{Var}(\hat{\bm{X}}_{LS})=(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} =\begin{bmatrix} 2.619401 & -2.271051 & -1.617138 & -2.487159\\ -2.271051 & 6.785561 & 4.615868 & 5.258552\\ -1.617138 & 4.615868 & 5.179831 & 4.284823\\ -2.487159 & 5.258552 & 4.284823 & 4.983638 \end{bmatrix}\mathrm{m}^2\] 接收机三维坐标的中误差分别为：\(\sigma_{\hat{X}}=1.62\,\mathrm{m}\)，\(\sigma_{\hat{Y}}=2.60\,\mathrm{m}\)，\(\sigma_{\hat{Z}}=2.27\,\mathrm{m}\)，接收机钟差的中误差为 \(\sigma_{\widehat{\delta t}}=2.23\,\mathrm{m}\)。

验后的单位权方差为 \[\hat{\sigma}_0^2=\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}/(7-4)=1.758\,\mathrm{m}^2\] 验后估计方差 \[\hat{\sigma}_0^2(\bm{H}^{\mathrm{T}}\bm{D}^{-1}\bm{H})^{-1} =\begin{bmatrix} 4.6049 & -3.9925 & -2.8429 & -4.3724\\ -3.9925 & 11.9290 & 8.1147 & 9.2445\\ -2.8429 & 8.1147 & 9.1061 & 7.5327\\ -4.3724 & 9.2445 & 7.5327 & 8.7612 \end{bmatrix}\mathrm{m}^2\]

将残差代入式， \[\Delta\hat{\bm{\rho}}=\Delta\bm{\rho}+\bm{v}\] 进而可以得到估计改正后的观测值 \[\hat{\bm{\rho}}=\Delta\hat{\bm{\rho}}+\bm{\rho}^{*}\]

从迭代过程看，由于第一次迭代估计时，对接收机的位置没有任何了解，所以取参数的近似值为 \(\bm{X}^{*(1)}=\left[\begin{array}{llll}0 & 0 & 0 & 0\end{array}\right]^{\mathrm{T}}\)，得到的 \(\hat{\bm{x}}^{(1)}\) 的绝对值和 \(\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}\) 都非常大，这是因为 \(\bm{X}^{*(1)}\) 与用户实际位置相差很大。在接下来的迭代计算中，\(\left|\hat{\bm{x}}^{(k)}\right|\) 和 \(\bm{v}^{\mathrm{T}}\bm{D}^{-1}\bm{v}\) 的数值迅速减小，表明迭代计算快速收敛，共迭代了 5 次后终止了迭代。在此问题中选取 \(1\,\mathrm{m}\) 作为迭代计算的门限值，是因为根据经验已知伪距观测值的精度为米级，所以导航定位精度也在米级，选取更小的迭代门限值进行迭代是徒劳的。

在最小二乘估计前，如图 2.8 所示，由于观测误差的干扰，观测值并不交汇于一点。在估计后，改正后的观测值（如虚线所示）交于一点，这就是 \(\hat{\bm{X}}_{LS}\) 所在的位置。此外，验前估计方差与验后的方差略有不同。这里的验前单位权方差默认为 \(1\,\mathrm{m}^2\)，验后的单位权方差由观测值计算得到，所以验后的估计方差更能客观地体现估计精度。

需要注意的是，不是所有迭代计算都可以如此例题中将近似值取为零，当近似值取值偏差太大，不能满足收敛条件时，会出现迭代发散，这时就不能得到正确的估计了。所以在解决实际问题时，如果对待估计参数有所了解，如已知 GPS 用户在武汉市的某个位置，就可以取任何已知的武汉某一点的坐标作为近似值，这样不仅可以保证收敛，也可以减少迭代计算的次数。

例 2.7用附有约束条件的最小二乘法估计例 2.3 中 \(P\) 点的坐标。

解：解决此问题的数学模型如例 2.3 所述，共有 5 个观测值和 1 个约束条件，估计 \(P_1\) 和 \(P_2\) 平面坐标的必要观测值数为 3，所以多余观测值数为 2。设参数的近似值为 \[\begin{aligned} \bm{X}^{*}&=\left[\begin{array}{llll}X_{P_1}^{*} & Y_{P_1}^{*} & X_{P_2}^{*} & Y_{P_2}^{*}\end{array}\right]^{\mathrm{T}}\\ &=\left[\begin{array}{llll}79.00 & 121.00 & 120.00 & 52.00\end{array}\right]^{\mathrm{T}}\mathrm{m} \end{aligned}\] 近似值可以用任意从已知基站观测的两条距离观测值交汇得到，如 \(P_1\) 的近似值可以用 \(Z_{A_1P_1}\) 和 \(Z_{A_2P_1}\) 交汇得到，也可以用 \(Z_{A_2P_1}\) 和 \(Z_{A_3P_1}\) 交汇得到。

设验前的单位权方差为 \(\sigma_0^2=5\,\mathrm{cm}^2\)，那么观测值权矩阵为单位矩阵。

在例 2.4 中已将观测方程和约束条件进行了线性化，按照附有约束条件的最小二乘估计方法进行估计并迭代计算，当 \(\max(|\bm{x}|)<0.01\,\mathrm{m}\) 时停止迭代。本例共进行了三次迭代计算，现以表格（表 2.5 和表 2.6）的形式给出第一次和最后一次迭代计算时的线性化观测方程、线性化约束条件和估计结果。由于表 2.5 和表 2.6 只将计算结果表示到千分位，所以当显示为 \(0.000\) 时，表示数值接近于零，但实际结果并不等于零。

第三次（最后一次）迭代得到参数估计的协因数矩阵为 \[\bm{Q}_{\hat{\bm{X}}}=\left(\bm{N}^{-1}-\bm{N}^{-1}\bm{C}^{\mathrm{T}}\bm{N}_{CC}^{-1}\bm{C}\bm{N}^{-1}\right) =\begin{bmatrix} 0.880 & -0.243 & 0.020 & -0.035\\ -0.257 & 0.629 & -0.021 & 0.037\\ 0.064 & -0.112 & 1.917 & -1.257\\ -0.069 & 0.122 & -1.299 & 1.812 \end{bmatrix}\] 验后的协方差矩阵为 \[\mathrm{Var}(\hat{\bm{x}}_{LS})=\hat{\sigma}_0^2\,\bm{Q}_{\hat{\bm{X}}} =\begin{bmatrix} 0.001 & -0.000 & 0.000 & 0.000\\ -0.000 & 0.001 & -0.000 & 0.000\\ 0.000 & -0.0004 & 0.003 & -0.002\\ -0.000 & 0.000 & -0.002 & 0.003 \end{bmatrix}\] 从上式中可以得到 \(P_1\) 和 \(P_2\) 点平面坐标的验后估计精度为 \[\begin{aligned} \sigma_{\hat{X}_{P_1}}&=0.04\,\mathrm{m}\\ \sigma_{\hat{Y}_{P_1}}&=0.04\,\mathrm{m}\\ \sigma_{\hat{X}_{P_2}}&=0.06\,\mathrm{m}\\ \sigma_{\hat{Y}_{P_2}}&=0.06\,\mathrm{m} \end{aligned}\]

在本例的模型中，估计的参数一定满足约束条件 \[\sqrt{\left(\hat{X}_{P_1}-\hat{X}_{P_2}\right)^2+\left(\hat{Y}_{P_1}-\hat{Y}_{P_2}\right)^2}-80.50=0\] 读者可自行对此进行验证（忽略计算误差）。

递推最小二乘估计

本章 2.2 节介绍了最小二乘的基本原理，本节将介绍在最小二乘估计基础上的“递推”实现形式。

在现实应用中，观测量并不是一次采样完成的，而是在不同时间点上获得的，我们可以等所有的观测完成后，利用所有观测值一次对参数进行估计，这也被称为“批处理”的最小二乘估计。批处理方法占用了计算机的大量内存，不能实时对数据进行处理，解决这个问题的方法是使用实时数据对参数估计不断进行更新，即递推最小二乘的算法。递推最小二乘算法与最小二乘批处理的关系如图 2.9 所示。假定在第 \(k\) 次观测后，得到第 \(k\) 次和所有之前的观测值，组成了观测值向量 \(\bm{Z}_k=\left[\begin{array}{llll}\bm{z}_1^{\mathrm{T}} & \bm{z}_2^{\mathrm{T}} & \cdots & \bm{z}_k^{\mathrm{T}}\end{array}\right]^{\mathrm{T}}\)，用观测值向量 \(\bm{Z}_k\) 得到最小二乘估计 \(\hat{\bm{X}}_{B(k)}\)，下标 \(B\) 表示当前所有观测值批处理得到的最小二乘估计。在第 \(k+1\) 次观测后，又得到了一组观测值向量 \(\bm{z}_{k+1}\)，集合所有的观测值 \(\bm{Z}_{k+1}=\left[\begin{array}{ll}\bm{Z}_k^{\mathrm{T}} & \bm{z}_{k+1}^{\mathrm{T}}\end{array}\right]^{\mathrm{T}}\)，估计得到的参数估计 \(\hat{\bm{X}}_{B(k+1)}\)。按照批处理方法，每次得到新的观测值都需要重新集合所有的历史观测值进行参数估计，占用了较多的内存，而递推最小二乘只需要用新观测值 \(\bm{z}_{k+1}\) 对 \(\hat{\bm{X}}_{B(k)}\) 进行更新就得到了 \(\hat{\bm{X}}_{B(k+1)}\)。

与第 4 章介绍的 Kalman 滤波比较，可以看出递推的最小二乘估计是当状态方程“静止”情况下的 Kalman 滤波。下面推导如何用 \(\bm{z}_{k+1}\) 对已有的估计 \(\hat{\bm{X}}_{B(k)}\) 进行更新修正得到 \(\hat{\bm{X}}_{B(k+1)}\)。在推导中，为了简单起见，最小二乘的批处理结果不再带有下标 \(B\)。

递推最小二乘估计推导

1. 参数的递推估计

在第 \(k\) 次观测后，观测值向量 \(\bm{Z}_k\) 的观测方程为 \[\bm{Z}_k=\bm{H}_k\bm{X}+\bm{\Delta}_k \tag{2.3.1}\] 测值向量 \(\bm{Z}_k\) 的方差阵为 \(\bm{D}_k\)，权矩阵为 \(\bm{W}_k=\sigma_0^2\bm{D}_k^{-1}\)。最小二乘估计 \(\hat{\bm{X}}_k\) \[\hat{\bm{X}}_k=\left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k\right)^{-1}\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{Z}_k \tag{2.3.2}\] \(\hat{\bm{X}}_k\) 的协因数矩阵为 \[\bm{Q}_{\hat{\bm{X}}_k}=\left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k\right)^{-1} \tag{2.3.3}\] 残差向量为 \[\bm{V}_k=\bm{H}_k\hat{\bm{X}}_k-\bm{Z}_k \tag{2.3.4}\] 验后单位权中误差为 \[\hat{\sigma}_{0,\ k}^2=\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{V}_k/(\ell_k-n) \tag{2.3.5}\] 其中 \(\ell_k\) 是 \(\bm{Z}_k\) 中的观测值个数。

\(k+1\) 次观测值向量为 \(\bm{z}_{k+1}\)，观测方程为 \[\bm{z}_{k+1}=\bm{h}_{k+1}\bm{X}+\bm{\Delta}_{k+1} \tag{2.3.6}\] 观测值 \(\bm{z}_{k+1}\) 的方差阵为 \(\bm{d}_{k+1}\)，权矩阵为 \(\bm{w}_{k+1}=\sigma_0^2\bm{d}_{k+1}^{-1}\)。集合所有的观测值构成观测值向量 \(\bm{Z}_{k+1}=\left[\begin{array}{ll}\bm{Z}_k^{\mathrm{T}} & \bm{z}_{k+1}^{\mathrm{T}}\end{array}\right]^{\mathrm{T}}\)，\(\bm{Z}_{k+1}\) 的观测方程为 \[\bm{Z}_{k+1}=\begin{bmatrix}\bm{Z}_k\\ \bm{z}_{k+1}\end{bmatrix} =\begin{bmatrix}\bm{H}_k\\ \bm{h}_{k+1}\end{bmatrix}\bm{X} +\begin{bmatrix}\bm{\Delta}_k\\ \bm{\Delta}_{k+1}\end{bmatrix} \tag{2.3.7}\] 观测值向量 \(\bm{Z}_{k+1}\) 的权矩阵为 \[\bm{W}_{k+1}=\begin{bmatrix}\bm{W}_k & \\ & \bm{w}_{k+1}\end{bmatrix} \tag{2.3.8}\] 设 \[\bm{H}_{k+1}=\begin{bmatrix}\bm{H}_k\\ \bm{h}_{k+1}\end{bmatrix}\qquad \bm{Z}_{k+1}=\begin{bmatrix}\bm{Z}_k\\ \bm{z}_{k+1}\end{bmatrix} \tag{2.3.9}\] 那么，批处理最小二乘估计为 \[\begin{aligned} &\hat{\bm{X}}_{k+1}=\left(\bm{H}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{H}_{k+1}\right)^{-1} \bm{H}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{Z}_{k+1}\\ =&\left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k+\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1}\right)^{-1} \left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{Z}_k+\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{z}_{k+1}\right) \end{aligned} \tag{2.3.10}\] 残差向量为 \[\bm{V}_{k+1}=\begin{bmatrix}\overline{\bm{V}}_k\\ \bm{v}_{k+1}\end{bmatrix} =\begin{bmatrix}\bm{H}_k\\ \bm{h}_{k+1}\end{bmatrix}\hat{\bm{X}}_{k+1} -\begin{bmatrix}\bm{Z}_k\\ \bm{z}_{k+1}\end{bmatrix} \tag{2.3.11}\] 注意，这里的 \(\overline{\bm{V}}_k\) 与式 (2.3.4) 中的 \(\bm{V}_k\) 并不相同。

\(\hat{\bm{X}}_{k+1}\) 的协因数矩阵为 \[\begin{aligned} \bm{Q}_{\hat{\bm{X}}_{k+1}}&=\left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k+\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1}\right)^{-1}\\ &=\left(\bm{Q}_{\hat{\bm{X}}_k}^{-1}+\bm{h}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{h}_{k+1}\right)^{-1} \end{aligned} \tag{2.3.12}\] 式 (2.3.10) 可表示为 \[\hat{\bm{X}}_{k+1}=\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{Z}_k +\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{z}_{k+1} \tag{2.3.13}\] 由于 \[\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{Z}_k=\bm{Q}_{\hat{\bm{X}}_k}^{-1}\hat{\bm{X}}_k \tag{2.3.14}\] 式 (2.3.13) 为 \[\hat{\bm{X}}_{k+1}=\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{Q}_{\hat{\bm{X}}_k}^{-1}\hat{\bm{X}}_k +\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{z}_{k+1} \tag{2.3.15}\] 从式 (2.3.12) 可以得到 \[\bm{Q}_{\hat{\bm{X}}_k}^{-1}=\bm{Q}_{\hat{\bm{X}}_{k+1}}^{-1}-\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1} \tag{2.3.16}\] 将式 (2.3.16) 代入式 (2.3.15) \[\begin{aligned} \hat{\bm{X}}_{k+1}&=\bm{Q}_{\hat{\bm{X}}_{k+1}}\left(\bm{Q}_{\hat{\bm{X}}_{k+1}}^{-1}-\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1}\right)\hat{\bm{X}}_k +\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{z}_{k+1}\\ &=\hat{\bm{X}}_k-\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1}\hat{\bm{X}}_k +\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{z}_{k+1}\\ &=\hat{\bm{X}}_k+\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1} \left(\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k\right) \end{aligned} \tag{2.3.17}\] 令 \[\bm{K}_{k+1}=\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1} \tag{2.3.18}\] 那么式 (2.3.17) 为 \[\hat{\bm{X}}_{k+1}=\hat{\bm{X}}_k+\bm{K}_{k+1}\left(\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k\right) \tag{2.3.19}\] 若设 \[\Delta\bm{z}_{k+1}=\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k \tag{2.3.20}\] 式 (2.3.19) 成为 \[\hat{\bm{X}}_{k+1}=\hat{\bm{X}}_k+\bm{K}_{k+1}\Delta\bm{z}_{k+1} \tag{2.3.21}\]

式 (2.3.19) 是从批处理结果式 (2.3.10) 推导得到，它将 \(\hat{\bm{X}}_{k+1}\) 表示为观测值 \(\bm{z}_{k+1}\) 对 \(\hat{\bm{X}}_k\) 的更新，其中 \(\bm{K}_{k+1}\) 将 \(\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k\) 映射到 \(\hat{\bm{X}}_{k+1}\)，被称为增益矩阵。由于增益矩阵 \(\bm{K}_{k+1}\) 由 \(\bm{Q}_{\hat{\bm{X}}_{k+1}}\) 计算得到，这给计算带来不便，下面推导由 \(\bm{Q}_{\hat{\bm{X}}_k}\) 和 \(\bm{w}_{k+1}\) 来计算 \(\bm{Q}_{\hat{\bm{X}}_{k+1}}\) 和 \(\bm{K}_{k+1}\) 的递推式。

将式 (2.3.12) 代入式 (2.3.18) \[\begin{aligned} \bm{K}_{k+1}&=\bm{Q}_{\hat{\bm{X}}_{k+1}}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\\ &=\left(\bm{Q}_{\hat{\bm{X}}_k}^{-1}+\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{h}_{k+1}\right)^{-1}\bm{h}_{k+1}^{\mathrm{T}}\bm{w}_{k+1} \end{aligned} \tag{2.3.22}\] 由矩阵的恒等式 (A-53) 可得到 \[\bm{K}_{k+1}=\bm{Q}_{\hat{\bm{X}}_k}\bm{h}_{k+1}^{\mathrm{T}} \left(\bm{w}_{k+1}^{-1}+\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k}\bm{h}_{k+1}^{\mathrm{T}}\right)^{-1} \tag{2.3.23}\] 上式即为由 \(\bm{Q}_{\hat{\bm{X}}_k}\) 和 \(\bm{w}_{k+1}\) 计算 \(\bm{K}_{k+1}\) 的递推式。

由附录中式 (A-52)，式 (2.3.12) 可表示为 \[\bm{Q}_{\hat{\bm{X}}_{k+1}}=\bm{Q}_{\hat{\bm{X}}_k}-\bm{Q}_{\hat{\bm{X}}_k}\bm{h}_{k+1}^{\mathrm{T}} \left(\bm{w}_{k+1}^{-1}+\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_l}\bm{h}_{k+1}^{\mathrm{T}}\right)^{-1}\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k} \tag{2.3.24}\] 根据式 (2.3.23)，式 (2.3.24) 为 \[\bm{Q}_{\hat{\bm{X}}_{k+1}}=\bm{Q}_{\hat{\bm{X}}_k}-\bm{K}_{k+1}\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k} \tag{2.3.25}\] 上式即为由 \(\bm{Q}_{\hat{\bm{X}}_k}\)、\(\bm{h}_{k+1}\) 和 \(\bm{w}_{k+1}\) 计算 \(\bm{Q}_{\hat{\bm{X}}_{k+1}}\) 的递推式。在得到 \(\bm{Q}_{\hat{\bm{X}}_{k+1}}\) 后，进而得到 \[\mathrm{Var}\left(\hat{\bm{X}}_{k+1}\right)=\sigma_0^2\left(\bm{Q}_{\hat{\bm{X}}_k}-\bm{K}_{k+1}\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k}\right) \tag{2.3.26}\] 式 (2.3.19)、(2.3.23)、(2.3.25) 和式 (2.3.26) 一起构成了对 \(\hat{\bm{X}}_{k+1}\) 和 \(\mathrm{Var}\left(\hat{\bm{X}}_{k+1}\right)\) 的递推公式。

2. 验后单位权方差的递推

将估计得到的 \(\hat{\bm{X}}_{k+1}\) 代入误差方程 (2.3.11) 得到残差 \(\bm{V}_{k+1}\)，从而可以计算验后单位权方差 \[\hat{\sigma}_{0,\ k+1}^2=\frac{\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1}}{\ell_{k+1}-n} \tag{2.3.27}\] 这里的 \(\ell_{k+1}\) 为观测值向量 \(\bm{Z}_{k+1}=\left[\begin{array}{ll}\bm{Z}_k^{\mathrm{T}} & \bm{z}_{k+1}^{\mathrm{T}}\end{array}\right]^{\mathrm{T}}\) 的观测值的个数。上式中的残差加权平方和为 \[\begin{aligned} &\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1} =\begin{bmatrix}\overline{\bm{V}}_k^{\mathrm{T}} & \bm{v}_{k+1}^{\mathrm{T}}\end{bmatrix} \begin{bmatrix}\bm{W}_k & \\ & \bm{w}_{k+1}\end{bmatrix} \begin{bmatrix}\overline{\bm{V}}_k\\ \bm{v}_{k+1}\end{bmatrix}\\ =&\overline{\bm{V}}_k^{\mathrm{T}}\bm{W}_k\overline{\bm{V}}_k+\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1} \end{aligned} \tag{2.3.28}\] 根据式 (2.3.11)、式 (2.3.19) 和式 (2.3.4) 得到 \[\begin{aligned} \overline{\bm{V}}_k&=\bm{H}_k\hat{\bm{X}}_{k+1}-\bm{Z}_k\\ &=\bm{H}_k\left(\hat{\bm{X}}_k+\bm{K}_{k+1}\Delta\bm{z}_{k+1}\right)-\bm{Z}_k\\ &=\bm{V}_k+\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1} \end{aligned} \tag{2.3.29}\] 将式 (2.3.29) 代入式 (2.3.28) \[\begin{aligned} &\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1} =\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{V}_k+\left(\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1}\right)^{\mathrm{T}}\bm{W}_k\bm{V}_k\\ &\qquad+\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1} +\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_k}\bm{K}_{k+1}\Delta\bm{z}_{k+1} +\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1} \end{aligned} \tag{2.3.30}\] 将式 (2.3.4) 代入上式中 \(\left(\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1}\right)^{\mathrm{T}}\bm{W}_k\bm{V}_k\)，得到 \[\begin{aligned} \left(\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1}\right)^{\mathrm{T}}\bm{W}_k\bm{V}_k &=\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{H}_k^{\mathrm{T}}\bm{W}_k\left(\bm{H}_k\hat{\bm{X}}_k-\bm{Z}_k\right)\\ &=\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\left(\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k\hat{\bm{X}}_k-\bm{H}_k^{\mathrm{T}}\bm{W}_k\bm{Z}_k\right) \end{aligned} \tag{2.3.31}\] 由式 (2.3.2) 可知，上式等于零，同理，\(\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{H}_k\bm{K}_{k+1}\Delta\bm{z}_{k+1}\) 也等于零，所以 \[\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1} =\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{V}_k+\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_k}\bm{K}_{k+1}\Delta\bm{z}_{k+1} +\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1} \tag{2.3.32}\] 其中 \(\bm{V}_k^{\mathrm{T}}\bm{W}_k\bm{V}_k\) 是由观测值向量 \(\bm{Z}_k\) 估计得到的残差加权平方和。将式 (2.3.5) 代入上式 \[\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1} =\hat{\sigma}_{0,\ k}^2(\ell_k-n)+\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_k}\bm{K}_{k+1}\Delta\bm{z}_{k+1} +\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1} \tag{2.3.33}\] 将上式代入式 (2.3.27)，得到 \[\begin{aligned} \hat{\sigma}_{0,\ k+1}^2&=\frac{1}{(\ell_{k+1}-n)}\bm{V}_{k+1}^{\mathrm{T}}\bm{W}_{k+1}\bm{V}_{k+1}\\ &=\frac{1}{(\ell_{k+1}-n)}\left[\hat{\sigma}_{0,\ k}^2(\ell_k-n)+\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_k}\bm{K}_{k+1}\Delta\bm{z}_{k+1} +\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1}\right] \end{aligned} \tag{2.3.34}\] 式 (2.3.34) 即为递推最小二乘估计的验后单位权方差。

计算流程

递推的最小二乘算法流程如图 2.10 所示。。在时刻 \(t_1\)，有的观测值向量 \(\bm{Z}_1\)，观测方程为 \[\bm{Z}_1=\bm{H}_1\bm{X}+\bm{\Delta}_1 \tag{2.3.35}\] 观测值向量 \(\bm{Z}_1\) 的方差阵为 \(\bm{D}_1\)，权矩阵为 \(\bm{W}_1=\sigma_0^2\bm{D}_1^{-1}\)，这时参数的最小二乘解为 \[\hat{\bm{X}}_1=\left(\bm{H}_1^{\mathrm{T}}\bm{W}_1\bm{H}_1\right)^{-1}\bm{H}_1^{\mathrm{T}}\bm{W}_1\bm{Z}_1 \tag{2.3.36}\] 协因数矩阵为 \[\bm{Q}_{\hat{\bm{X}}_1}=\left(\bm{H}_1^{\mathrm{T}}\bm{W}_1\bm{H}_1\right)^{-1} \tag{2.3.37}\] 残差向量为 \[\bm{V}_1=\bm{h}_1\hat{\bm{X}}_1-\bm{Z}_1 \tag{2.3.38}\] 然后，验后单位权中误差 \[\hat{\sigma}_{0,\ 1}^2=\frac{1}{(\ell_1-n)}\bm{V}_1^{\mathrm{T}}\bm{W}_1\bm{V}_1 \tag{2.3.39}\] 保存 \(\hat{\bm{X}}_1\)，\(\bm{Q}_{\hat{\bm{X}}_1}\) 和 \(\hat{\sigma}_{0,\ 1}^2\)，接下来，在 \(t_{k+1}\)（\(k=1,\ 2,\ \cdots,\ N\)）时刻获得新的观测值 \(\bm{z}_{k+1}\)，观测值 \(\bm{z}_{k+1}\) 的方差阵为 \(\bm{d}_{k+1}\)，权矩阵为 \(\bm{w}_{k+1}=\sigma_0^2\bm{d}_{k+1}^{-1}\)，对 \(k\) 时刻的估计进行更新 \[\begin{aligned} \hat{\bm{X}}_{k+1}&=\hat{\bm{X}}_k+\bm{K}_{k+1}\Delta\bm{z}_{k+1} \tag{2.3.40}\\ \bm{K}_{k+1}&=\bm{Q}_{\hat{\bm{X}}_k}\bm{h}_{k+1}^{\mathrm{T}}\left(\bm{w}_{k+1}^{-1}+\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k}\bm{h}_{k+1}^{\mathrm{T}}\right)^{-1} \tag{2.3.41}\\ \Delta\bm{z}_{k+1}&=\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k \tag{2.3.42}\\ \bm{Q}_{\hat{\bm{X}}_{k+1}}&=\bm{Q}_{\hat{\bm{X}}_k}-\bm{K}_{k+1}\bm{h}_{k+1}\bm{Q}_{\hat{\bm{X}}_k} \tag{2.3.43} \end{aligned}\] \[\hat{\sigma}_{0,\ k+1}^2=\frac{1}{(\ell_{k+1}-n)}\left[\hat{\sigma}_{0,\ k}^2(\ell_k-n)+\Delta\bm{z}_{k+1}^{\mathrm{T}}\bm{K}_{k+1}^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_k}\bm{K}_{k+1}\Delta\bm{z}_{k+1}+\bm{v}_{k+1}^{\mathrm{T}}\bm{w}_{k+1}\bm{v}_{k+1}\right] \tag{2.3.44}\] \[\mathrm{Var}\left(\hat{\bm{X}}_{k+1}\right)=\bm{Q}_{k+1}\hat{\sigma}_{0,\ k+1}^2 \tag{2.3.45}\] 在 \(t_{k+2}\) 时刻获得 \(\bm{z}_{k+2}\) 后，\(\bm{z}_{k+2}\) 再对 \(t_{k+1}\) 的估计结果进行更新，如此进行下去，直到所有观测值对参数更新完毕。

若观测值是参数的非线性函数 \[\bm{z}_{k+1}=\bm{f}_{k+1}(\bm{X}_{k+1})+\bm{\Delta}_{k+1} \tag{2.3.46}\] 首先应该将其线性化 \[\bm{z}_{k+1}=\bm{f}_{k+1}\left(\bm{X}_{k+1}^{*}\right)+\bm{h}_{k+1}\left(\bm{X}_{k+1}-\bm{X}_{k+1}^{*}\right)+\bm{\Delta}_{k+1} \tag{2.3.47}\] \(\bm{h}_{k+1}\) 为 \[\bm{h}_{k+1}=\left[\left(\frac{\partial\bm{f}_{k+1}}{\partial\bm{X}_{k+1}}\right)\right]_{\bm{X}_{k+1}=\bm{X}_{k+1}^{*}} \tag{2.3.48}\] 式 (2.3.47) 可以表示为 \[\bm{z}_{k+1}-\left[\,\bm{f}_{k+1}\left(\bm{X}_{k+1}^{*}\right)-\bm{h}_{k+1}\bm{X}_{k+1}^{*}\,\right]=\bm{h}_{k+1}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \tag{2.3.49}\] 记 \[\bm{z}'_{k+1}=\bm{z}_{k+1}-\left[\,\bm{f}_{k+1}\left(\bm{X}_{k+1}^{*}\right)-\bm{h}_{k+1}\bm{X}_{k+1}^{*}\,\right] \tag{2.3.50}\] 得到线性化后的观测方程为 \[\bm{z}'_{k+1}=\bm{h}_{k+1}\bm{X}_{k+1}+\bm{\Delta}_{k+1} \tag{2.3.51}\] 这里需要注意的是，在 \(t_{k+1}\) 时刻的每一次迭代中都要更新 \(\bm{h}_{k+1}\)、\(\bm{Q}_{\hat{\bm{X}}_{k+1}}\)、\(\bm{K}_{k+1}\) 和 \(\Delta\bm{z}_{k+1}\)。其中 \(\Delta\bm{z}_{k+1}\) 的迭代计算为 \[\Delta\bm{z}_{k+1}=\bm{z}'_{k+1}-\bm{h}_{k+1}\bm{X}_{k+1}^{*} \tag{2.3.52}\] 将式 (2.3.50) 代入式 (2.3.52)，得到 \[\Delta\bm{z}_{k+1}=\bm{z}_{k+1}-\bm{f}_{k+1}\left(\bm{X}_{k+1}^{*}\right) \tag{2.3.53}\] 在第一次迭代可以取近似值为 \(\hat{\bm{X}}_k\)，即 \(\bm{X}_{k+1}^{*}=\hat{\bm{X}}_k\)。从第二次迭代开始，每次取上一次迭代得到的参数估计作为近似值，直到 \(\hat{\bm{X}}_{k+1}\) 没有明显的变化，终止迭代。

从上面的计算过程可以看出，每次最小二乘估计只需存储估计结果 \(\hat{\bm{X}}_k\)、\(\bm{Q}_{\hat{\bm{X}}_k}\) 和 \(\hat{\sigma}_{0,\ k}^2\)，当有新的观测值 \(\bm{z}_{k+1}\) 后，无需存储旧的观测值，只需对 \(\hat{\bm{X}}_k\)、\(\bm{Q}_{\hat{\bm{X}}_k}\) 和 \(\hat{\sigma}_{0,\ k}^2\) 进行更新计算，这给大样本的观测数据处理带来了便利。

在递推的最小二乘估计中，若将 \(\hat{\bm{X}}_k\) 看作是 \(t_{k+1}\) 时刻参数估计的先验信息，\(\sigma_0^2\bm{Q}_{\hat{\bm{X}}_k}\) 为先验方差，那么上面的递推实际上是 \(t_{k+1}\) 时刻的观测值 \(\bm{z}_{k+1}\) 对先验信息的更新，所以递推的最小二乘也被称为“最小二乘滤波”。如果被估计目标是静态的，这样获得的参数估计与将集合所有观测对参数进行一次估计的结果等价；如果系统是运动状态的，那么每一次更新得到的 \(\hat{\bm{X}}_k\) 就是对当前状态的估计值。

算例分析

例 2.8为了估计例 2.5 中匀速运动质点的轨迹，除了已有的 6 个观测值外（第一期观测值见表 2.1），现又观测了如表 2.7 的观测值（第二期观测值）。请利用表 2.1 和表 2.7 的所有观测值重新估计运动质点的轨迹。

解：第一期观测值有 6 个，参数有 2 个，多余观测数为 \(\ell_1-n=4\)。第二期又增加了 3 个观测值，可以构成 9 个观测方程估计运动质点的轨迹，多余观测有 \(\ell_2-n=7\)。

现用本节中介绍的递推最小二乘算法来解决此问题。在例 2.5 中用已有 6 个观测值估计了轨迹参数为 \[\hat{\bm{X}}_1=\begin{bmatrix}\hat{\alpha}\\ \hat{\beta}\end{bmatrix}=\begin{bmatrix}2.23\\ 1.20\end{bmatrix}\mathrm{m}\] 协因数矩阵为 \[\bm{Q}_{\hat{\bm{X}}_1}=\begin{bmatrix}0.1842 & -0.046\\ -0.046 & 0.0133\end{bmatrix}\] 验后单位权中误差为 \(\hat{\sigma}_{0,1}=1.6\,\mathrm{m}\)。

新增观测值的观测方程为 \[\bm{z}_2=\bm{h}_2\bm{X}+\bm{\Delta}_2\] 其中 \[\bm{h}_2=\begin{bmatrix}1 & 7\\ 1 & 8\\ 1 & 9\end{bmatrix}\qquad \bm{z}_2=\begin{bmatrix}11.29\\ 12.81\\ 13.99\end{bmatrix}\] 单位权方差与第一期观测值进行最小二乘估计的单位权方差一致，\(\sigma_0^2=1\)，那么新增观测值的权阵为 \[\bm{w}_2=1\times\begin{bmatrix}0.16 & 0 & 0\\ 0 & 0.16 & 0\\ 0 & 0 & 0.16\end{bmatrix}^{-1}\]

现用新增观测值对参数估计进行更新修正。首先计算 \[\Delta\bm{z}_2=\bm{z}_2-\bm{h}_2\hat{\bm{X}}_1=\begin{bmatrix}0.70\\ 1.03\\ 1.01\end{bmatrix}\] 将 \(\bm{Q}_{\hat{\bm{X}}_1}\)，\(\bm{w}_2\) 和 \(\bm{h}_2\) 代入式 (2.3.41)，得到增益矩阵 \[\bm{K}_2=\begin{bmatrix}-0.0441 & -0.1416 & -0.2391\\ 0.0340 & 0.0528 & 0.0715\end{bmatrix}\] 由式 (2.3.40) 得到 \[\hat{\bm{X}}_2=\begin{bmatrix}1.81\\ 1.34\end{bmatrix}\mathrm{m}\] 利用式 (2.3.43)，得到 \[\bm{Q}_{\hat{\bm{X}}_2}=\begin{bmatrix}0.1021 & -0.0156\\ -0.0156 & 0.0030\end{bmatrix}\] \(\bm{z}_2\) 的残差为 \[\begin{aligned} \bm{v}_2&=\bm{h}_2\hat{\bm{X}}_2-\bm{z}_2\\ &=\begin{bmatrix}v_7\\ v_8\\ v_9\end{bmatrix}=\begin{bmatrix}-0.07\\ -0.24\\ -0.07\end{bmatrix}\mathrm{m} \end{aligned}\] 残差平方和为 \[\bm{v}_2^{\mathrm{T}}\bm{w}_2\bm{v}_2=0.4305\] 将以上结果代入式 (2.3.44)，得到验后单位权方差 \[\begin{aligned} \hat{\sigma}_{0,\ 2}^2&=\frac{1}{(\ell_2-n)}\left[\hat{\sigma}_{0,\ 1}^2(\ell_1-n)+\Delta\bm{z}_2^{\mathrm{T}}\bm{K}_2^{\mathrm{T}}\bm{Q}_{\hat{\bm{X}}_1}\bm{K}_2\Delta\bm{z}_2+\bm{v}_2^{\mathrm{T}}\bm{w}_2\bm{v}_2\right]\\ &=12.2929/7\\ &=1.75\,\mathrm{m}^2 \end{aligned}\] 验后单位权方差估计为 \[\hat{\sigma}_{0,\ 2}=1.4\,\mathrm{m}\] 在得到验后的验后单位权方差后，代入式 (2.3.45)，即可计算 \(\hat{\bm{X}}_2\) 的验后方差。

在此例中，新增观测值通过递推最小二乘算法对已有的参数估计 \(\hat{\bm{X}}_1\) 进行更新重新得到的参数估计 \(\hat{\bm{X}}_2\)、方差矩阵和验后单位权方差估计，结果与集合所有观测值（批处理）进行估计的结果一致，读者可自行进行计算验证。

极大似然估计

极大似然估计是遗传学家和统计学家罗纳德·费雪爵士（Ronald Aylmer Fisher）在 1912 年至 1922 年间提出并开始使用的。

极大似然估计

极大似然估计提供了一种给定观测值来评估模型参数的方法，即“模型已定，参数未知”。例如，我们已经知道观测值服从某一分布 \(p(\bm{z}\mid\bm{x})\)，其中 \(\bm{z}=\left[\begin{array}{llll}z_1, & z_2, & \cdots, & z_{\ell}\end{array}\right]^{\mathrm{T}}\) 和 \(\bm{x}=\left[\begin{array}{llll}x_1, & x_2, & \cdots, & x_n\end{array}\right]^{\mathrm{T}}\)，如观测值服从正态分布：\(\bm{Z}\sim N(\bm{\mu}_Z,\ \bm{D})\)，但是该分布函数中的参数 \(\bm{\mu}_Z\) 和 \(\bm{D}\) 未知，就可以通过采样，即观测值 \(\bm{Z}=\left[\begin{array}{llll}Z_1, & Z_2, & \cdots, & Z_{\ell}\end{array}\right]\)，来求分布函数中的未知参数 \(\bm{\mu}_Z\) 和 \(\bm{D}\)。为了能估计得到分布中的未知参数 \(\bm{x}=\left[\begin{array}{llll}x_1, & x_2, & \cdots, & x_n\end{array}\right]^{\mathrm{T}}\)，一个合理想法是参数的估值 \(\hat{x}_1,\ \hat{x}_2,\ \cdots,\ \hat{x}_n\) 使得概率密度函数 \(p(\bm{z}\mid\bm{x})\) 最大 \[p(\bm{z}\mid\hat{\bm{x}})=\max \tag{2.4.1}\] \(p(\bm{z}\mid\bm{x})\) 也称为似然函数。记满足使似然函数最大的估计为 \(\hat{\bm{X}}_{ML}\)，即 \[\hat{\bm{X}}_{ML}=\arg\max_{\hat{\bm{x}}}p(\bm{z}\mid\hat{\bm{x}}) \tag{2.4.2}\] 极大似然估计可以理解为：“在什么样的状态下，最可能产生现在的观测数据”，在这样的准则下，得到的参数估计是最符合观测值的估计。Fisher 学派认为参数 \(\bm{x}\) 虽然未知，但是固定的常数，为非随机量，此时的 \(p(\bm{z}\mid\bm{x})\) 是观测值 \(\bm{z}\) 的概率密度函数，所以 \[p(\bm{z}\mid\bm{x})=p_z(\bm{z}) \tag{2.4.3}\] 在这里为了突出 \(\bm{x}\) 是未知待求解的参数，仍然用 \(p(\bm{z}\mid\bm{x})\) 来表示的 \(\bm{z}\) 分布。

为了求得似然函数的最大值，对似然函数求导并求解 \[\left.\frac{\partial p(\bm{z}\mid\bm{x})}{\partial\bm{x}}\right|_{\bm{x}=\hat{\bm{X}}_{ML}}=0 \tag{2.4.4}\] 得到 \(\hat{\bm{X}}_{ML}\)。由于很多时候似然函数含有指数函数，为了方便解得 \(\hat{\bm{X}}_{ML}\)，可先对似然函数求对数，然后求导解得 \(\hat{\bm{X}}_{ML}\) \[\left.\frac{\partial\ln p(\bm{z}\mid\bm{x})}{\partial\bm{x}}\right|_{\bm{x}=\hat{\bm{X}}_{ML}}=0 \tag{2.4.5}\] 上式中的 \(\ln p(\bm{z}\mid\bm{x})\) 称为对数似然函数。

从以上过程来看，极大似然估计 \(\hat{\bm{X}}_{ML}\) 就是使似然函数 \(p(\bm{z}\mid\bm{x})\) 或者对数似然函数 \(\ln p(\bm{z}\mid\bm{x})\) 最大的估计。

算例分析

例 2.9设观测值 \(\bm{Z}=\left[\begin{array}{llll}Z_1, & Z_2, & \cdots, & Z_n\end{array}\right]\) 相互独立，\(Z_i(i=1,\ 2,\ \cdots,\ n)\) 服从 \(N(\mu,\ \sigma^2)\)，其中 \(\mu\)，\(\sigma^2\) 未知。试求 \(\mu\)，\(\sigma^2\) 的极大似然估计。

解：由于观测值相互独立，\(\left[\begin{array}{llll}Z_1, & Z_2, & \cdots, & Z_n\end{array}\right]\) 的联合分布为各个观测值分布函数的乘积，所以似然函数为： \[\begin{aligned} L(\mu,\ \sigma^2)&=\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(z_i-\mu)^2}{2\sigma^2}}\\ &=(2\pi\sigma^2)^{-n/2}e^{-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(z_i-\mu)^2} \end{aligned} \tag{2.4.6}\] 其对数似然函数为 \[l(\mu,\ \sigma^2)=-\frac{n}{2}\ln(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(z_i-\mu)^2 \tag{2.4.7}\] 由极值条件，有 将 \(l(\mu,\ \sigma^2)\) 分别对 \(\mu\)、\(\sigma^2\) 求偏导，并令它们都为 \(0\)，得似然方程组为： \[\begin{cases} \dfrac{\partial l(\mu,\ \sigma^2)}{\partial\mu}=\dfrac{1}{\sigma^2}\displaystyle\sum_{i=1}^{n}(z_i-\mu)=0\\[10pt] \dfrac{\partial l(\mu,\ \sigma^2)}{\partial\sigma^2}=-\dfrac{n}{2\sigma^2}+\dfrac{1}{2\sigma^4}\displaystyle\sum_{i=1}^{n}(z_i-\mu)^2=0 \end{cases} \tag{2.4.8}\] 解似然方程组得： \[\begin{aligned} \hat{\mu}_{ML}&=\frac{1}{n}\sum_{i=1}^{n}z_i=\overline{z}\\ \hat{\sigma}_{ML}^2&=\frac{1}{n}\sum_{i=1}^{n}(z_i-\overline{z})^2 \end{aligned} \tag{2.4.9}\] 上述过程对一切样本成立，故用观测值 \(Z_i\) 代替 \(z_i\)，\(\mu\) 和 \(\sigma^2\) 的极大似然估计分别为： \[\hat{\mu}_{ML}=\overline{Z}\ ,\quad\hat{\sigma}_{ML}^2=\frac{1}{n}\sum_{i=1}^{n}(Z_i-\overline{Z})^2 \tag{2.4.10}\] 上述参数估计的期望为 \[E(\hat{\mu}_{ML})=\frac{1}{n}\sum_{k=1}^{n}E(Z_k)=\mu \tag{2.4.11}\]

\[E(\hat{\sigma}_{ML}^2)=\frac{1}{n}E\left|\sum_{k=1}^{n}(Z_k-\hat{\mu}_{ML})^2\right|=\frac{n-1}{n}\sigma^2 \tag{2.4.12}\] \(\hat{\mu}_{ML}\) 是 \(\mu\) 的无偏估计量，而 \(\hat{\sigma}_{ML}^2\) 是 \(\sigma^2\) 的有偏估计量，但当 \(n\rightarrow\infty\) 时，\(E(\hat{\sigma}_{ML}^2)\rightarrow\sigma^2\)，因此，\(\hat{\sigma}_{ML}^2\) 是 \(\sigma^2\) 的渐近无偏估计量。

从结果看，极大似然估计结果并不总是线性估计，而且也不总是无偏估计。

例 2.10设有观测方程 \(\bm{Z}=\bm{H}\bm{X}+\bm{\Delta}\)，\(\bm{\Delta}\) 为正态随机向量，\(\bm{\Delta}\sim N(\bm{0},\ \bm{D}_{\Delta})\)，\(\bm{X}\) 为常量，求参数 \(\bm{X}\) 的最大似然估计 \(\hat{\bm{X}}_{ML}\)。

解：由于 \(\bm{X}\) 为常量，有 \(\bm{Z}\sim N(\bm{H}\bm{X},\ \bm{D}_{\Delta})\)，观测值向量的概率密度函数为 \[p_{Z\mid X}(\bm{z}\mid\bm{x})=\frac{1}{(2\pi)^{\ell/2}\left|\bm{D}_{\Delta}\right|^{1/2}} \exp\left|-\frac{1}{2}(\bm{z}-\bm{H}\bm{x})^{\mathrm{T}}\bm{D}_{\Delta}^{-1}(\bm{z}-\bm{H}\bm{x})\right| \tag{2.4.13}\] 当上式中的指数部分最小时，\(p_{Z\mid X}(\bm{z}\mid\bm{x})\) 可以得到最大值，即 \[(\bm{z}-\bm{H}\bm{x})^{\mathrm{T}}\bm{D}_{\Delta}^{-1}(\bm{z}-\bm{H}\bm{x})=\min \tag{2.4.14}\] 将上式对 \(\bm{x}\) 的导数，并令导数等于 \(0\) \[2\bm{H}^{\mathrm{T}}\bm{D}_{\Delta}^{-1}(\bm{z}-\bm{H}\bm{x})=0 \tag{2.4.15}\] 设极大似然解为 \(\hat{\bm{x}}_{ML}\)，那么 \[\hat{\bm{x}}_{ML}=\left(\bm{H}^{\mathrm{T}}\bm{D}_{\Delta}^{-1}\bm{H}\right)^{-1}\bm{H}^{\mathrm{T}}\bm{D}_{\Delta}^{-1}\bm{z} \tag{2.4.16}\] 用样本代替以上变量 \[\hat{\bm{X}}_{ML}=\left(\bm{H}^{\mathrm{T}}\bm{D}_{\Delta}^{-1}\bm{H}\right)^{-1}\bm{H}^{\mathrm{T}}\bm{D}_{\Delta}^{-1}\bm{Z} \tag{2.4.17}\] 上例说明，对于线性模型 \(\bm{Z}=\bm{H}\bm{X}+\bm{\Delta}\)，且观测值误差 \(\bm{\Delta}\) 在正态分布的情况下，极大似然估计与最小二乘估计等价。

例 2.11设有观测方程 \(\bm{Z}=\bm{H}\bm{X}+\bm{\Delta}\)，随机误差 \(\bm{\Delta}\) 服从拉普拉斯分布 \[p(\bm{\Delta})=(1/2a)e^{-\|\bm{\Delta}\|_1/a} \tag{2.4.18}\] 上式中 \(a(a>0)\) 为常数；\(\|\cdot\|_1\) 表示向量的 1-范数。求参数 \(\bm{X}\) 的极大似然估计。

解：已知的 \(\bm{\Delta}\) 分布，可求得 \(\bm{Z}\) 的分布为

\[p(\bm{Z}\mid\bm{X})=(1/2a)e^{-\|\bm{Z}-\bm{H}\bm{X}\|_1/a} \tag{2.4.19}\] 将上式取对数 \[\ln p(\bm{Z}\mid\bm{X})=\ln(1/2a)-\|\bm{Z}-\bm{H}\bm{X}\|_1/a \tag{2.4.20}\] 上式的 \(\|\bm{Z}-\bm{H}\bm{X}\|_1\) 最小的时候，有最大的 \(\ln p(\bm{Z}\mid\bm{X})\)。所以 \[\hat{\bm{X}}_{ML}=\arg\,\min_{\hat{\bm{x}}}\|\bm{Z}-\bm{H}\bm{X}\|_1 \tag{2.4.21}\] 从结果来看，\(\hat{\bm{X}}_{ML}\) 是使 \(\bm{Z}-\bm{H}\bm{X}\) 的 1-范数最小的解，由于无法给出 \(\hat{\bm{X}}_{ML}\) 的解析解，只能通过数值方法求得。

将极大似然估计与最小二乘估计比较分析可以看到：

(1) 极大似然估计和最小二乘估计都不考虑参数的先验分布，也就是说，在极大似然估计和最小二乘估计中，都将 \(\bm{X}\) 视为非随机变量；

(2) 极大似然估计需要已知观测值的概率密度函数，即通过已知观测值的分布来建立似然函数，而最小二乘估计只需要知道观测值与参数的函数关系和观测值的特征值。

(3) 极大似然估计不总是无偏估计，也并不是线性估计，有时甚至不能得到解析解，但在例 2.10 中的特殊情况下：观测值误差服从正态分布，并且观测值与参数之间是线性关系时，极大似然估计与最小二乘估计等价。

极大验后估计

极大似然估计是以“\(p(\bm{z}\mid\bm{x})=\max\)”为准则的估计，在其估计中并没有考虑参数 \(\bm{X}\) 的先验信息。如果我们事先从经验和历史资料中已经知道了参数 \(\bm{X}\) 的信息，就希望在对参数进行估计的时候能够利用这些已知的先验信息，以得到对参数更为准确的估计和判断。例如，在递推的最小二乘估计中，从递推式上看，在对参数的第 \(k+1\) 次估计中，就利用了第 \(k\) 次估计的信息。如果我们已知 \(\bm{X}\) 的先验随机信息，在获得观测值后，应该利用观测和先验信息一并对 \(\bm{X}\) 进行估计，这也是贝叶斯学派的观点。

极大验后估计

设 \(p(\bm{x})\) 是 \(\bm{x}\) 的先验概率密度函数，\(p_z(\bm{z})\) 是观测值 \(\bm{Z}\) 的概率密度函数，\(p(\bm{z}\mid\bm{x})\) 是观测值的条件概率密度函数，那么有 \[p(\bm{x}\mid\bm{z})=\frac{p(\bm{z}\mid\bm{x})p_x(\bm{x})}{p_z(\bm{z})} \tag{2.5.1}\] \(p(\bm{x}\mid\bm{z})\) 是验后条件概率密度函数。极大验后估计就是使 \[p(\bm{x}\mid\bm{z})=\max \tag{2.5.2}\] 它的含义是：给定了观测值 \(\bm{Z}=\left[\begin{array}{llll}Z_1, & Z_2, & \cdots, & Z_{\ell}\end{array}\right]\) 的条件下使得验后分布的 \(\bm{x}\) 有最大的概率。记极大验后估计为 \(\hat{\bm{X}}_{MAP}\) \[\hat{\bm{X}}_{MAP}=\arg\,\max_{\hat{\bm{x}}}p(\hat{\bm{x}}\mid\bm{z}) \tag{2.5.3}\] 由于式 (2.5.2) 是求 \(\bm{x}\) 使得 \(p(\bm{x}\mid\bm{z})\) 最大，与 \(p(\bm{z})\) 没有关系，所以 \[p(\bm{x}\mid\bm{z})\propto p(\bm{z}\mid\bm{x})p_x(\bm{x}) \tag{2.5.4}\] 上式的 \(\propto\) 表示“正比例于”，因此，极大验后估计也为 \[\hat{\bm{X}}_{MAP}=\arg\,\max_{\hat{\bm{x}}}p(\bm{z}\mid\hat{\bm{x}})p_x(\hat{\bm{x}}) \tag{2.5.5}\] 式 (2.5.4) 和 (2.5.5) 表明：求解最大验后概率相当于最大化似然概率 \(p(\bm{z}\mid\bm{x})\) 和先验概率 \(p_x(\bm{x})\) 的乘积。当直接求解后验概率分布 \(p(\bm{x}\mid\bm{z})\) 困难时，就可以通过式 (2.5.5) 来求解 \(\hat{\bm{X}}_{MAP}\)。

在很多时候，\(p(\bm{x}\mid\bm{z})\) 和 \(p(\bm{z}\mid\bm{x})p_x(\bm{x})\) 中有指数函数，所以先对其取自然对数后再求极大值更加方便，即 \[\left.\frac{\partial\ln\left[\,p(\bm{x}\mid\bm{z})\,\right]}{\partial\bm{x}}\right|_{\bm{x}=\hat{\bm{x}}_{MAP}}=0 \tag{2.5.6}\] 或者 \[\left.\left[\frac{\partial\ln\left[\,p(\bm{z}\mid\bm{x})\,\right]}{\partial\bm{x}}+\frac{\partial\ln\left[\,p_x(\bm{x})\,\right]}{\partial\bm{x}}\right]\right|_{\bm{x}=\hat{\bm{x}}_{MAP}}=0 \tag{2.5.7}\]

算例分析

例 2.12某随机变量 \(z\) 的概率密度函数为 \[p_{z/\theta}(z/\theta)=\begin{cases}\theta e^{-\theta z}, & z\geq 0\\ 0, & \text{其他}\end{cases}\] 其中参数 \(\theta\) 为随机量，概率密度函数为 \[p_\theta(\theta)=\begin{cases}\dfrac{1}{\theta}, & 0<\theta\leq 1\\[6pt] 0, & \text{其他}\end{cases}\] 现在对随机变量 \(z\) 进行独立观测，\(Z_1\ \cdots\ Z_n\) 为观测值，且观测值 \(Z_1\ \cdots\ Z_n\) 条件随机独立。求参数 \(\theta\) 的极大验后估计。

解：由于各观测值相互随机独立，所以 \(\bm{Z}=\left[\begin{array}{lll}Z_1 & \cdots & Z_{\ell}\end{array}\right]^{\mathrm{T}}\) 的概率密度函数为 \[p(z_1,\ z_2\cdots z_{\ell}/\theta)=\begin{cases}\theta^{\ell}e^{-\theta\sum\limits_1^{\ell}z_i}, & z\geq 0\\[6pt] 0, & \text{其他}\end{cases} \tag{2.5.8}\] 设 \(p(\bm{z},\ \theta)\) 为 \(\bm{Z}\) 和 \(\theta\) 的联合概率密度函数 \[p(\bm{z},\ \theta)=p(z_1,\ z_2\cdots z_{\ell}/\theta)\,p_\theta(\theta) =\begin{cases}\theta^{\ell-1}e^{-\theta\sum\limits_1^{\ell}z_i}, & 0<\theta\leq 1,\ z\geq 0\\[6pt] 0, & \text{其他}\end{cases} \tag{2.5.9}\] \(p(\bm{z},\ \theta)=\max\) 等价于 \(\ln\left(\,p(\bm{z},\ \theta)\,\right)=\max\) \[\ln\left(\,P(\bm{z},\ \theta)\,\right)=(\ell-1)\ln\theta-\theta\sum_1^{\ell}z_i \tag{2.5.10}\] 为了得到极大值，将 \(\ln\left(\,p(\bm{z},\ \theta)\,\right)\) 对 \(\theta\) 求导，并等于零，得到 \[\hat{\theta}=\frac{\ell-1}{\displaystyle\sum_1^{\ell}z_i} \tag{2.5.11}\] 用样本 \(Z_i\) 代替上式中的 \(z_i\)，\(\theta\) 的极大验后估计为 \[\hat{\theta}=\frac{\ell-1}{\displaystyle\sum_1^{\ell}Z_i}\quad (0<\theta\leq 1,\ z\geq 0)\]

例 2.13设有观测方程 \(\bm{Z}=\bm{H}\bm{X}+\bm{\Delta}\)，且 \(\bm{\Delta}\sim N_{\ell}(\bm{0},\ \bm{D}_{\Delta})\) 和 \(\bm{X}\sim N_n(\bm{\mu}_x,\ \bm{D}_X)\)，\(\bm{\Delta}\) 与 \(\bm{X}\) 相互独立，求参数 \(\bm{X}\) 的极大验后估计和它的方差。

解：由于 \(\bm{\Delta}\sim N_{\ell}(\bm{0},\ \bm{D}_{\Delta})\) 和 \(\bm{X}\sim N_n(\bm{\mu}_x,\ \bm{D}_X)\)，\(\bm{\Delta}\) 与 \(\bm{X}\) 相互，所以 \(\bm{Z}\) 也服从正态分布分布，即 \(\bm{Z}\sim N_{\ell}(\bm{\mu}_Z,\ \bm{D}_Z)\)，其中 \[\bm{\mu}_Z=E(\bm{Z})=\bm{H}\bm{\mu}_x \tag{2.5.12}\] \[\bm{D}_Z=\bm{H}\bm{D}_X\bm{H}^{\mathrm{T}}+\bm{D}_{\Delta} \tag{2.5.13}\] 并且容易求得 \[\bm{D}_{XZ}=\bm{D}_X\bm{H}^{\mathrm{T}} \tag{2.5.14}\] 根据 1.4.6 节可知，已知 \(\bm{X}\sim N_n(\bm{\mu}_x,\ \bm{D}_X)\) 和 \(\bm{Z}\sim N_{\ell}(\bm{\mu}_Z,\ \bm{D}_Z)\)，那么以 \(\bm{Z}\) 为条件，\(\bm{X}\) 的概率密度函数为 \[p(\bm{X}\mid\bm{Z})=(2\pi)^{-\frac{n}{2}}\left|\bm{D}_{x/z}\right|^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}\left(\bm{X}-\bm{\mu}_{x/z}\right)^{\mathrm{T}}\bm{D}_{x/z}^{-1}\left(\bm{X}-\bm{\mu}_{x/z}\right)\right\} \tag{2.5.15}\] 当 \(p(\bm{X}\mid\bm{Z})\) 中的指数部分 \(\left(\bm{X}-\bm{\mu}_{x/z}\right)^{\mathrm{T}}\bm{D}_{x/z}^{-1}\left(\bm{X}-\bm{\mu}_{x/z}\right)\) 最小时，\(p(\bm{X}\mid\bm{Z})\) 有最大值，所以极大验后估计为 \[\hat{\bm{X}}_{MAP}=\bm{\mu}_{x/z} \tag{2.5.16}\] \(\hat{\bm{X}}_{MAP}\) 的估计误差为 \[\begin{aligned} \Delta\hat{\bm{X}}_{MAP}&=\bm{X}-\hat{\bm{X}}_{MAP}\\ &=\bm{X}-\bm{\mu}_{x/z} \end{aligned} \tag{2.5.17}\] 那么 \(\hat{\bm{X}}_{MAP}\) 的方差为 \[\begin{aligned} \mathrm{Var}(\hat{\bm{X}}_{MAP})&=E\left[\,\Delta\hat{\bm{X}}_{MAP}\Delta\hat{\bm{X}}_{MAP}^{\mathrm{T}}\,\right]\\ &=E\left[\,\left(\bm{X}-\bm{\mu}_{x/z}\right)\left(\bm{X}-\bm{\mu}_{x/z}\right)^{\mathrm{T}}\,\right] \end{aligned} \tag{2.5.18}\] 注意到 \(E\left[\,\left(\bm{X}-\bm{\mu}_{x/z}\right)\left(\bm{X}-\bm{\mu}_{x/z}\right)^{\mathrm{T}}\,\right]\) 为 \(\bm{X}\) 的验后方差的定义 \[\bm{D}_{x/z}=E\left[\,\left(\bm{X}-\bm{\mu}_{x/z}\right)\left(\bm{X}-\bm{\mu}_{x/z}\right)^{\mathrm{T}}\,\right] \tag{2.5.19}\] 所以 \[\mathrm{Var}(\hat{\bm{X}}_{MAP})=\bm{D}_{x/z} \tag{2.5.20}\] 式 (2.5.16) 和式 (2.5.20) 说明，当 \(\bm{X}\) 和 \(\bm{Z}\) 都是正态分布时，\(\bm{X}\) 的极大验后估计即为验后期望 \(\bm{\mu}_{x/z}\)，其方差为验后方差 \(\bm{D}_{x/z}\)。

根据 1.4.6 节可知 \[\begin{aligned} \hat{\bm{X}}_{MAP}&=\bm{\mu}_{x/z}\\ &=\bm{\mu}_x+\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-\bm{\mu}_Z\right) \end{aligned} \tag{2.5.21}\] \[\begin{aligned} \mathrm{Var}(\hat{\bm{X}}_{MAP})&=\bm{D}_{x/z}\\ &=\bm{D}_X-\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_{ZX} \end{aligned} \tag{2.5.22}\] 将式 (2.5.12)、式 (2.5.13) 和式 (2.5.14) 代入式 (2.5.21) 和式 (2.5.22)，得到 \[\hat{\bm{X}}_{MAP}=\bm{\mu}_x+\bm{D}_X\bm{H}^{\mathrm{T}}\left(\bm{H}\bm{D}_X\bm{H}^{\mathrm{T}}+\bm{D}_{\Delta}\right)^{-1}\left(\bm{Z}-\bm{H}\bm{\mu}_x\right) \tag{2.5.23}\] \[\mathrm{Var}(\hat{\bm{X}}_{MAP})=\bm{D}_X-\bm{D}_X\bm{H}^{\mathrm{T}}\left(\bm{H}\bm{D}_X\bm{H}^{\mathrm{T}}+\bm{D}_{\Delta}\right)^{-1}\bm{H}\bm{D}_X \tag{2.5.24}\]

若在 \(t_k\) 时刻对 \(\bm{X}\) 的估计为 \(\hat{\bm{X}}_k\) 和 \(\bm{D}_{\hat{X}_k}\)，在 \(t_{k+1}\) 时刻有观测方程 \[\bm{z}_{k+1}=\bm{h}_{k+1}\bm{X}+\bm{\Delta}_{k+1} \tag{2.5.25}\] 将在 \(t_k\) 时刻的估计视为先验信息，并设 \(\bm{X}\sim N(\hat{\bm{X}}_k,\ \bm{D}_{\hat{X}_k})\)，根据式 (2.5.23) 和 (2.5.24)，可得到 \(t_{k+1}\) 时刻的极大验后估计为 \[\hat{\bm{X}}_{k+1}=\hat{\bm{X}}_k+\bm{D}_{\hat{X}_k}\bm{h}_{k+1}^{\mathrm{T}} \left(\bm{h}_{k+1}\bm{D}_{\hat{X}_k}\bm{h}_{k+1}^{\mathrm{T}}+\bm{D}_{\Delta_k}\right)^{-1} \left(\bm{z}_{k+1}-\bm{h}_{k+1}\hat{\bm{X}}_k\right) \tag{2.5.26}\] \[\bm{D}_{\hat{X}_{k+1}}=\bm{D}_{\hat{X}_k}-\bm{D}_{\hat{X}_k}\bm{h}_{k+1}^{\mathrm{T}} \left(\bm{h}_{k+1}\bm{D}_{\hat{X}_k}\bm{h}_{k+1}^{\mathrm{T}}+\bm{D}_{\Delta_k}\right)^{-1}\bm{h}_{k+1}\bm{D}_{\hat{X}_k} \tag{2.5.27}\] 观察式 (2.5.26) 和式 (2.5.27) 可知它们即为递推的最小二乘公式。

比较已经介绍的几种估计方法，可以总结得到：

(1) 极大验后估计将待估计参数 \(\bm{X}\) 视为随机量，在估计时考虑参数的先验随机信息，最大化似然概率 \(p(\bm{z}\mid\hat{\bm{x}})\) 和先验概率 \(p(\hat{\bm{x}})\) 的乘积：\(\hat{\bm{X}}_{MAP}=\arg\max\limits_{\hat{x}}p(\bm{z}\mid\bm{x})p_x(\bm{x})\)。如果先验概率 \(p_x(\bm{x})\) 未知，或者简单地认为先验概率为均匀分布，那么极大验后估计就退化为极大似然估计了。

(2) 极大验后估计和极大似然估计一样，不总是观测值的线性函数；也没有如最小二乘估计那样固定的解析表达式。

(3) 当 \(\bm{X}\) 和 \(\bm{Z}\) 都是正态分布时，\(\bm{X}\) 极大验后估计为验后期望 \(\bm{\mu}_{x\mid z}\)，其方差为验后方差 \(\bm{D}_{x/z}\)。

最小方差估计

最小方差估计

准确地讲，最小方差估计指的是均方差最小估计。均方差指的是估计值与其真值之间的密集程度或者估计值的真误差在零附近的密集程度，它是评价估计质量的重要指标。由于这里的最小方差估计得到的是无偏估计，均方差就是其方差，所以也被称为最小方差估计。

1. 最小方差估计

最小方差估计的准则为 \[L_0(\hat{\bm{X}}(\bm{Z}))=E\left[\,(\hat{\bm{X}}-\bm{X})(\hat{\bm{X}}-\bm{X})^{\mathrm{T}}\,\right]=\min \tag{2.6.1}\] 上式中 \(\hat{\bm{X}}\) 是通过观测值 \(\bm{Z}\) 得到的对 \(\bm{X}\) 的估计，它是 \(\bm{Z}\) 的函数，所以 \((\hat{\bm{X}}-\bm{X})(\hat{\bm{X}}-\bm{X})^{\mathrm{T}}\) 是 \(\bm{X}\) 和 \(\bm{Z}\) 的函数，假设 \(p(\bm{x},\ \bm{z})\) 为 \(\bm{X}\) 和 \(\bm{Z}\) 的联合分布，那么 \[\begin{aligned} L_0(\hat{\bm{X}})&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})(\hat{\bm{x}}-\bm{x})^{\mathrm{T}} p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{x}\mathrm{d}\bm{z}\\ &=\int_{-\infty}^{\infty}\left\{\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})(\hat{\bm{x}}-\bm{x})^{\mathrm{T}} p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}p_z(\bm{z})\,\mathrm{d}\bm{z} \end{aligned} \tag{2.6.2}\] 设上式花括号部分为 \[L(\hat{\bm{X}})=\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})(\hat{\bm{x}}-\bm{x})^{\mathrm{T}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\] 那么 \[L_0(\hat{\bm{X}})=\int_{-\infty}^{\infty}L(\hat{\bm{X}})p_z(\bm{z})\,\mathrm{d}\bm{z}\] 上式中的 \(p(\bm{z})\) 是非负的，所以 \(L_0(\hat{\bm{X}}(\bm{Z}))=\min\) 等价于 \[L(\hat{\bm{X}})=\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})(\hat{\bm{x}}-\bm{x})^{\mathrm{T}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=\min \tag{2.6.3}\] 将 \(L(\hat{\bm{X}})\) 展开为 \[\begin{aligned} L(\hat{\bm{X}})=&\int_{-\infty}^{\infty}\left\{(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z})+E(\bm{X}\mid\bm{Z})-\bm{x}) (\hat{\bm{x}}-E(\bm{X}\mid\bm{Z})+E(\bm{X}\mid\bm{Z})-\bm{x})^{\mathrm{T}}\right\}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ =&\int_{-\infty}^{\infty}(E(\bm{X}\mid\bm{Z})-\bm{x})(E(\bm{X}\mid\bm{Z})-\bm{x})^{\mathrm{T}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ &+(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))^{\mathrm{T}}\int_{-\infty}^{\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ &+\left\{\int_{-\infty}^{\infty}(E(\bm{X}\mid\bm{Z})-\bm{x})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))^{\mathrm{T}}\\ &+(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))\int_{-\infty}^{\infty}(\bm{x}-E(\bm{X}\mid\bm{Z}))^{\mathrm{T}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \end{aligned} \tag{2.6.4}\] 上式中的 \((\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))\) 之所以可以放在积分外，是因为 \(\hat{\bm{x}}\) 和 \(E(\bm{X}\mid\bm{Z})\) 都是 \(\bm{z}\) 的函数，与 \(\bm{x}\) 无关；上式中的第三项为 \[\begin{aligned} &\left\{\int_{-\infty}^{\infty}(E(\bm{X}\mid\bm{Z})-\bm{x})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))^{\mathrm{T}}\\ =&\left\{\int_{-\infty}^{\infty}E(\bm{X}\mid\bm{Z})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} -\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}(E(\bm{X}\mid\bm{Z})-\hat{\bm{x}})^{\mathrm{T}}\\ =&\left\{E(\bm{X}\mid\bm{Z})-E(\bm{X}\mid\bm{Z})\right\}(E(\bm{X}\mid\bm{Z})-\hat{\bm{x}})^{\mathrm{T}}\\ =&0 \end{aligned} \tag{2.6.5}\]

同理，第四项也为零，所以 \[\begin{aligned} L(\hat{\bm{X}})=&\int_{-\infty}^{\infty}(E(\bm{X}\mid\bm{Z})-\bm{x})(E(\bm{X}\mid\bm{Z})-\bm{x})^{\mathrm{T}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ &+(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))(\hat{\bm{x}}-E(\bm{X}\mid\bm{Z}))^{\mathrm{T}}\int_{-\infty}^{\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \end{aligned} \tag{2.6.6}\] 上式中的第一项为非负矩阵，第二项中的 \(\int_{-\infty}^{\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=1\)，所以使 \(L(\hat{\bm{X}})\) 最小就是使第二项为零矩阵，那么 \[E(\bm{X}\mid\bm{Z})-\hat{\bm{x}}=0 \tag{2.6.7}\] 若记在最小方差准则下的参数估计为 \(\hat{\bm{X}}_{MV}\)，那么 \[\hat{\bm{X}}_{MV}=E(\bm{X}\mid\bm{Z}) \tag{2.6.8}\]

由于方差矩阵迹最小的估计与方差矩阵最小的估计完全相同，所以也可以用方差矩阵迹最小的准则来代替方差最小准则得到 \(\hat{\bm{X}}_{MV}\)，下面给出证明过程： \[\mathrm{tr}\left(L_0(\hat{\bm{X}}(\bm{Z}))\right)=\min \tag{2.6.9}\] 那么， \[\mathrm{tr}\left(L_0(\hat{\bm{X}}(\bm{Z}))\right)=E\left[\,(\hat{\bm{X}}-\bm{X})^{\mathrm{T}}(\hat{\bm{X}}-\bm{X})\,\right]=\min \tag{2.6.10}\] 根据期望的定义 \[\begin{aligned} \mathrm{tr}\left(L_0(\hat{\bm{X}}(\bm{Z}))\right) &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})^{\mathrm{T}}(\hat{\bm{x}}-\bm{x}) p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{x}\mathrm{d}\bm{z}\\ &=\int_{-\infty}^{\infty}\left\{\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})^{\mathrm{T}}(\hat{\bm{x}}-\bm{x}) p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}p_z(\bm{z})\,\mathrm{d}\bm{z} \end{aligned} \tag{2.6.11}\] 这等价于 \[\int_{-\infty}^{\infty}(\hat{\bm{x}}-\bm{x})^{\mathrm{T}}(\hat{\bm{x}}-\bm{x})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=\min \tag{2.6.12}\] 展开上式为 \[\hat{\bm{x}}^{\mathrm{T}}\hat{\bm{x}}+\int_{-\infty}^{\infty}\bm{x}^{\mathrm{T}}\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} -2\hat{\bm{x}}^{\mathrm{T}}\int_{-\infty}^{\infty}\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=\min \tag{2.6.13}\] 将上式对 \(\hat{\bm{x}}\) 求导并令其为零得到 \[2\hat{\bm{x}}-2\int_{-\infty}^{\infty}\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=0 \tag{2.6.14}\] 注意到 \(\int_{-\infty}^{\infty}\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}=E(\bm{X}\mid\bm{Z})\)，所以从方差矩阵迹最小也证明了最小方差估计为 \[\hat{\bm{X}}_{MV}=E(\bm{X}\mid\bm{Z}) \tag{2.6.15}\]

2. 最小方差估计的统计特性

式 (2.6.15) 表明，最小方差估计是在 \(\bm{Z}\) 的条件下 \(\bm{X}\) 的条件期望，它是 \(\bm{Z}\) 的函数，所以 \(\hat{\bm{X}}_{MV}\) 的期望为 \[\begin{aligned} E(\hat{\bm{X}}_{MV})&=\int_{-\infty}^{\infty}E(\bm{x}\mid\bm{z})p_z(\bm{z})\,\mathrm{d}\bm{z}\\ &=\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}\bm{x}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right]p_z(\bm{z})\,\mathrm{d}\bm{z}\\ &=\int_{-\infty}^{\infty}\bm{x}\left[\int_{-\infty}^{\infty}p(\bm{x}\mid\bm{z})p_z(\bm{z})\,\mathrm{d}\bm{z}\right]\mathrm{d}\bm{x}\\ &=\int_{-\infty}^{\infty}\bm{x}\left[\int_{-\infty}^{\infty}p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{z}\right]\mathrm{d}\bm{x}\\ &=\int_{-\infty}^{\infty}\bm{x}p_x(\bm{x})\,\mathrm{d}\bm{x}\\ &=E(\bm{X}) \end{aligned} \tag{2.6.16}\] 显然，最小方差估计是无偏估计。\(\hat{\bm{X}}_{MV}\) 的方差为 \[D(\hat{\bm{X}}_{MV})=E\left\{(\hat{\bm{X}}_{MV}-\bm{X})(\hat{\bm{X}}_{MV}-\bm{X})^{\mathrm{T}}\right\} \tag{2.6.17}\] 将式 (2.6.15) 代入上式 \[D(\hat{\bm{X}}_{MV})=E\left\{(E(\bm{X}\mid\bm{Z})-\bm{X}))(E(\bm{X}\mid\bm{Z})-\bm{X})^{\mathrm{T}}\right\} \tag{2.6.18}\] \(E\left\{(E(\bm{X}\mid\bm{Z})-\bm{X}))(E(\bm{X}\mid\bm{Z})-\bm{X})^{\mathrm{T}}\right\}\) 即为 \(\bm{X}\) 的验后方差的定义，所以 \[D(\hat{\bm{X}}_{MV})=D(\bm{X}\mid\bm{Z}) \tag{2.6.19}\] 从上面的分析看出，无论何种分布，最小方差估计 \(\hat{\bm{X}}_{MV}\) 都为其分布的验后期望 \(E(\bm{X}\mid\bm{Z})\)，\(\hat{\bm{X}}_{MV}\) 的方差为其验后方差 \(D(\bm{X}\mid\bm{Z})\)。在上一节的分析中我们看到，在正态分布情况下，极大验后估计也为 \(E(\bm{X}\mid\bm{Z})\)，也就是说，在正态分布情况下，极大验后估计与最小方差估计等价，如果 \(\bm{X}\) 和 \(\bm{Z}\) 不是正态分布随机向量，极大似然估计就不一定是其验后期望，也就不与最小方差估计等价了。

线性最小方差估计

最小方差估计是其分布的验后期望，不同的分布函数有不同的验后期望，它可能是关于观测值的线性函数，也可能是非线性函数。如果估计是关于观测值的线性函数，并且满足均方差最小，这样的估计称为线性最小方差估计。

线性最小方差估计是一种特殊的最小方差估计，它是指估计值 \(\hat{\bm{X}}_{MV}\) 是观测量 \(\bm{Z}\) 的线性函数，并使得估计的均方差最小的估计，估计量具有如下形式 \[\underbrace{\hat{\bm{X}}}_{n\times 1}=\underbrace{\bm{a}_L}_{n\times 1}+\underbrace{\bm{B}_L}_{n\times\ell}\underbrace{\bm{Z}}_{\ell\times 1} \tag{2.6.20}\] 且 \(\hat{\bm{X}}\) 的均方误差最小 \[L_0(\hat{\bm{X}})=\mathrm{MSE}(\hat{\bm{X}})=E\left[\,(\bm{X}-\hat{\bm{X}})(\bm{X}-\hat{\bm{X}})^{\mathrm{T}}\,\right]=\min \tag{2.6.21}\] 从上一节分析知道 \(L_0(\hat{\bm{X}})=\min\) 等价于 \[\mathrm{tr}\left[L_0(\hat{\bm{X}})\right]=\mathrm{tr}\left\{E\left[\,(\bm{X}-\hat{\bm{X}})(\bm{X}-\hat{\bm{X}})^{\mathrm{T}}\,\right]\right\}=\min \tag{2.6.22}\] 将式 (2.6.20) 代入式 (2.6.22) 有 \[\mathrm{tr}\left(L_0(\hat{\bm{X}})\right)=E\left[\,(\bm{X}-\bm{a}_L-\bm{B}_L\bm{Z})^{\mathrm{T}}(\bm{X}-\bm{a}_L-\bm{B}_L\bm{Z})\,\right]=\min \tag{2.6.23}\] 根据极值理论将 \(\mathrm{tr}\left(L_0(\hat{\bm{X}})\right)\) 对 \(\bm{a}_L\) 和 \(\bm{B}_L\) 分别求偏导并令其为零 \[E(\bm{X}-\bm{a}_L-\bm{B}_L\bm{Z})=0 \tag{2.6.24}\] \[E((\bm{X}-\bm{a}_L-\bm{B}_L\bm{Z})\bm{Z}^{\mathrm{T}})=0 \tag{2.6.25}\] 将式 (2.6.24) 和式 (2.6.25) 联立，解得 \[\left. \begin{aligned} \bm{a}_L&=E(\bm{X})-\bm{D}_{XZ}\bm{D}_Z^{-1}E(\bm{Z})\\ \bm{B}_L&=\bm{D}_{XZ}\bm{D}_Z^{-1} \end{aligned} \quad\right\} \tag{2.6.26}\] 将 \(\bm{a}_L\) 和 \(\bm{B}_L\) 代入式 (2.6.20)，得到最小方差估计 \[\begin{aligned} \hat{\bm{X}}_L&=\bm{a}_L+\bm{B}_L\bm{Z}\\ &=E(\bm{X})+\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right) \end{aligned} \tag{2.6.27}\] 由于线性最小方差估计是最小方差估计的特殊形式，所以具有最小方差估计的无偏性。\(\hat{\bm{X}}_L\) 的估计误差为 \[\begin{aligned} \Delta\hat{\bm{X}}_L&=\bm{X}-\hat{\bm{X}}_L\\ &=\bm{X}-E(\bm{X})-\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right) \end{aligned} \tag{2.6.28}\] \(\hat{\bm{X}}_L\) 的方差为 \[\begin{aligned} \mathrm{Var}(\hat{\bm{X}}_L)=&E\left[\,\left(\hat{\bm{X}}_L-\bm{X}\right)\left(\hat{\bm{X}}_L-\bm{X}\right)^{\mathrm{T}}\,\right]\\ =&E\left\{\left[\bm{X}-E(\bm{X})-\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\right.\\ &\qquad\left.\cdot\left[\bm{X}-E(\bm{X})-\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\}\\ =&E\left\{\left[\bm{X}-E(\bm{X})\right]\left[\bm{X}-E(\bm{X})\right]^{\mathrm{T}}\right\} -E\left\{\left[\bm{X}-E(\bm{X})\right]\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\}\\ &-E\left\{\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\left[\bm{X}-E(\bm{X})\right]^{\mathrm{T}}\right\}\\ &+E\left\{\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\} \end{aligned} \tag{2.6.29}\] 上式最后一个等号后的第一项为 \(\bm{X}^{\mathrm{T}}\) 的方差

\[E\left\{\left[\bm{X}-E(\bm{X})\right]\left[\bm{X}-E(\bm{X})\right]^{\mathrm{T}}\right\}=\bm{D}_X \tag{2.6.30}\] 第二项为 \[\begin{aligned} &E\left\{\left[\bm{X}-E(\bm{X})\right]\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\}\\ =&E\left\{\left[\bm{X}-E(\bm{X})\right]\left[\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\}\bm{D}_Z^{-1}\bm{D}_{ZX}\\ =&\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_{ZX} \end{aligned} \tag{2.6.31}\] 同样，第三项也为 \[E\left\{\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\left[\bm{X}-E(\bm{X})\right]^{\mathrm{T}}\right\} =\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_{ZX} \tag{2.6.32}\] 第四项为 \[\begin{aligned} &E\left\{\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]^{\mathrm{T}}\right\}\\ =&E\left\{\left[\bm{D}_{XZ}\bm{D}_Z^{-1}\left(\bm{Z}-E(\bm{Z})\right)\right]\left[\left(\bm{Z}-E(\bm{Z})\right)^{\mathrm{T}}\bm{D}_Z^{-1}\bm{D}_{ZX}\right]\right\}\\ =&\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_Z\bm{D}_Z^{-1}\bm{D}_{ZX}\\ =&\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_{ZX} \end{aligned} \tag{2.6.33}\] 综合上述， \[\bm{D}_{\hat{\bm{X}}_L}=\bm{D}_X-\bm{D}_{XZ}\bm{D}_Z^{-1}\bm{D}_{ZX} \tag{2.6.34}\]

贝叶斯估计

贝叶斯估计和贝叶斯相关理论由英国神甫托马斯·贝叶斯（1702–1761 年）提出，它对统计学界和估计理论产生了深远的影响，广泛地应用于模式识别和人工智能等领域。本节将介绍贝叶斯理论中的贝叶斯估计方法。

在参数估计中，不考虑参数的先验信息，这是统计学中频率学派的观点，他们认为参数虽然未知，但参数是固定的常数。在已经学习的估计方法中，最小二乘估计和极大似然估计都不考虑参数的先验信息，最小二乘估计甚至不需要任何随机变量的分布，所以在工程中有最广泛的应用。贝叶斯学派认为参数不是常数，它是变化的随机量，它的变化可以用一个概率分布来描述，并且在参数估计时应该利用参数的先验信息。前面介绍的极大验后、最小方差估计都利用了参数的先验随机信息或分布，它们都属于贝叶斯估计。

在估计某个量时，随机误差的干扰使估计产生误差 \[\Delta\bm{X}=\bm{X}-\hat{\bm{X}}(\bm{Z}) \tag{2.7.1}\] 这种差异造成估计的“损失”，我们可以定义损失函数对其进行量化 \[L(\Delta\bm{X})=L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right) \tag{2.7.2}\] 根据需要可以定义不同的损失函数。一般而言，估计误差越大，损失就越大，典型的损失函数有平方损失函数、绝对值损失函数和均值损失函数，这三种损失函数如图 2.11 所示。

(1) 平方损失函数： \[L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)=\left(\bm{X}-\hat{\bm{X}}(\bm{Z})\right)\left(\bm{X}-\hat{\bm{X}}(\bm{Z})\right)^{\mathrm{T}} \tag{2.7.3}\]

(2) 绝对值损失函数： \[L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)=\left|\bm{X}-\hat{\bm{X}}(\bm{Z})\right| \tag{2.7.4}\]

(3) 均值损失函数： \[L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)=\begin{cases}0, & |\Delta\hat{\bm{X}}|\leq\Delta/2\\ 1, & |\Delta\hat{\bm{X}}|>\Delta/2\end{cases} \tag{2.7.5}\]

损失函数 \(L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)\) 是 \(\bm{X}\) 和 \(\bm{Z}\) 的函数，所以损失函数的期望为 \[E\left[\,L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)\,\right] =\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{x}\mathrm{d}\bm{z} \tag{2.7.6}\] 损失函数的期望 \(E\left[\,L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)\,\right]\) 即为贝叶斯风险（Bayes Risk）， \[R_B\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}L\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{x}\mathrm{d}\bm{z} \tag{2.7.7}\] 它表示损失函数的平均值。在决策论中，使贝叶斯风险最小的决策是最优决策。贝叶斯估计就是得到使贝叶斯风险值最小的估计 \(\hat{\bm{X}}_B(\bm{Z})\)： \[\hat{\bm{X}}_B(\bm{Z})=\mathrm{argmin}_{\hat{\bm{X}}}\left(\,R_B\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)\,\right) \tag{2.7.8}\] \(\hat{\bm{X}}_B(\bm{Z})\) 表示参数估计是 \(\bm{Z}\) 的函数，在下面的推导中，为了简单起见，记 \(\hat{\bm{X}}_B(\bm{Z})\) 为 \(\hat{\bm{X}}_B\)。

式 (2.7.7) 可以表示为 \[R_B(\bm{X},\ \hat{\bm{X}})=\int_{-\infty}^{+\infty}\left\{\int_{-\infty}^{+\infty}L(\bm{X},\ \hat{\bm{X}})p_{X/Z}(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\right\}p_Z(\bm{z})\,\mathrm{d}\bm{z} \tag{2.7.9}\] 记上式中的花括号部分为 \[r(\hat{\bm{X}}\mid\bm{Z})=\int_{-\infty}^{+\infty}L(\bm{X},\ \hat{\bm{X}})\,p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \tag{2.7.10}\] 所以 \[R_B(\bm{X},\ \hat{\bm{X}})=\int_{-\infty}^{+\infty}r(\hat{\bm{X}}\mid\bm{Z})\,p_z(\bm{z})\,\mathrm{d}\bm{z} \tag{2.7.11}\] 上式中的 \(r(\hat{\bm{X}}\bm{z})\) 是损失函数的验后条件分布的期望，所以也称为“验后风险”或者“验后期望损失”。由于 \(p_Z(\bm{z})\) 非负，并且与 \(\hat{\bm{X}}\) 无关，所以，贝叶斯风险最小等价于 \[r(\hat{\bm{X}}\mid\bm{Z})=\min \tag{2.7.12}\] 这说明验后风险最小估计和贝叶斯风险最小估计是等价的 \[\hat{\bm{X}}_B(\bm{Z})=\mathrm{argmin}_{\hat{\bm{X}}}\left(\,R_B\left(\bm{X},\ \hat{\bm{X}}(\bm{Z})\right)\,\right) =\mathrm{argmin}_{\hat{\bm{X}}}\left(\,r(\hat{\bm{X}}\mid\bm{Z})\,\right) \tag{2.7.13}\]

下面根据不同的损失函数来推导贝叶斯估计。

(1) 平方损失函数的贝叶斯估计：

平方损失函数的贝叶斯风险 \[R_B(\bm{X},\ \hat{\bm{X}})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\left[\,\bm{x}-\hat{\bm{x}}\,\right]\left[\,\bm{x}-\hat{\bm{x}}\,\right]^{\mathrm{T}}p(\bm{x},\ \bm{z})\,\mathrm{d}\bm{x}\mathrm{d}\bm{z} \tag{2.7.14}\] 将上式与式 (2.6.2) 比较发现，平方损失函数的贝叶斯风险就是其均方差，所以当损失函数为平方损失函数时，贝叶斯风险最小的参数估计 \(\hat{\bm{X}}_B\) 就是方差最小估计，有 \[\hat{\bm{X}}_B=\hat{\bm{X}}_{MV}=E(\bm{X}\mid\bm{Z}) \tag{2.7.15}\]

(2) 绝对值损失函数的贝叶斯估计：

绝对值损失函数的验后风险为 \[\begin{aligned} r_{abs}(\hat{\bm{X}}\mid\bm{Z})&=\int_{-\infty}^{+\infty}\left|\bm{x}-\hat{\bm{x}}\right|p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ &=\int_{-\infty}^{\hat{\bm{x}}}(\hat{\bm{x}}-\bm{x})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} +\int_{\hat{\bm{x}}}^{+\infty}(\bm{x}-\hat{\bm{x}})p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \end{aligned} \tag{2.7.16}\] 将上式对 \(\hat{\bm{x}}\) 求导 \[\begin{aligned} \frac{\mathrm{d}r_{abs}(\hat{\bm{X}}/\bm{Z})}{\mathrm{d}\hat{\bm{X}}} &=\int_{-\infty}^{\hat{x}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} +\hat{\bm{x}}p(\hat{\bm{x}}\mid\bm{z})-\hat{\bm{x}}p(\hat{\bm{x}}\mid\bm{z})\\ &\quad-\hat{\bm{x}}p(\hat{\bm{x}}\mid\bm{z})-\int_{\hat{x}}^{\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} +\hat{\bm{x}}p(\hat{\bm{x}}\mid\bm{z})\\ &=\int_{-\infty}^{\hat{x}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} -\int_{\hat{x}}^{+\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \end{aligned} \tag{2.7.17}\] 使 \(\dfrac{\mathrm{d}r_{abs}(\hat{\bm{X}}\mid\bm{Z})}{\mathrm{d}\hat{\bm{X}}}=0\)，得到 \[\int_{-\infty}^{\hat{x}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} =\int_{\hat{x}}^{+\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \tag{2.7.18}\] 记绝对值损失函数的贝叶斯估计为 \(\hat{\bm{X}}_{abs}\)，那么 \[\int_{-\infty}^{\hat{X}_{abs}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} =\int_{\hat{X}_{abs}}^{+\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \tag{2.7.19}\] \(\hat{\bm{X}}_{abs}\) 两侧的积分相等，也为积分中数 \(\hat{\bm{X}}_{med}\)，即 \(\hat{\bm{X}}_{abs}=\hat{\bm{X}}_{med}\)。

(3) 均匀损失价函数的贝叶斯估计：

当损失函数为均匀损失函数时，验后风险为 \[\begin{aligned} r_{unf}(\hat{\bm{X}}\mid\bm{Z})&=\int_{-\infty}^{\hat{X}-\frac{\Delta}{2}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} +\int_{\hat{X}+\frac{\Delta}{2}}^{+\infty}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x}\\ &=1-\int_{\hat{X}-\frac{\Delta}{2}}^{\hat{X}+\frac{\Delta}{2}}p(\bm{x}\mid\bm{z})\,\mathrm{d}\bm{x} \end{aligned} \tag{2.7.20}\] 将 \(r_{unf}(\hat{\bm{X}}\mid\bm{Z})\) 对 \(\hat{\bm{X}}\) 求导并使导数为零，得到 \[p\left(\left.\left(\hat{\bm{X}}+\frac{\Delta}{2}\right)\right|\bm{z}\right) -p\left(\left.\left(\hat{\bm{X}}-\frac{\Delta}{2}\right)\right|\bm{z}\right)=0 \tag{2.7.21}\] 显然，只有当 \(p(\bm{x}\mid\bm{z})\) 在 \(\hat{\bm{X}}\) 处有极大值时，才会有上式成立（如图 2.12 所示），即 \[\left.p(\bm{x}\mid\bm{z})\,\right|_{x=\hat{x}}=\max \tag{2.7.22}\] 这与极大验后估计的准则一致，这说明当损失函数为均匀损失价函数时，贝叶斯估计就是极大验后估计 \(\hat{\bm{X}}_{MAP}\)。

参数估计方法的相互关系

本章学习了统计论中的经典估计方法：最小二乘估计、极大似然估计、极大验后估计、最小方差估计和贝叶斯估计。这些估计方法都有各自的估计准则：最小二乘估计和最小方差估计是使损失函数最小的估计；极大似然和极大验后是以其相关的分布函数最大的估计。每种估计准则下得到的估计不同，但在某些情况下，不同估计准则下结果又是等价的。表 2.8 给出了这五种估计方法各自的准则，估计时需要的随机变量的先验信息和估计结果。从表中的各项比较可以看出：

(1) 最小二乘估计是使以残差定义的损失函数最小的估计，在估计时，将参数 \(\bm{X}\) 视为非随机量，只需要已知观测值与参数的函数关系和观测值误差的方差矩阵。

(2) 极大似然估计为：\(\hat{\bm{X}}_{ML}=\arg\max\limits_{\hat{\bm{x}}}p(\bm{z}\mid\hat{\bm{x}})\)；当观测值与参数呈线性关系并且正态分布时，极大似然估计与最小二乘估计等价。

(3) 极大验后估计为：\(\hat{\bm{X}}_{MAP}=\arg\max\limits_{\hat{\bm{x}}}p(\hat{\bm{x}}\mid\bm{z})\)；当 \(\bm{X}\) 和 \(\bm{Z}\) 都是正态分布时，估计值为 \(\hat{\bm{X}}_{MAP}=E(\bm{X}\mid\bm{Z})\)。

(4) 最小方差估计是使以估计误差定义的损失函数最小的估计，其估计为：\(\hat{\bm{X}}_{MV}=E(\bm{X}\mid\bm{Z})\)；线性最小方差估计是最小方差估计的特殊形式。

(5) 贝叶斯估计使估计误差造成的损失平均值最小，即使损失函数的验后期望最小的估计。不同的损失函数得到不同的估计，如当损失函数为平方损失函数时，贝叶斯估计即为最小方差估计等价；当损失函数为均匀损失函数时，贝叶斯估计与最大验后估计等价。

(6) 在所有的估计中，最小二乘估计不需要参数的任何先验信息；极大似然估计也不需要参数的任何先验信息，但需要已知观测值的概率分布函数；除了最小二乘估计和极大似然估计外，其他估计方法都将 \(\bm{X}\) 视为随机变量，在估计时需要已知 \(\bm{X}\) 的分布或者先验随机信息。

(7) 从是否为线性估计的角度来说，最小二乘估计和线性最小方差估计是线性估计，其他估计方法不一定是线性估计，也不一定有解析解，有时需要用数值方法得到满足其准则的估计。

本章介绍的经典最优估计方法之间的关系如图 2.13 所示。从本章的学习中也可以看出，贝叶斯估计并不是某一种估计方法，而是利用验后概率分布的一系列估计方法，它也可以推演出极大似然估计和最小二乘估计。当对 \(\bm{X}\) 毫无所知，没有任何先验信息的时候，可以简单地认为 \(\bm{X}\) 在其数值域内为均匀分布，这时极大验后估计准则 \[p(\bm{z}\mid\bm{x})p_x(\bm{x})=\max \tag{2.8.1}\] 就退化为 \[p(\bm{z}\mid\bm{x})=\max \tag{2.8.2}\] 满足上式的估计即为极大似然估计。

当 \(\bm{X}\) 为正态分布时，根据例 2.10，上式等价于 \[(\bm{z}-\bm{H}\bm{x})^{\mathrm{T}}\bm{D}_{\Delta}^{-1}(\bm{z}-\bm{H}\bm{x})=\min \tag{2.8.3}\] 即为最小二乘估计准则。