SciPy 曲线拟合不符合幂律

Question

因此，我正在尝试使用以下类型的幂律拟合一组数据：

def f(x,N,a): # Power law fit
    if a >0:
        return N*x**(-a)
    else:
        return 10.**300

par,cov = scipy.optimize.curve_fit(f,data,time,array([10**(-7),1.2]))

其中 else 条件只是强制 a 为正。使用 scipy.optimize.curve_fit 产生 an awful fit (green line)，N 和 a 的返回值分别为 1.2e+04 和 1.9e0-7，与数据绝对没有交集。从我手动输入的拟合来看，N 和 a 的值应该分别在 1e-07 和 1.2 左右，尽管将它们放入 curve_fit 作为初始参数不会改变结果。删除 a 为正的条件会导致更差的拟合，因为它选择负，这会导致拟合错误的符号斜率。

我不知道如何从这个例程中得到可信的，更不用说可靠的拟合了，但我找不到任何其他好的 Python 曲线拟合例程。我需要编写自己的最小二乘算法还是我做错了什么？

Answer 1

更新

在最初的 post 中，我展示了一个使用 lmfit 的解决方案，它允许为您的参数分配边界。从 0.17 版开始，scipy 还允许直接为您的参数分配边界（参见 documentation）。请在 EDIT 之后找到下面的解决方案，它有望作为一个最小示例，说明如何将 scipy 的 curve_fit 与参数边界一起使用。

原创post

正如@Warren Weckesser 所建议的，您可以使用 lmfit 来完成此任务，这允许您为参数分配界限并避免此 'ugly' if 子句。

由于您没有提供任何数据，我创建了一些数据，如下所示：

他们遵守法律f(x) = 10.5 * x ** (-0.08)

我适合它们 - 正如@roadrunner66 所建议的 - 通过将幂律转换为线性函数：

y = N * x ** a
ln(y) = ln(N * x ** a)
ln(y) = a * ln(x) + ln(N)

所以我先用np.log对原始数据做拟合。当我现在使用 lmfit 时，我得到以下输出：

[[Variables]]
    lN:   2.35450302 +/- 0.019531 (0.83%) (init= 1.704748)
    a:   -0.08035342 +/- 0.005158 (6.42%) (init=-0.5)

所以 a 非常接近原始值，np.exp(2.35450302) 给出 10.53 也非常接近原始值。

剧情如下；如您所见，拟合很好地描述了数据：

这是带有一些内联注释的完整代码：

import numpy as np
import matplotlib.pyplot as plt
from lmfit import minimize, Parameters, Parameter, report_fit

# generate some data with noise
xData = np.linspace(0.01, 100., 50.)
aOrg = 0.08
Norg = 10.5
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData))
plt.plot(xData, yData, 'bo')
plt.show()

# transform data so that we can use a linear fit
lx = np.log(xData)
ly = np.log(yData)
plt.plot(lx, ly, 'bo')
plt.show()

def decay(params, x, data):

    lN = params['lN'].value
    a = params['a'].value

    # our linear model
    model = a * x + lN
    return model - data # that's what you want to minimize

# create a set of Parameters
params = Parameters()
params.add('lN', value=np.log(5.5), min=0.01, max=100)  # value is the initial value
params.add('a', value=-0.5, min=-1, max=-0.001)  # min, max define parameter bounds

# do fit, here with leastsq model
result = minimize(decay, params, args=(lx, ly))

# write error report
report_fit(params)

# plot data
xnew = np.linspace(0., 100., 5000.)
# plot the data
plt.plot(xData, yData, 'bo')
plt.plot(xnew, np.exp(result.values['lN']) * xnew ** (result.values['a']), 'r')
plt.show()

编辑

假设您安装了 scipy 0.17，您还可以使用 curve_fit 执行以下操作。我展示了您对幂律的原始定义（下图中的红线）以及对数数据（下图中的黑线）。数据的生成方式与上述相同。情节如下：

如您所见，数据描述的非常好。如果您打印 popt 和 popt_log，您将分别获得 array([ 10.47463426, 0.07914812]) 和 array([ 2.35158653, -0.08045776])（注意：对于字母 1，您必须取第一个参数的指数 - np.exp(popt_log[0]) = 10.502接近原始数据）。

完整代码如下：

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

# generate some data with noise
xData = np.linspace(0.01, 100., 50)
aOrg = 0.08
Norg = 10.5
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData))

# get logarithmic data
lx = np.log(xData)
ly = np.log(yData)

def f(x, N, a):
    return N * x ** (-a)

def f_log(x, lN, a):
    return a * x + lN

# optimize using the appropriate bounds
popt, pcov = curve_fit(f, xData, yData, bounds=(0, [30., 20.]))
popt_log, pcov_log = curve_fit(f_log, lx, ly, bounds=([0, -10], [30., 20.]))

xnew = np.linspace(0.01, 100., 5000)

# plot the data
plt.plot(xData, yData, 'bo')
plt.plot(xnew, f(xnew, *popt), 'r')
plt.plot(xnew, f(xnew, np.exp(popt_log[0]), -popt_log[1]), 'k')
plt.show()

SciPy 曲线拟合不符合幂律

SciPy Curve Fit Fails Power Law

python

numpy

optimization

scipy

curve-fitting