使用 scipy curve_fit 不合适

Question

我正在尝试将一些数据拟合到具有指数截断的幂律函数。我用 numpy 生成了一些数据，我试图用 scipy.optimization 来拟合这些数据。这是我的代码：

import numpy as np
from scipy.optimize import curve_fit

def func(x, A, B, alpha):
    return A * x**alpha * np.exp(B * x)

xdata = np.linspace(1, 10**8, 1000)
ydata = func(xdata, 0.004, -2*10**-8, -0.75)
popt, pcov = curve_fit(func, xdata, ydata)
print popt

我得到的结果是：[1, 1, 1] 与数据不符。 ¿我做错了什么吗？

Answer 1

显然您的初始猜测（默认为 [1,1,1]，因为您没有给出 - 参见 the docs）与实际参数相差太大，无法使算法收敛。主要问题可能出在 B 上，如果它是正数，它会将您的指数函数发送到您提供的 xdata.

的非常大的值

尝试提供一些更接近实际参数的东西，它起作用了：

p0 = 0.01, -5e-7, -0.4    # Initial guess for the parameters
popt, pcov = curve_fit(func, xdata, ydata, p0)
print popt

输出：

[  4.00000000e-03  -2.00000000e-08  -7.50000000e-01]

Answer 2

虽然 xnx 给了你为什么 curve_fit 在这里失败的答案，但我想我会建议一种不同的方法来解决不依赖于梯度下降的函数形式的拟合问题（和因此初步猜测是合理的）

请注意，如果您对正在拟合的函数取对数，您将得到以下形式

每个未知参数（log A、alpha、B）都是线性的

因此，我们可以使用线性代数的机制来解决这个问题，将方程写成矩阵形式

log y = M p

其中 log y 是 y 数据点对数的列向量，p 是未知参数的列向量，M 是矩阵 [[1], [log x], [x]]

或明确

然后可以使用 np.linalg.lstsq

高效地找到最合适的参数向量

您的代码示例问题可以写成

import numpy as np

def func(x, A, B, alpha):
    return A * x**alpha * np.exp(B * x)

A_true = 0.004
alpha_true = -0.75
B_true = -2*10**-8

xdata = np.linspace(1, 10**8, 1000)
ydata = func(xdata, A_true, B_true, alpha_true)

M = np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T

logA, alpha, B = np.linalg.lstsq(M, np.log(ydata))[0]

print "A =", np.exp(logA)
print "alpha =", alpha
print "B =", B

它很好地恢复了初始参数：

A = 0.00400000003736
alpha = -0.750000000928
B = -1.9999999934e-08

另请注意，对于手头的问题，此方法比使用 curve_fit 快 20 倍左右

In [8]: %timeit np.linalg.lstsq(np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T, np.log(ydata))
10000 loops, best of 3: 169 µs per loop


In [2]: %timeit curve_fit(func, xdata, ydata, [0.01, -5e-7, -0.4])
100 loops, best of 3: 4.44 ms per loop

使用 scipy curve_fit 不合适

Wrong fit using scipy curve_fit

python

numpy

curve-fitting

scipy