Scipy curve_fit 没有给出合理的结果答案

【问题标题】：Scipy's curve_fit not giving reasonable resultScipy curve_fit 没有给出合理的结果
【发布时间】：2014-03-23 21:47:53
【问题描述】：

至少乍一看，我有一个简单的x,y 数据集适合。问题是scipy.optimize.curve_fit 为拟合的参数之一返回了一个 very 大的值，我不知道这在数学上是否正确，或者我拟合数据的方式是否有问题.

下图以蓝色显示了数据点和获得的最佳拟合。使用的曲线（下面MWE 中的func）有四个参数a, b, c, d 正在拟合：

a 大约给出了曲线达到一半最大值时的 x 值。
b 表示曲线稳定处的x 值。这个func 值由d 参数给出，即：func(b) = d
c 与曲线在原点的最大值有关：func(0) = c*constant + d
d 是曲线稳定的地方（图中黑线）。

b 参数是我遇到的问题（见问题结尾），它也是我最有兴趣分配合理值的参数。

MWE 显示正在拟合的函数和结果：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Function to be fitted.
def func(x, a, b, c, d):
    return c * (1 / np.sqrt(1 + (np.asarray(x) / a) ** 2) -
        1 / np.sqrt(1 + (b / a) ** 2)) ** 2 + d

# Define x,y data.    
x_list = [12.5, 37.5, 62.5, 87.5, 112.5, 137.5, 162.5, 187.5, 212.5, 237.5,
    262.5, 287.5, 312.5, 337.5, 362.5, 387.5, 412.5, 437.5, 462.5, 487.5,
    512.5]
y_list = [0.008, 0.0048, 0.0032, 0.00327, 0.0023, 0.00212, 0.00187,
    0.00086, 0.00070, 0.00100, 0.00056, 0.00076, 0.00052, 0.00077, 0.00067,
    0.00048, 0.00078, 0.00067, 0.00069, 0.00061, 0.00047]

# Initial guess for the 4 parameters.
guess = (50., 200., 80. / 10000., 6. / 10000.)

# Fit curve to x,y data.
f_prof, f_err = curve_fit(func, x_list, y_list, guess)

# Values for the a,b,c,d fitted parameters.
print f_prof

# Errors (standard deviations) for the fitted parameters.
print np.sqrt(f_err[0][0]), np.sqrt(f_err[1][1]), np.sqrt(f_err[2][2]),\
    np.sqrt(f_err[3][3])

# Generate plot.
plt.scatter(x_list, y_list)
plt.plot(x_list, func(x_list, f_prof[0], f_prof[1], f_prof[2], f_prof[3]))
plt.hlines(y=f_prof[3], xmin=0., xmax=max(x_list))
plt.show()

我得到的结果是：

# a, b, c, d
 52.74, 2.52e+09, 7.46e-03, 5.69e-04

# errors
11.52, 1.53e+16, 0.0028, 0.00042

b 参数值很大，错误也很大。通过查看图中绘制的数据，可以通过肉眼估计b 的值（即：数据集稳定的x 值）应该在x=300 附近。为什么b 和它的错误值这么大？

【问题讨论】：

恐怕问题出在您正在拟合的数据（或函数）上。首先，您的错误中有一个错字（应该是 f_err[1][1]）。我已经使用全局求解器检查了您的示例，它在其他参数中给出了相同的结果，但我什至得到 b = 3.56564242e+18。但毫不奇怪，点 x = 300 附近的导数是一个真正的蜗牛——它接近最小值但非常慢。您可以手动检查 - 通过参数对数据计算偏导数，并尝试求解 4 个非线性方程组。
@Martin 感谢您指出错字，我已经修复它（虽然值很好）
我知道。我已经检查过了。

标签： python numpy scipy curve-fitting curve

【解决方案1】：

我不知道这是故意的还是错误的，但在我看来，“b”将与“a”和“d”密切相关，并且与自变量“x”没有“相互作用”。如果 b/a 足够大，您可以将 1/np.sqrt(1 + (b / a) ** 2)) ** 2 近似为 a/b，这样您的函数就变为 c * function_of(x, a) - a/b + d

你的 'a' 和 'x' 值足够大，这变得非常接近 c*a/x - a/b + d。

正如 behzad.nouri 所指出的，与其他最小化器相比，curve_fit 可能稍微不稳定，并且总是最小化最小二乘。但它确实返回完整的协方差矩阵，包括变量之间的相关性（f_err 的非对角元素）。用这些！！

如果您确定 'b' 的值约为 300，或者有兴趣在 fmin 和 levenberg-marquardt 算法之间轻松切换，您可能会发现 lmfit 包 (http://lmfit.github.io/lmfit-py/) 很有用。它允许您对参数设置界限，在拟合算法之间轻松切换，还可以对参数的置信区间进行更强力的探索。

【讨论】：

感谢lmfitMatt 的提示。不幸的是，我尝试应用此处显示的所有最小化方法lmfit.github.io/lmfit-py/fitting.html#fit-engines-label，结果是相同的：参数b 非常大。即使通过最大值限制它也只会导致返回最大值，所以它没有用。再次感谢您 +1 的软件包推荐！
变量 a 和 b 之间有什么相关性？如果所有拟合方法都给出相同的结果，并且如果约束值导致限制总是被击中，那么这不是向您暗示问题是不适定的吗？同样，你基本上有 c*a/x - a/b + d。除非 'a' 和 'd' 以非常高的准确度已知，否则 'b' 中的不确定性将非常大......这就是你所看到的。

【解决方案2】：

您可以对参数的范数使用惩罚值，并使用fmin：

from scipy.optimize import fmin

def func(x, a, b, c, d):
    return c * (1 / np.sqrt(1 + (x / a) ** 2) - 1 / np.sqrt(1 + (b / a) ** 2)) ** 2 + d

def errfn(params, xs, ys, lm, ord=1):
    '''
    lm: penalty maltiplier
    ord: order in norm calculation
    '''
    from numpy.linalg import norm
    a, b, c, d = params
    err = func(xs, a, b, c, d) - ys
    return norm(err) + lm * norm(params, ord)

params = fmin(errfn, guess, args=(xs, ys, 1e-6, 2))

上面我使用了1e-6的小罚分，拟合结果是

[6.257e+01   3.956e+02   9.926e-03   7.550e-04]

合身：

编辑：使用惩罚函数和范数顺序，它非常适合

params = [  1.479e+01  -3.344e+00  -8.781e-03   8.347e-03]

【讨论】：

第二行中的b 参数是否为负值？这绝对不行，请阅读问题上b 参数的定义/基本原理，它永远不会是负数。不过，您的第一行价值观非常合理，我会尝试一下。
这种方法的问题是，通过调整 penalty 我可以让b 返回几乎任何我想要的值，这不好。
@Gabriel 这基本上意味着您的模型过度参数化。也就是说，有太多方法可以很好地拟合数据点，而没有一种方法是真正的拟合。考虑减少模型中的参数数量。

【解决方案3】：

看一眼，好像大号b会淘汰func()的第二个词：

当b/a 趋于无穷时，1 / np.sqrt(1 + (b / a) ** 2)) ** 2 趋于零。

这向我表明，这部分功能在模型中是不需要的，而且弊大于利。

只需将func 设置为：

c * (1 / np.sqrt(1 + (np.asarray(x) / a) ** 2) + d

【讨论】：

不幸的是，这对我不起作用。 b 参数是该函数中最重要的参数。 +1 的答案。谢谢！