SciPy曲线拟合失败幂定律 - python

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

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

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

参考方案

更新

在原始帖子中，我展示了一个使用lmfit的解决方案，该解决方案允许为您的参数分配界限。从版本0.17开始，scipy还允许直接为参数分配范围（请参见documentation）。请在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也非常接近原始值。

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

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

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

在返回'Response'(Python)中传递多个参数 - python

我在Angular工作，正在使用Http请求和响应。是否可以在“响应”中发送多个参数。角度文件：this.http.get("api/agent/applicationaware").subscribe((data:any)... python文件：def get(request): ... return Response(seriali…

Python exchangelib在子文件夹中读取邮件 - python

我想从Outlook邮箱的子文件夹中读取邮件。Inbox ├──myfolder 我可以使用account.inbox.all()阅读收件箱，但我想阅读myfolder中的邮件我尝试了此页面folder部分中的内容，但无法正确完成https://pypi.python.org/pypi/exchangelib/ 参考方案 您需要首先掌握Folder的myfo…

python JSON对象必须是str，bytes或bytearray，而不是'dict - python

在Python 3中，要加载以前保存的json，如下所示:json.dumps(dictionary)输出是这样的{"('Hello',)": 6, "('Hi',)": 5}当我使用json.loads({"('Hello',)": 6,…

R'relaimpo'软件包的Python端口 - python

我需要计算Lindeman-Merenda-Gold（LMG）分数，以进行回归分析。我发现R语言的relaimpo包下有该文件。不幸的是，我对R没有任何经验。我检查了互联网，但找不到。这个程序包有python端口吗？如果不存在，是否可以通过python使用该包？ python参考方案 最近，我遇到了pingouin库。

Python ThreadPoolExecutor抑制异常 - python

from concurrent.futures import ThreadPoolExecutor, wait, ALL_COMPLETED def div_zero(x): print('In div_zero') return x / 0 with ThreadPoolExecutor(max_workers=4) as execut…