问题描述
假设我们有一组观测值： \[ \underline{y_0} = \begin{bmatrix}y_{01}\\...\\y_{0i}\\...\\y_{0m}\end{bmatrix} \in \mathbb{R}^m \] 观测值之间相互独立。
我们想从中估算出一组未知参数： \[ \underline{x} = \begin{bmatrix}x_1\\...\\x_i\\...\\x_n\end{bmatrix} \in \mathbb{R}^n \quad (n \le m \text{ for redundancy}) \]
最小二乘估计（Least Squares Estimation）
最小二乘法，又称最小平方法，是一种数学优化建模方法。它通过最小化误差的平方和寻找数据的最佳函数匹配。 利用最小二乘法可以简便的求得未知的数据，并使得求得的数据与实际数据之间误差的平方和为最小。 “最小二乘法”是对线性方程组，即方程个数比未知数更多的方程组，以回归分析求得近似解的标准方法。
—— 维基百科
数据模型
假设 \(\underline{y_0}\) 是 \(\underline{y}\) 的一组观测值，而 \(\underline{y}\) 和 \(\underline{x}\) 之间为线性关系： \[\begin{aligned} \underline{y} = A\underline{x}+\underline{b} \\ \text{其中, } A: \text{ design matrix} \quad \underline{b}: \text{ known term} \end{aligned}\]
假设 \(\underline{y_0}\) 和 \(\underline{y}\) 之间存在误差 \(\underline{\varepsilon}\)，即： \[\begin{aligned} \underline{y_0} = \underline{y} + \underline{\varepsilon} = A\underline{x}+\underline{b} + \underline{\varepsilon}\\ \text{其中, }\underline{\varepsilon}: \text{ unknown error} \\ \end{aligned}\]
根据中心极限定理以及大数定理，假设： \[\begin{aligned} \underline{\varepsilon} \backsim N(0, \sigma_0^2Q) \\ \text{其中, } \sigma_0^2: \text{a prior variance} \quad Q: \text{ cofactor matrix} \end{aligned}\]
因此，问题转换为，求得 \(\hat{\underline{x}}\) 和 \(\hat{\underline{y}}\)，使其满足：
\[\begin{aligned} \left\{\begin{matrix} \hat{\underline{y}}=A\hat{\underline{x}}+\underline{b}\\ \hat{\underline{y}} \textbf{ at minimal distance from } \underline{y_0} \Rightarrow (\hat{\underline{y}}-\underline{y_0})^TQ^{-1}(\hat{\underline{y}}-\underline{y_0}): minimum \quad (\text{stochastic distance}) \end{matrix}\right. \end{aligned}\]
非线性场景
如果 \(\underline{y}\) 和 \(\underline{x}\) 之间为非线性关系，即 \(y = g(x)\)，则首先，我们要将其线性化：
- 通过求解观测值等式来获取未知参数的近似值 \(\underline{\widetilde{x}}\)
- 计算雅可比矩阵(Jacobian Matrix)： \[ J(\widetilde{x}) = \begin{bmatrix} \frac{\partial g_1}{\partial x_1} & \frac{\partial g_1}{\partial x_2} & ... & \frac{\partial g_1}{\partial x_n}\\ \frac{\partial g_2}{\partial x_1} & \frac{\partial g_2}{\partial x_2} & ... & \frac{\partial g_2}{\partial x_n}\\ ... & ... & ... & ...\\ \frac{\partial g_m}{\partial x_1} & \frac{\partial g_m}{\partial x_2} & ... & \frac{\partial g_m}{\partial x_n} \end{bmatrix}_{\underline{\widetilde{x}}} \\ \]
- 模型转换成： \[\begin{aligned} y_0 = f(\widetilde{x}) + J(\widetilde{x})(x-\widetilde{x}) \\ \text{其中， } A = J, b = f(\widetilde{x}) \end{aligned}\]
求解
说明：求最小值问题可以转换为求一阶导数为 0 的问题。这里省略推导过程，直接给出计算步骤.
- 计算正规矩阵（Normal Matrix）：\(N = A^TQ^{-1}A\)
- 计算正规已知项（Normal Known Term）：\(\underline{t} = A^TQ^{-1}(\underline{y}_0 - \underline{b})\)
- 计算未知参数的估算解： \[N\hat{\underline{x}} = A^TQ^{-1}A\hat{\underline{x}} = A^TQ^{-1}(\underline{y}_0 - \underline{b}) = \underline{t} \quad \Rightarrow \quad \underline{\hat{x}} = N^{-1}\underline{t} = N^{-1}A^TQ^{-1}(\underline{y}_0 - \underline{b}) \qquad (\text{the normal system})\]
- 计算观测值的估算值：\(\underline{\hat{y}} = A\underline{\hat{x}} + \underline{b}\)
- 计算误差的估算值：\(\underline{\hat{\varepsilon}} = \underline{y}_0 - \underline{\hat{y}}\)
注意： 1. \(N\) 必须是可逆的，这就要求 \(A\) 满秩 ## 精确度估算 ### 方差 上一步中，我们得到了误差的估算值，因此可以计算后验方差（posterior variance）： \[ \hat{\sigma}^2 = \frac{\underline{\hat{\varepsilon}}^TQ^{-1}\underline{\hat{\varepsilon}}}{m-n} \]
协方差
根据 协方差传播（covariance propagation）：
\[\begin{aligned} \textbf{在线性的场景下：} y = Ax+b \\ \text{则 }\left\{\begin{matrix}\text{均值 } \mu_y = A\mu_x + b \\ \text{协方差 } C_{yy} = AC_{xx}A^T \end{matrix}\right.\\ \quad \\ \textbf{在非线性的场景下：} y = g(x) \\ \text{则 }\left\{\begin{matrix}\text{均值 } \mu_y \cong g(\mu_x) \\ \text{协方差 } C_{yy} = JC_{xx}J^T \end{matrix}\right.\\ \end{aligned}\]
因此观测值的协方差为： \[\begin{aligned} C_{y0y0} = IC_{\varepsilon\varepsilon}I^T = C_{\varepsilon\varepsilon} = \sigma_0^2Q \end{aligned}\]
由于前面已经得出一个后验方差，因此，上面又可以写成： \[\begin{aligned} C_{y0y0} = \hat{\sigma}^2Q \end{aligned}\]
未知参数的估算解的协方差
\[\begin{aligned} \because \underline{\hat{x}} = N^{-1}A^TQ^{-1}(\underline{y}_0 - \underline{b}) \\ \therefore C_{\hat{x}\hat{x}} = N^{-1}A^TQ^{-1}C_{y_0y_0}(N^{-1}A^TQ^{-1})^T \\ = N^{-1}A^TQ^{-1}\hat{\sigma}^2Q(N^{-1}A^TQ^{-1})^T \\ = \hat{\sigma}^2N^{-1} \end{aligned}\]
观测值的估算值的协方差
\[\begin{aligned} \because \underline{\hat{y}} = A\underline{\hat{x}} + \underline{b} \\ \therefore C_{\hat{y}\hat{y}} = AC_{\hat{x}\hat{x}}A^T \\ = A\hat{\sigma}^2N^{-1}A^T \end{aligned}\]
误差的估算值的协方差
\[\begin{aligned} \because \underline{\hat{\varepsilon}} = \underline{y}_0 - \underline{\hat{y}} \\ \therefore C_{\hat{\varepsilon}\hat{\varepsilon}} = C_{y0y0} - C_{\hat{y}\hat{y}} \\ = \hat{\sigma}^2Q - A\hat{\sigma}^2N^{-1}A^T \\ = \hat{\sigma}^2(Q-AN^{-1}A^T) \end{aligned}\]
置信度评估
全局测试
- 做出假设 \[H_0: \hat{\sigma}^2=\sigma^2\]
- 计算：\(\frac{\hat{\sigma}^2}{\sigma^2}(m-n)\backsim\chi^2_{exp}\)
- 查询卡方分布表：\(\chi^2_{lim}=\chi^2_{(m-n)}(\alpha)\)
- 如果 \(\chi^2_{exp}\le\chi^2_{lim}\)，则接受假设；如果 \(\chi^2_{exp}\gt\chi^2_{lim}\)，则拒绝假设。
局部测试
用以识别异常值（outliers）。
- 做出假设 \[\frac{\hat{\varepsilon_i}}{\sigma_{\varepsilon_i}} = \tau_{exp} \sim \tau_{(m-n)}\]
- 计算 \(\tau_{exp}\)，其中 \(\sigma_{\varepsilon_i}=\hat{\sigma}\sqrt{C_{\hat{\varepsilon}\hat{\varepsilon}}(i,i)}\)
- 查询 T 分布表：\(\tau_{lim}=\tau_{(m-n)}(\frac{\alpha}{2})\)
- 如果 \(|\tau_{exp}|\le\tau^2_{lim}\)，则接受假设；如果 \(|\tau_{exp}|\gt\tau^2_{lim}\)，则拒绝假设。
- 不满足假设的观测值即为异常者。将其移除后，重复整个过程。
