Han*_*off 31 python curve scipy
我一直在scipy.optimize.leastsq用来拟合一些数据.我想在这些估计值上得到一些置信区间,所以我查看cov_x输出但是文档很清楚这是什么以及如何从中得到我的参数的协方差矩阵.
首先,它说它是一个雅可比行列式,但在笔记中它也说" cov_x是雅各比式的雅可比近似",所以它实际上不是雅可比,而是一个使用雅可比的近似的Hessian.以下哪些陈述是正确的？
其次这句话对我来说很困惑:
该矩阵必须乘以残差方差,以得到参数估计的协方差 - 见
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
Han*_*off 27
好的,我想我找到了答案.首先解决方案:cov_x*s_sq只是参数的协方差,这是你想要的.取对角元素的sqrt将给出标准偏差(但要小心协方差!).
残差方差=减小的chi square = s_sq = sum [(f(x)-y)^ 2] /(Nn),其中N是数据点的数量,n是拟合参数的数量.卡方减少.
我混淆的原因是,由minimalsq给出的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
一旦找到最佳值,这似乎相当简单,至少对于线性最小二乘法.
我在寻找类似问题的过程中找到了这个解决方案,而HansHarhoff的答案我只有一点点改进.来自leastsq的完整输出提供了一个返回值infodict,其中包含infodict ['fvec'] = f(x)-y.因此,计算减小的卡方=(在上面的表示法中)
s_sq = (infodict['fvec']**2).sum()/ (N-n)
BTW.感谢HansHarhoff为解决这个问题做了大部分繁重的工作.
首先我们从线性回归开始。在许多统计问题中，我们假设变量具有一些具有某些未知参数的潜在分布，并且我们估计这些参数。在线性回归中，我们假设因变量 y i与自变量 x ij具有线性关系：
\ny i = x i1 \xce\xb2 1 + ... + x ip \xce\xb2 p + \xcf\x83\xce\xb5 i , i = 1, ..., n。
\n其中 \xce\xb5 i具有独立的标准正态分布， \xce\xb2 j是 p 个未知参数， \xcf\x83 也是未知的。我们可以将其写成矩阵形式：
\nY = X \xce\xb2 + \xcf\x83\xce\xb5,
\n其中 Y、\xce\xb2 和 \xce\xb5 是列向量。为了找到最好的\xce\xb2，我们最小化平方和
\nS = (Y - X \xce\xb2) T (Y - X \xce\xb2)。
\n我只是写出解决方案，即
\n\xce\xb2^ = (X T X) -1 X T Y。
\n如果我们将 Y 视为特定的观测数据，则 \xce\xb2^ 是该观测下 \xce\xb2 的估计。另一方面，如果我们将 Y 视为随机变量，则估计器 \xce\xb2^ 也成为随机变量。这样我们就可以看出\xce\xb2^的协方差是多少。
\n由于 Y 具有多元正态分布，并且 \xce\xb2^ 是 Y 的线性变换，因此 \xce\xb2^ 也具有多元正态分布。\xce\xb2^ 的协方差矩阵为
\nCov(\xce\xb2^) = (X T X) -1 X T Cov(Y) ((X T X) -1 X T ) T = (X T X) -1 \xcf\x83 2。
\n但这里\xcf\x83是未知的，所以我们还需要估计它。如果我们让
\nQ = (Y - X \xce\xb2^) T (Y - X \xce\xb2^),
\n可以证明Q / \xcf\x83 2具有n - p 自由度的卡方分布（而且Q 与\xce\xb2^ 无关）。这使得
\n\xcf\x83^ 2 = Q / (n - p)
\n\xcf\x83 2的无偏估计量。所以 Cov(\xce\xb2^) 的最终估计量是
\n(X T X) -1 Q / (n - p)。
\ncurve_fit最方便的是，第二个返回值pcov只是 \xce\xb2^ 协方差的估计，即上面最终的结果 (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非线性地依赖于参数：
\ny i = f( x i , \xce\xb2 1 , ... , \xce\xb2 p ) + \xcf\x83\xce\xb5 i。
\n我们可以计算 f 对 \xce\xb2 j的偏导数，因此它变得近似线性。那么计算基本上与线性回归相同，只是我们需要迭代地逼近最小值。实际上，该算法可以是一些更复杂的算法，例如 Levenberg\xe2\x80\x93Marquardt 算法，它是 的默认算法curve_fit。
本节介绍中的sigma和参数。对于基本用法，当您对 Y 的协方差没有先验知识时，可以忽略本节。absolute_sigmacurve_fitcurve_fit
在上面的线性回归中， y i的方差为 \xcf\x83 并且未知。如果你知道方差。curve_fit您可以通过参数提供sigma并设置absolute_sigma=True。
假设您提供的sigma矩阵是 \xce\xa3。IE
Y ~ N(X \xce\xb2, \xce\xa3)。
\nY 具有均值 X \xce\xb2 和协方差 \xce\xa3 的多元正态分布。我们想要最大化 Y 的可能性。根据 Y 的概率密度函数，这相当于最小化
\nS = (Y - X \xce\xb2) T \xce\xa3 -1 (Y - X \xce\xb2)。
\n解决办法是
\n\xce\xb2^ = (X T \xce\xa3 -1 X) -1 X T \xce\xa3 -1 Y。
\n和
\nCov(\xce\xb2^) = (X T \xce\xa3 -1 X) -1。
\ncurve_fit上面的 \xce\xb2^ 和 Cov(\xce\xb2^) 是with的返回值absolute_sigma=True。
在某些情况下，您不知道 y i的确切方差，但知道不同 y i之间的相对关系，例如 y 2的方差是 y 1方差的 4 倍。然后就可以通过sigma并设置了absolute_sigma=False。
这次
\nY ~ N(X \xce\xb2, \xce\xa3\xcf\x83)
\n提供已知矩阵 \xce\xa3 和未知数 \xcf\x83。最小化的目标函数与绝对 sigma 相同，因为 \xcf\x83 是常数，因此估计器 \xce\xb2^ 是相同的。但协方差
\nCov(\xce\xb2^) = (X T \xce\xa3 -1 X) -1 \xcf\x83 2 ,
\n其中有未知的\xcf\x83。要估计 \xcf\x83，让
\nQ = (Y - X \xce\xb2^) T \xce\xa3 -1 (Y - X \xce\xb2^)。
\n同样，Q / \xcf\x83 2具有 n - p 自由度的卡方分布。
\nCov(\xce\xb2^) 的估计为
\n(X T \xce\xa3 -1 X) -1 Q / (n - p)。
\n这是curve_fitwith的第二个返回值absolute_sigma=False。
| 归档时间： |
|
| 查看次数： |
18435 次 |
| 最近记录： |