误差传播定律与最小二乘法


误差传播定律

误差传播(Error Propagation)定律指的是已知自变量之间的协方差矩阵,如何求得因变量之间的协方差矩阵。举例来讲,若 $\mathbf Y = \mathbf F(\mathbf X)$,其中 $\mathbf X$,$\mathbf Y$ 均为矢量,$\mathbf F$ 为矢量函数。已知自变量 $\mathbf X$ 的各个分量之间的协方差矩阵为 $\Sigma_X$,那么因变量 $\mathbf Y$ 的各个分量之间的协方差矩阵是什么?如果记因变量 $\mathbf Y$ 的各个分量之间的协方差矩阵为 $\Sigma_Y$,则
$$\Sigma_Y = \frac {\partial \mathbf F}{\partial \mathbf X}\Sigma_X \left[ \frac {\partial \mathbf F}{\partial \mathbf X} \right ]^T\;.$$

误差传播定律的推导

已知随机变量 $\boldsymbol \epsilon$ 的期望和协方差矩阵差满足 $$ E[\boldsymbol \epsilon]= 0, V[\boldsymbol \epsilon] = \Sigma_{\boldsymbol \epsilon}\;.$$ 通过协方差矩阵的定义,我们可得 $$V[\boldsymbol \epsilon] = E[(\boldsymbol \epsilon-E(\boldsymbol \epsilon))(\boldsymbol \epsilon-E(\boldsymbol \epsilon))^T] = E[\boldsymbol \epsilon {\boldsymbol \epsilon}^T]\;.$$ 变量 $\mathbf X+\boldsymbol \epsilon$ 的期望与协方差矩阵分别为 $$E[\mathbf X+\boldsymbol \epsilon] = E[\mathbf X]+E[\boldsymbol \epsilon] = \mathbf X \\ V[\mathbf X+\boldsymbol \epsilon] = E[(\mathbf X+\boldsymbol \epsilon-E[\mathbf X+\boldsymbol \epsilon])(\mathbf X+\boldsymbol \epsilon-E[\mathbf X+\boldsymbol \epsilon])^T] = E[\boldsymbol \epsilon {\boldsymbol \epsilon}^T]\;.$$ 变量 $\mathbf F(\mathbf X+\boldsymbol \epsilon)$ 的期望与协方差矩阵分别为 $$E[\mathbf F(\mathbf X+\boldsymbol \epsilon)] = E[\mathbf F(\mathbf X)+\dfrac{\partial \mathbf F}{\partial \mathbf X}\boldsymbol \epsilon] = \mathbf F(\mathbf X) \\ V[\mathbf F(\mathbf X+\boldsymbol \epsilon)] = E[(\mathbf F(\mathbf X+\boldsymbol \epsilon)-E[\mathbf F(\mathbf X+\boldsymbol \epsilon)])(\mathbf F(\mathbf X+\boldsymbol \epsilon)-E[\mathbf F(\mathbf X+\boldsymbol \epsilon)])^T] = \\ E[\dfrac{\partial \mathbf F}{\partial \mathbf X}\boldsymbol \epsilon (\dfrac{\partial \mathbf F}{\partial \mathbf X}\boldsymbol \epsilon)^T] = \dfrac{\partial \mathbf F}{\partial \mathbf X}\Sigma_{\boldsymbol \epsilon}(\dfrac{\partial \mathbf F}{\partial \mathbf X})^T\;.$$

用 uncertainties 计算变量的不确定度

uncertainties 是一个利用误差传播定律来计算变量不确定度、变量间协方差矩阵的 Python package。安装 uncertaintiespip install uncertainties。下面的代码给出了一些用 uncertainties 计算变量不确定度和协方差矩阵的实例。已知变量 $x =2\pm0.25$,求解 $x^2$:

输出结果为

计算 $sin(1+x^2)$ 在 $x =2\pm0.25$ 处的值:

输出结果为 -0.96+/-0.28。计算两个随机变量 $1\pm0.1$ 与 $2\pm0.2$ 之和、平均值以及余弦函数值:

输出结果为

已知两个随机变量 $u$ 和 $v$,构建一个新变量 $f(u,v) = u+2v$;若变量 $u$ 和 $v$ 分别取 $1\pm0.1$ 和 $10\pm0.2$,那么求 $f(u,v)$ 的值、不确定度以及偏导数:

输出结果为

求随机变量 $u$、$v$ 和新变量 $f(u,v)$ 的协方差矩阵:

输出结果为

从上面的结果可以看出,变量 $u$ 和 $v$ 之间的协方差为 0,而变量 $u$ 和 $f(u,v)$ 之间的协方差、变量 $v$ 和 $f(u,v)$ 之间的协方差不为零,这是因为变量 $f(u,v)$ 是变量 $u$ 和 $v$ 的函数。

最小二乘法

在简述最小二乘法之前,先回顾一下矩阵的微分运算。

Let the scalar α be defined by $$\alpha = \mathbf y^T \mathbf {Ax}\;,$$ where $\mathbf y$ is m×1, $\mathbf x$ is n×1, $\mathbf A$ is m×n, and $\mathbf A$ is independent of $\mathbf x$ and $\mathbf y$, then $$\frac{\partial \alpha}{\partial \mathbf x }=\mathbf y^T \mathbf A \\ \frac{\partial \alpha}{\partial \mathbf y }=\mathbf x^T \mathbf A^T\;.$$

最小二乘法简介

最小二乘法常用在通过观测量来拟合模型的参数上,而模型一般可用 $\mathbf Y = \mathbf F(\mathbf X,\boldsymbol \beta)$ 来表示,其中 $\mathbf Y$ 为观测量,$\mathbf X$ 为模型中涉及的变量,$\boldsymbol \beta$ 为模型中需要拟合的参数。以板块运动模型为例,$\mathbf Y$ 为 GPS 测站的东向和北向速度,$\mathbf X$ 为测站在地球上的位置(经纬度),$\boldsymbol \beta$ 为板块的欧拉矢量(待拟合参数)。

如何通过模型和观测量来求解最佳的模型参数?一种最直观的方法就是使观测值与模型值的残差平方和最小,此即最小二乘的基本原理。
令 $$S_i = \dfrac{1}{2}(\mathbf Y_i-\mathbf F(\mathbf X_i,\boldsymbol \beta))^T\Sigma_i^{-1}(\mathbf Y_i-\mathbf F(\mathbf X_i,\boldsymbol \beta))\;,$$ 其中 $\mathbf Y_i$ 代表第 $i$ 组观测量(每组观测量有 $k$ 个分量),$\Sigma_i$ 代表第 $i$ 组观测量的协方差矩阵,加权矩阵 $W_i = \Sigma_i^{-1}$,那么要得到最佳的模型参数 $\boldsymbol \beta$,则需要满足条件 $$\sum \limits_{i=1}^n S_i = min \;.$$ 上式等价于 $$\dfrac{\partial \sum \limits_{i=1}^n S_i}{\partial \boldsymbol \beta} = 0\;,$$ 根据矩阵的微分运算,上式化简为 $$\sum \limits_{i=1}^n (\mathbf Y_i-\mathbf F(\mathbf X_i,\boldsymbol \beta))^T \Sigma_i^{-1} \dfrac{\partial \mathbf F}{\partial \boldsymbol \beta} = 0$$。若模型是线性的,即 $$F(\mathbf X_i,\boldsymbol \beta) = A(\mathbf X_i)\boldsymbol \beta\;,$$ 那么将其带入上式并化简得 $$\widehat{\boldsymbol \beta} = \left[ \sum \limits_{i=1}^n A(\mathbf X_i)^T \Sigma_i^{-1} A(\mathbf X_i)\right]^{-1} \left[ \sum \limits_{i=1}^n A(\mathbf X_i)^T \Sigma_i^{-1} \mathbf Y_i\right]\;.$$ 通过误差传播定律可知模型参数之间的协方差矩阵 $\Sigma_{\boldsymbol \beta}$ 为 $$\Sigma_{\boldsymbol \beta} = \left[ \sum \limits_{i=1}^n A(\mathbf X_i)^T \Sigma_i^{-1} A(\mathbf X_i)\right]^{-1}\;,$$ 若第 $i$ 组观测量的协方差矩阵 $\Sigma_i$ 并未给出,则 $\Sigma_i$ 可统一用 $\sigma^2 \mathbf I[k \times k]$)代替,而 $\sigma^2$ 为下文中 WRMSE 的平方,即 $$\sigma^2 = \dfrac{\sum \limits_{i=1}^n \mathbf R_i^T \mathbf R_i}{kn-m} \;.$$

观测数据与模型估计值之间的残差(Residual)表示为 $$\mathbf R_i = \mathbf Y_i-\mathbf F(\mathbf X_i,\widehat{\boldsymbol \beta})\;,$$ 加权均方根差 WRMSE(Weighted Root Mean Square Error) 的表达式为 $$WRMSE = \sqrt{\dfrac{\sum \limits_{i=1}^n \mathbf R_i^T W_i \mathbf R_i}{\frac{kn-m}{kn}\sum \limits_{i=1}^n tr(W_i)}}\;,$$ 其中 $tr(\bullet)$、$m$ 分别为矩阵的迹和待拟合参数的个数。注意:这里的加权均方根差指的是分量的综合标准差,若估算总量的标准差,则需要乘以 $\sqrt k$。

最小二乘法实例

给定一组数据,拟合函数 $f(x) = a_0 + a_1 \dfrac{x}{T} + a_2\cos(\dfrac{2\pi}{T}x)+a_3\sin(\dfrac{2\pi}{T}x) $ 中的参数 $a_0$、$a_1$、$a_2$、$a_3$,其中 $T = 12$。曲线拟合的方法有很多,scipy.optimize 中提供了一种简单的方法 curve_fit

首先定义一个标准函数 $f(x) = 2 + \dfrac{x}{T} + 3\cos(\dfrac{2\pi}{12}x)+4\sin(\dfrac{2\pi}{12}x)$,然后根据标准函数生成 (x,y) 数据,并在 y 数据上添加随机扰动以及在特定位置添加较大的异常值(outlier)。

不迭代的最小二乘法

不迭代的最小二乘法指的是对观测数据中的异常值不做任何处理,直接对所有的观测数据进行拟合。

用不迭代的最小二乘法求解参数并作图

输出结果为

非迭代的最小二乘

因为异常值会对拟合的结果造成较大的偏差(观测数据中存在异常值会导致较大的残差;为了使残差的平方和最小,拟合的模型参数往往会更加偏向于异常值以减小残差的平方和),因此很有必要剔除观测数据中的异常值。剔除异常值的方法有多种,下面简单介绍两种常用的方法。

迭代的最小二乘法

迭代的最小二乘法指用观测数据对模型进行拟合得到模型参数,然后根据模型估计值和观测值之间的残差对异常值进行判定,例如残差大于 3 个 WRMSE 的观测数据则认为是异常值,因此需要对该异常值进行剔除。对剔除异常值后的观测数据重新进行最小二乘拟合,得到新的模型参数,并重新对所有的观测数据计算残差;依然对残差大于 3 个 WRMSE 的观测数据进行剔除,然后对剔除异常值后的观测数据进行拟合。不断重复上面的过程直到第 n 次保留的观测数据与第 n+1 次保留的观测数据一致时迭代结束,那么第 n 次拟合的模型参数为最终的模型参数。

用迭代的最小二乘法求解参数并作图

输出结果为

迭代的最小二乘法

通过与不迭代的最小二乘法作比较,用迭代最小二乘法求解的参数更加接近标准函数中的参数,并且标准差更小。

稳健的最小二乘法

损失函数(Loss Function)度量了观测值与模型估计值的偏离程度。上文描述的最小二乘法所采用的损失函数为 $$LF = \mathbf R_i^T \mathbf R_i\;.$$ 正如前面所讲,这种以残差平方形式的损失函数受异常值的影响很大,因此为了减小异常值对模型拟合的影响,我们需要采用一种降低异常值影响的损失函数。绝对值形式的损失函数为 $$LF = |\mathbf R_i|\;,$$ 尽管这种损失函数减小了异常值对模型拟合的影响,但是它难以微分,因此我们需要一种折中的损失函数既能保证异常值对模型拟合的影响较小又能保证它的可微分性。scipy.optimize.least_squares 提供了多种损失函数,例如 Huber 损失函数,其表达式为 $$LF = \left \lbrace \begin{array}{ll} \mathbf R_i^T \mathbf R_i & |\mathbf R_i| < 1 \\ 2 |\mathbf R_i|-1 & otherwise \end{array} \right. \;,$$;soft_l1 损失函数的表达式为 $$LF = 2(\sqrt{1+\mathbf R_i^T \mathbf R_i}-1) \;,$$ 这两种损失函数均能有效地降低异常值对模型参数的影响。通过引入损失函数的概念,我们引入稳健的最小二乘法(Robust Least Square)。稳健的最小二乘法指的是异常值对模型参数影响较小的一种拟合方法。下面采用 soft_l1 损失函数来进行稳健的最小二乘拟合。

输出结果为

通过比较非迭代和迭代最小二乘估计的结果,稳健最小二乘估计的结果与迭代最小二乘估计的结果基本一致,这说明稳健最小二乘法很大程度上降低了异常值对估计模型参数的影响。注意:两者定义函数的方式不一样,用 curve_fit 定义的函数是自变量在前,参数在后,且多个参数不能用数组表示;而 least_squares 定义的函数恰恰相反,参数在前,且多个参数必须用数组表示,自变量在后,且因变量也需要紧跟自变量之后。参考

考虑数据不确定度的最小二乘法

给定一组数据,拟合函数 $f(x) = a_0 + a_1 \dfrac{x}{T} $ 中的参数 $a_0$、$a_1$,其中 $T = 12$。

加权最小二乘法

构建函数 $f(x) = 2 + \dfrac{x}{T} $ ,对 $y$ 值添加随机波动和异常,设置其不确定度(标准差)服从正态分布 $\mathcal N(0.5,0.01)$。

定义加权最小二乘法函数

用加权最小二乘法求解参数并作图

输出结果为

迭代加权最小二乘

定义迭代加权最小二乘法函数。

对上面的数据用迭代加权最小二乘法进行拟合。

输出结果为

结果表明对于含有不确定度的数据,迭代的加权最小二乘法比非迭代的加权最小二乘法稍好,但并非绝对。同样,加权最小二乘法未必比非加权的最小二乘法拟合结果好,这跟数据的不确定度分布有密切关系。

病态的最小二乘估计与 Tikhonov 正规化

对于线性模型 $$\mathbf y = \mathbf {Ax}+ \boldsymbol \epsilon\;,$$待估参数 $$\mathbf x= [\mathbf A^T\mathbf A]^{-1}\mathbf A^T \mathbf y\;,$$其中 $\mathbf y$ 为观测量,$\mathbf A$ 为设计矩阵,$\boldsymbol \epsilon$ 为满足正态分布的误差量。给定一组数据,用 numpy.linalg.lstsq 拟合 $$ \mathbf y = a_1 \mathbf x + a_0 $$ 中的参数 $a_1$、$a_0$。

输出结果为

病态的最小二乘拟合

很显然对于病态的的数据分布,常规的最小二乘拟合得到了不正确的拟合结果,针对这种情况,这里采用 Tikhonov 正规化法对常规的最小二乘拟合进行改进。Tikhonov 正规化指通过添加拉格朗日乘子 $\lambda$,使得$$\parallel \mathbf A \mathbf x-\mathbf y\parallel^2+\lambda\parallel \mathbf x \parallel^2=min\;,$$待估参数$$\mathbf x= [\mathbf A^T\mathbf A + \lambda \mathbf I]^{-1}\mathbf A^T \mathbf y\;,$$其中 $\mathbf I$ 为单位矩阵。$\lambda$ 可通过 L-Curve 法求解,其基本原理如下:

  • 绘制 L-Curve 图:$z = f(u)$
  • 绘制 L-Curve 的曲率图:$\kappa=\frac{\left|z^{\prime \prime}\right|}{\left(1+z^{\prime 2}\right)^{3 / 2}}$
  • 计算最大曲率所对应的 $\lambda$

用自定义函数 L_Curve(A,y) 解算最佳的拉格朗日乘子 $\lambda$ 与待估参数 $\mathbf x$,并作 L-Curve 图和 L-Curve 的曲率图。

输出结果为

L-Curve 及其曲率
Tikhonov 正规化的最小二乘拟合

上图表明用 Tikhonov 正规化能够很好地解决病态的最小二乘问题。当 $\lambda$ 为零时,病态的最小二乘问题转换为常规的最小二乘问题;$\lambda$ 偏离 0 越远则表明数据越病态;上面的结果显示 $\lambda$ 的值为 2152.615,远大于零,说明数据存在很明显的异常,此点也可从图像得到验证。

多项式拟合

给定一组数据,用 numpy.polyfit 线性拟合 $ a_1x + a_0 $ 中的参数 $a_1$、$a_0$。

输出结果为 [2.04608068 0.10155647]

给定一组数据,用 numpy.polyfit 拟合多项式 $ a_0x^2 + a_1x +a_2 $ 中的参数 $a_0$、$a_1$、$a_2$。

输出结果为 [ 1.00686224 -3.41221002 -12.63385556]


参考

关于 “误差传播定律与最小二乘法” 的 2 个意见

发表评论

电子邮件地址不会被公开。 必填项已用*标注