python - 在 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

好的，我想我找到了答案。首先是解决方案： cov_x*s_sq 只是您想要的参数的协方差。取对角线元素的 sqrt 会给你标准偏差（但要小心协方差！）。

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

我感到困惑的原因是，由 minimumsq 给出的 cov_x 实际上并不是其他地方所谓的 cov(x)，而是简化的 cov(x) 或分数 cov(x)。它没有出现在任何其他参考资料中的原因是它是一种简单的重新缩放，在数值计算中很有用，但与教科书无关。

关于 Hessian 与 Jacobian，文档措辞不佳。在这两种情况下计算的都是 Hessian 矩阵，这是显而易见的，因为 Jacobian 矩阵至少为零。他们的意思是他们正在使用雅可比矩阵的近似值来找到黑森矩阵。

进一步说明。似乎curve_fit结果实际上并没有考虑到错误的绝对大小，而只是考虑了所提供的sigma的相对大小。这意味着即使错误栏改变了一百万倍，返回的 pcov 也不会改变。这当然是不对的，但似乎是标准做法，即。Matlab 在使用他们的曲线拟合工具箱时做同样的事情。此处描述了正确的程序：https ://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

一旦找到最优值，至少对于线性最小二乘而言，这样做似乎相当简单。

Answer 2

我在搜索类似问题的过程中找到了这个解决方案，我对 HansHarhoff 的回答只有一点点改进。leastsq 的完整输出提供了一个返回值 infodict，其中包含 infodict['fvec'] = f(x) -y。因此，要计算减少的卡方 =（在上述符号中）

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

顺便提一句。感谢 HansHarhoff 完成了大部分繁重的工作来解决这个问题。

Answer 3

数学

首先，我们从线性回归开始。在许多统计问题中，我们假设变量具有一些带有一些未知参数的潜在分布，并且我们估计这些参数。在线性回归中，我们假设因变量 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)

^{σ 2}的无偏估计量。所以 Cov(β^) 的最终估计量是

(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 的 Hessian（准确地说是 Hessian 的一半），这就是为什么文档说cov_x是 Hessian 的逆。要获得协方差，您需要乘以cov_xQ / (n - p)。

非线性回归

在非线性回归中，y _i非线性地取决于参数：

y _i = f( x _i , β ₁ , ... , β _p ) + σε _i。

我们可以计算 f 对 β _j的偏导数，所以它变成近似线性的。然后计算与线性回归基本相同，只是我们需要迭代地逼近最小值。在实践中，该算法可以是一些更复杂的算法，例如默认的 Levenberg-Marquardt 算法curve_fit。

更多关于提供 Sigma

本节介绍 中的sigmaandabsolute_sigma参数curve_fit。对于curve_fit当您对 Y 的协方差没有先验知识时的基本用法，您可以忽略此部分。

绝对西格玛

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

假设您提供的sigma矩阵是 Σ。IE

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。

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

相对西格玛

在某些情况下，您不知道 y _{i的确切方差，但您知道不同 y i}之间的相对关系，例如 y _{2的方差是 y 1}的方差的 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_fitwith的第二个返回值absolute_sigma=False。