在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，因此它实际上不是Jacobian而是使用与Jacobian近似的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 。

我确实去看看 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 应该是派生的（雅可比或黑森州），所以我在想： b cov_x *协方差（参数）=错误总和（残差）其中， sigma（参数）是我想要的东西

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就是您想要的参数的协方差。取对角线元素的平方根会给您标准偏差（但要注意协方差！）。

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] /（Nn），其中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.

原因因为我的困惑是，由lesssqq给出的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最小为零。他们的意思是说他们在使用近似于雅可比行列式来找到黑森州。

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屋！