Scipy的curve_fit函数没有给出合理的结果

Question

我有一组简单的 x,y 数据需要拟合，至少乍一看是这样。问题在于 scipy.optimize.curve_fit 返回了一个非常大的拟合参数值，我不知道这是否正确或者我拟合数据的方法是否有问题。

下面的图显示了数据点和最佳拟合曲线（蓝色）。使用的曲线函数（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值和它的误差如此之大？

- Gabriel

---

很抱歉，问题可能出在您拟合的数据（或函数）上。首先，您的误差中有一个错别字（应该是f_err [1] [1]）。我已经使用全局求解器检查了您的示例，并在其他参数中得到了相同的结果，但我甚至得到了b = 3.56564242e + 18。但这并不奇怪，在x = 300点附近的导数非常缓慢 - 它接近最小值，但非常缓慢。您可以手动检查 - 计算关于数据的参数的偏导数，并尝试解决4个非线性方程组的系统。 - Martin

@Martin 感谢指出笔误，我已经修正了（不过数值本身是正确的）。 - Gabriel

我知道。我已经检查过了。 - Martin

Answer 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/）很有用。它允许您对参数设置边界，轻松地在拟合算法之间切换，并且还可以对参数的置信区间进行更加粗略的探索。

- user3059642

---

感谢您提供 lmfit 的提示，Matt。不幸的是，我尝试了这里显示的所有最小化方法 http://lmfit.github.io/lmfit-py/fitting.html#fit-engines-label，结果都一样：参数 b 是巨大的。即使通过最大值限制它也只会返回该最大值，所以没有用。再次感谢您的推荐！ - Gabriel

你看到变量a和b之间有什么相关性吗？如果所有的拟合方法都给出相同的结果，并且限制值会导致极限总是被触发，那么这不是说明问题不明确吗？再次强调，你基本上有c*a/x - a/b + d。除非'a'和'd'已知非常准确，否则'b'的不确定性将会非常大...这就是你所看到的。 - user3059642

Answer 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]

- behzad.nouri

---

在第二行中，您是否得到了b参数的负值？这绝对是不正确的，请阅读问题中b参数的定义/原理，它永远不能为负数。尽管如此，您的第一行数值非常合理，我会尝试一下。 - Gabriel

这种方法的问题在于，通过调整“惩罚”参数，我可以使“b”返回几乎任何我想要的值，这是不好的。 - Gabriel

@Gabriel，这基本上意味着你的模型过度参数化了。也就是说，有太多方法可以使数据点拟合良好，而没有一个真正的拟合。考虑减少模型中的参数数量。 - behzad.nouri

Answer 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

- Bitwise

---

很遗憾，这对我不起作用。在该函数中，参数b是最重要的一个。对于这个答案加一分。谢谢！ - Gabriel