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

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

我无法弄清楚如何获得一个可信的,更不用说可靠的,适合这个例程,但我找不到任何其他好的Python曲线拟合例程.我是否需要编写自己的最小二乘算法,或者我在这里做错了什么？

Answer 1

UPDATE

在原帖中,我展示了一个使用的解决方案,lmfit它允许为参数分配边界.从版本0.17开始,scipy还允许直接为参数分配边界(参见文档).请在EDIT之后找到下面的解决方案,这可以作为如何使用scipy curve_fit与参数边界的最小例子.

原帖

正如@Warren Weckesser所建议的那样,你可以使用lmfit来完成这项任务,这允许你为你的参数分配边界并避免这个'丑陋'的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,这也非常接近原始值.

情节如下:你可以看到fit非常好地描述了数据:

以下是包含几个内联注释的完整代码:

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]),分别为(注:字母之一,你将不得不采取的第一个参数的exponantial - 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()