SciPy曲线拟合失败幂定律 [英] SciPy Curve Fit Fails Power Law

问题描述

因此，我正在尝试使用以下幂定律拟合一组数据:

So, I'm trying to fit a set of data with a power law of the following kind:

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会产生可怕的拟合(绿线)，返回值1.2e + N和a分别为04和1.9e0-7，与数据绝对没有交集.根据我手动输入的拟合，N和a的值应分别落在1e-07和1.2左右，尽管将它们放入curve_fit中是因为初始参数不会改变结果.删除a为正数的条件会导致拟合度变差，因为它选择了负数，这会导致符号斜率的拟合错误.

where the else condition is just to force a to be positive. Using scipy.optimize.curve_fit yields an awful fit (green line), returning values of 1.2e+04 and 1.9e0-7 for N and a, respectively, with absolutely no intersection with the data. From fits I've put in manually, the values should land around 1e-07 and 1.2 for N and a, respectively, though putting those into curve_fit as initial parameters doesn't change the result. Removing the condition for a to be positive results in a worse fit, as it chooses a negative, which leads to a fit with the wrong sign slope.

我无法弄清楚如何从该例程中获得可信的，更不用说可靠的参数，但是我找不到任何其他良好的Python曲线拟合例程.我需要编写自己的最小二乘算法还是在这里做错了什么?

I can't figure out how to get a believable, let alone reliable, fit out of this routine, but I can't find any other good Python curve fitting routines. Do I need to write my own least-squares algorithm or is there something I'm doing wrong here?

推荐答案

更新

在原始帖子中，我展示了一个使用lmfit的解决方案，该解决方案允许为您的参数分配界限.从版本0.17开始，scipy还允许直接为参数分配范围(请参见

In the original post, I showed a solution that uses lmfit which allows to assign bounds to your parameters. Starting with version 0.17, scipy also allows to assign bounds to your parameters directly (see documentation). Please find this solution below after the EDIT which can hopefully serve as a minimal example on how to use scipy's curve_fit with parameter bounds.

原始帖子

如@Warren Weckesser所建议，您可以使用 lmfit 为了完成此任务，这使您可以为参数分配范围，并避免出现这种难看的" if子句.

As suggested by @Warren Weckesser, you could use lmfit to get this task done, which allows you to assign bounds to your parameters and avoids this 'ugly' if-clause.

由于您未提供任何数据，因此我创建了一些数据，如下所示:

Since you do not provide any data, I created some which are shown here:

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

按照@ roadrunner66的建议，我通过线性函数变换幂定律来拟合它们:

I fit them - as suggested by @roadrunner66 - by transforming the power law in a linear function:

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

因此，我首先在原始数据上使用np.log，然后进行拟合.现在使用lmfit时，将获得以下输出:

So I first use np.log on the original data and then do the fit. When I now use lmfit, I get the following output:

[[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也非常接近原始值.

So a is pretty close to the original value and np.exp(2.35450302) gives 10.53 which is also very close to the original value.

该图如下所示；如您所见，拟合度很好地描述了数据:

The plot then looks as follows; as you can see the fit describes the data very well:

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

Here is the entire code with a couple of inline comments:

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执行以下操作.我将其显示为您对幂律的原始定义(下图中的红线)以及对数数据(下图中的黑线).数据的生成方法与上述相同.该图的外观如下:

Assuming that you have scipy 0.17 installed, you can also do the following using curve_fit. I show it for your original definition of the power law (red line in the plot below) as well as for the logarithmic data (black line in the plot below). The data is generated in the same way as above. The plot the looks as follows:

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

As you can see, the data is described very well. If you print popt and popt_log, you obtain array([ 10.47463426, 0.07914812]) and array([ 2.35158653, -0.08045776]), respectively (note: for the letter one you will have to take the exponantial of the first argument - np.exp(popt_log[0]) = 10.502 which is close to the original data).

这是完整的代码:

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曲线拟合失败幂定律的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！