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 值和它的误差会这么大呢？

Answer 1

从快速观察来看，似乎一个很大的 b 会让 func() 中的第二项失效：

当 b/a 变得非常大时，1 / np.sqrt(1 + (b / a) ** 2)) ** 2 会变得接近于零。

这让我觉得这个函数的这一部分在模型中并不需要，反而会带来更多问题。

只需将 func 设置为：

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

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]

Answer 3

我不知道这是不是故意的还是个错误，但在我看来，'b' 和 'a' 以及 'd' 之间的关系很强，而 'b' 和独立变量 '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/）很有用。它允许你对参数设置范围，轻松切换拟合算法，还可以更深入地探索参数的置信区间。