在 Scipy 中，curve_fit 如何以及为什么计算参数估计的协方差 [英] In Scipy how and why does curve_fit calculate the covariance of the parameter estimates

问题描述

我一直在使用 scipy.optimize.leastsq 来拟合一些数据.我想获得这些估计的一些置信区间，因此我查看了 cov_x 输出，但文档非常不清楚这是什么以及如何从中获取参数的协方差矩阵.

I have been using scipy.optimize.leastsq to fit some data. I would like to get some confidence intervals on these estimates so I look into the cov_x output but the documentation is very unclear as to what this is and how to get the covariance matrix for my parameters from this.

首先它说它是雅可比行列式，但在 注释 它还说cov_x 是 Hessian 的雅可比近似值"，因此它实际上不是雅可比矩阵，而是使用雅可比矩阵的一些近似值的黑森矩阵.以下哪些说法是正确的?

First of all it says that it is a Jacobian, but in the notes it also says that "cov_x is a Jacobian approximation to the Hessian" so that it is not actually a Jacobian but a Hessian using some approximation from the Jacobian. Which of these statements is correct?

其次这句话让我很困惑:

Secondly this sentence to me is confusing:

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

This matrix must be multiplied by the residual variance to get the covariance of the parameter estimates – see curve_fit.

我确实去看看他们所做的 curve_fit 的源代码:

I indeed go look at the source code for curve_fit where they do:

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

对应于 cov_x 乘以 s_sq 但我在任何参考资料中都找不到这个方程.有人能解释一下为什么这个方程是正确的吗?我的直觉告诉我它应该是相反的，因为 cov_x 应该是衍生物(Jacobian 或 Hessian)所以我在想:cov_x * covariance(parameters) = sum of errors(residuals) 其中 sigma(parameters) 是我想要的东西.

which corresponds to multiplying cov_x by s_sq but I cannot find this equation in any reference. Can someone explain why this equation is correct? My intuition tells me that it should be the other way around since cov_x is supposed to be a derivative (Jacobian or Hessian) so I was thinking: cov_x * covariance(parameters) = sum of errors(residuals) where sigma(parameters) is the thing I want.

我如何将curve_fit正在做的事情与我所看到的联系起来，例如.维基百科:http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

How do I connect the thing curve_fit is doing with what I see at eg. wikipedia: http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

推荐答案

好的，我想我找到了答案.先说解决办法:cov_x*s_sq 只是您想要的参数的协方差.取对角线元素的 sqrt 会给你标准偏差(但要注意协方差！).

OK, I think I found the answer. First the solution: cov_x*s_sq is simply the covariance of the parameters which is what you want. Taking sqrt of the diagonal elements will give you standard deviation (but be careful about covariances!).

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

Residual variance = reduced chi square = s_sq = sum[(f(x)-y)^2]/(N-n), where N is number of data points and n is the number of fitting parameters. Reduced chi square.

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

The reason for my confusion is that cov_x as given by leastsq is not actually what is called cov(x) in other places rather it is the reduced cov(x) or fractional cov(x). The reason it does not show up in any of the other references is that it is a simple rescaling which is useful in numerical computations, but is not relevant for a textbook.

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

About Hessian versus Jacobian, the documentation is poorly worded. It is the Hessian that is calculated in both cases as is obvious since the Jacobian is zero at a minimum. What they mean is that they are using an approximation to the Jacobian to find the Hessian.

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

A further note. It seems that the curve_fit result does not actually account for the absolute size of the errors, but only take into account the relative size of the sigmas provided. This means that the pcov returned doesn't change even if the errorbars change by a factor of a million. This is of course not right, but seems to be standard practice ie. Matlab does the same thing when using their Curve fitting toolbox. The correct procedure is described here: https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

一旦找到最优值，就可以很容易地做到这一点，至少对于线性最小二乘法来说是这样.

It seems fairly straightforward to do this once the optimum has been found, at least for Linear Least squares.

这篇关于在 Scipy 中，curve_fit 如何以及为什么计算参数估计的协方差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！