在Scipy中，curve_fit如何及为何计算参数估计的协方差

Question

我一直在使用 scipy.optimize.leastsq 来拟合一些数据。我想要得到这些估计值的置信区间，所以我查看了 cov_x 的输出，但文档上对这个内容的解释非常不清楚，我不知道该如何从中得到参数的协方差矩阵。

首先，它说这是一个雅可比矩阵，但在说明中也提到“cov_x 是对海森矩阵的雅可比近似”，所以这实际上不是一个雅可比矩阵，而是用雅可比的一些近似值得到的海森矩阵。这两种说法哪个是正确的呢？

其次，这句话让我感到困惑：

这个矩阵必须乘以残差方差，才能得到参数估计的协方差——请参见 curve_fit。

我确实去查看了 curve_fit 的源代码，他们做了：

s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

这对应于将 cov_x 乘以 s_sq，但我在任何参考资料中都找不到这个公式。有人能解释一下这个公式为什么是正确的吗？我的直觉告诉我应该是反过来的，因为 cov_x 应该是一个导数（雅可比或海森），所以我在想：cov_x * covariance(parameters) = sum of errors(residuals)，其中 sigma(parameters) 是我想要的东西。

我该如何将 curve_fit 所做的与我在维基百科上看到的内容联系起来呢？http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

Answer 1

数学基础

我们先从线性回归开始。在很多统计问题中，我们假设变量有一些潜在的分布，里面有一些未知的参数，我们需要估计这些参数。在线性回归中，我们假设因变量 y_i 和自变量 x_ij 之间有线性关系：

y_i = x_i1β₁ + ... + x_ipβ_p + σε_i, i = 1, ..., n。

这里，ε_i 是独立的标准正态分布，β_j 是 p 个未知参数，σ 也是未知的。我们可以把这个写成矩阵的形式：

Y = X β + σε。

在这里，Y、β 和 ε 都是列向量。为了找到最好的 β，我们需要最小化平方和：

S = (Y - X β)^T (Y - X β)。

我直接写出解决方案：

β^ = (X^T X)^-1 X^T Y。

如果我们把 Y 看作具体的观测数据，那么 β^ 就是根据这个观测得到的 β 的估计。另一方面，如果我们把 Y 看作随机变量，那么估计值 β^ 也会变成随机变量。这样，我们就可以看到 β^ 的协方差。

因为 Y 服从多元正态分布，而 β^ 是 Y 的线性变换，所以 β^ 也服从多元正态分布。β^ 的协方差矩阵是：

Cov(β^) = (X^T X)^-1 X^T Cov(Y) ((X^T X)^-1 X^T)^T = (X ^T X)^-1 σ²。

但这里的 σ 是未知的，所以我们也需要估计它。如果我们设：

Q = (Y - X β^)^T (Y - X β^)，

可以证明 Q / σ² 服从自由度为 n - p 的卡方分布（而且，Q 和 β^ 是独立的）。这使得：

σ^² = Q / (n - p)

成为 σ² 的无偏估计。因此，β^ 的最终协方差估计是：

(X^T X)^-1 Q / (n - p)。

SciPy API

curve_fit 是最方便的，第二个返回值 pcov 就是 β^ 的协方差估计，也就是上面提到的最终结果 (X^T X)^-1 Q / (n - p)。

在 leastsq 中，第二个返回值 cov_x 是 (X^T X)^-1。从 S 的表达式中，我们看到 X^T X 是 S 的海森矩阵（准确来说是海森矩阵的一半），这就是文档中说 cov_x 是海森矩阵的逆的原因。要得到协方差，你需要把 cov_x 乘以 Q / (n - p)。

非线性回归

在非线性回归中，y_i 与参数之间的关系是非线性的：

y_i = f(x_i, β₁, ... , β_p) + σε_i。

我们可以计算 f 对 β_j 的偏导数，这样就可以把它近似为线性。然后计算基本上和线性回归是一样的，只不过我们需要迭代地近似最小值。在实际中，算法可以是一些更复杂的，比如 Levenberg–Marquardt 算法，这是 curve_fit 的默认算法。

关于提供 Sigma 的更多信息

这一部分讲的是 curve_fit 中的 sigma 和 absolute_sigma 参数。如果你在使用 curve_fit 时对 Y 的协方差没有任何先验知识，可以忽略这一部分。

绝对 Sigma

在上面的线性回归中，y_i 的方差是 σ，且是未知的。如果你知道方差，可以通过 sigma 参数提供给 curve_fit，并设置 absolute_sigma=True。

假设你提供的 sigma 矩阵是 Σ。即：

Y ~ N(X β, Σ)。

Y 服从均值为 X β 和协方差为 Σ 的多元正态分布。我们想要最大化 Y 的似然性。从 Y 的概率密度函数来看，这等价于最小化：

S = (Y - X β)^T Σ^-1 (Y - X β)。

解决方案是：

β^ = (X^T Σ^-1 X)^-1 X^T Σ^-1 Y。

而：

Cov(β^) = (X^T Σ^-1 X)^-1。

上面的 β^ 和 Cov(β^) 是 curve_fit 在 absolute_sigma=True 时的返回值。

相对 Sigma

在某些情况下，你不知道 y_i 的确切方差，但你知道不同 y_i 之间的相对关系，比如 y₂ 的方差是 y₁ 方差的 4 倍。这时你可以传递 sigma 并设置 absolute_sigma=False。

这时：

Y ~ N(X β, Σσ)

提供一个已知的矩阵 Σ 和一个未知的数 σ。要最小化的目标函数和绝对 sigma 的情况是一样的，因为 σ 是一个常数，因此估计值 β^ 也是一样的。但协方差：

Cov(β^) = (X^T Σ^-1 X)^-1 σ²，

里面有未知的 σ。为了估计 σ，我们设：

Q = (Y - X β^)^T Σ^-1 (Y - X β^)。

同样，Q / σ² 服从自由度为 n - p 的卡方分布。

Cov(β^) 的估计是：

(X^T Σ^-1 X)^-1 Q / (n - p)。

这就是 curve_fit 在 absolute_sigma=False 时的第二个返回值。

Answer 2

我在寻找类似问题的过程中发现了这个解决方案，对HansHarhoff的回答有一点小改进。使用leastsq函数的完整输出会返回一个叫做infodict的值，这里面包含了infodict['fvec'] = f(x) - y。这样，我们就可以计算简化的卡方值，公式是：

s_sq = (infodict['fvec']**2).sum()/ (N-n)

顺便说一下，感谢HansHarhoff为解决这个问题做了大部分的工作。

Answer 3

好的，我想我找到了答案。首先，解决方案是：cov_x*s_sq实际上就是你想要的参数的协方差。对对角线上的元素取平方根会得到标准差（但要小心协方差哦！）。

残差方差 = 减少的卡方 = s_sq = sum[(f(x)-y)^2]/(N-n)，这里的N是数据点的数量，n是拟合参数的数量。减少的卡方。

我困惑的原因是，leastsq给出的cov_x实际上在其他地方并不叫做cov(x)，而是减少的cov(x)或分数cov(x)。它在其他参考资料中没有出现的原因是，这只是一个简单的重新缩放，这在数值计算中很有用，但在教科书中并不相关。

关于Hessian和Jacobian，文档的表述不太清楚。实际上在这两种情况下计算的是Hessian，因为在最小值处Jacobian是零。他们的意思是，他们使用了一种近似的Jacobian来找到Hessian。

还有一点。似乎curve_fit的结果并没有真正考虑错误的绝对大小，而只是考虑了提供的sigma的相对大小。这意味着即使误差条的大小变化了一百万倍，返回的pcov也不会改变。当然，这样是不对的，但似乎是标准做法，比如Matlab在使用他们的曲线拟合工具箱时也是这样。正确的程序在这里描述：https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

一旦找到最佳值，这个过程似乎相当简单，至少对于线性最小二乘法来说是这样。