Python曲线拟合库,允许我为参数分配边界

Question

我希望能够执行允许我将任意曲线函数拟合到数据的拟合,并允许我在参数上设置任意边界,例如我想拟合函数:

f(x) = a1(x-a2)^a3\cdot\exp(-\a4*x^a5)

并说:

• a2 在以下范围内: (-1, 1)
• a3并且a5是积极的

有一个很好的scipy curve_fit 函数,但它不允许指定参数边界.还有很好的http://code.google.com/p/pyminuit/库可以进行通用最小化,它允许设置参数的边界,但在我的情况下,它没有覆盖.

Answer 1

注意:SciPy版本0.17中的新功能

假设您想要将模型拟合到如下所示的数据:

y=a*t**alpha+b

以及对alpha的约束

0<alpha<2

而其他参数a和b仍然是免费的.然后我们应该以下列方式使用curve_fit的bounds选项:

import numpy as np
from scipy.optimize import curve_fit
def func(t, a,alpha,b):
     return a*t**alpha+b
param_bounds=([-np.inf,0,-np.inf],[np.inf,2,np.inf])
popt, pcov = curve_fit(func, xdata, ydata,bounds=param_bounds)

来源就在这里.

Answer 2

正如已经提到的罗布法尔克，你可以使用，例如，SciPy的非线性优化程序在scipy.minimize减少的任意错误的功能，例如均方误差。

请注意，您提供的函数不一定具有实际值 - 也许这就是您在 pyminuit 中的最小化没有收敛的原因。您必须更明确地对待这一点，请参见示例 2。

下面的例子都使用了L-BFGS-B支持有界参数区域的最小化方法。我把这个答案分成两部分：

1. 具有真实codomain的函数，类似于您提供的函数。我添加了绝对值以确保您提供的函数返回域中的实数 [-3,3)
2. 你给出的实际函数，它有一个复杂的codomain

1.真正的codomain

下面的示例显示了对您函数的这个稍微修改过的版本的优化。

import numpy as np
import pylab as pl
from scipy.optimize import minimize

points = 500
xlim = 3.

def f(x,*p):
    a1,a2,a3,a4,a5 = p
    return a1*np.abs(x-a2)**a3 * np.exp(-a4 * np.abs(x)**a5)

# generate noisy data with known coefficients
p0 = [1.4,-.8,1.1,1.2,2.2]
x = (np.random.rand(points) * 2. - 1.) * xlim
x.sort()
y = f(x,*p0)
y_noise = y + np.random.randn(points) * .05

# mean squared error wrt. noisy data as a function of the parameters
err = lambda p: np.mean((f(x,*p)-y_noise)**2)

# bounded optimization using scipy.minimize
p_init = [1.,-1.,.5,.5,2.]
p_opt = minimize(
    err, # minimize wrt to the noisy data
    p_init, 
    bounds=[(None,None),(-1,1),(None,None),(0,None),(None,None)], # set the bounds
    method="L-BFGS-B" # this method supports bounds
).x

# plot everything
pl.scatter(x, y_noise, alpha=.2, label="f + noise")
pl.plot(x, y, c='#000000', lw=2., label="f")
pl.plot(x, f(x,*p_opt) ,'--', c='r', lw=2., label="fitted f")

pl.xlabel("x")
pl.ylabel("f(x)")
pl.legend(loc="best")
pl.xlim([-xlim*1.01,xlim*1.01])

pl.show()

2. 复杂codomain的扩展

可以通过显式转换为复数并调整误差函数来将上述最小化扩展到复数域：

首先，您将值 x 显式转换为复数值，以确保 f 返回复数值并且可以实际计算负数的小数指数。其次，我们计算实部和虚部的一些误差函数 - 一个直接的候选者是平方复绝对值的均值。

import numpy as np
import pylab as pl
from scipy.optimize import minimize

points = 500
xlim = 3.

def f(x,*p):
    a1,a2,a3,a4,a5 = p
    x = x.astype(complex) # cast x explicitly to complex, to ensure complex valued f
    return a1*(x-a2)**a3 * np.exp(-a4 * x**a5)

# generate noisy data with known coefficients
p0 = [1.4,-.8,1.1,1.2,2.2]
x = (np.random.rand(points) * 2. - 1.) * xlim
x.sort()
y = f(x,*p0)
y_noise = y + np.random.randn(points) * .05 + np.random.randn(points) * 1j*.05

# error function chosen as mean of squared absolutes
err = lambda p: np.mean(np.abs(f(x,*p)-y_noise)**2)

# bounded optimization using scipy.minimize
p_init = [1.,-1.,.5,.5,2.]
p_opt = minimize(
    err, # minimize wrt to the noisy data
    p_init, 
    bounds=[(None,None),(-1,1),(None,None),(0,None),(None,None)], # set the bounds
    method="L-BFGS-B" # this method supports bounds
).x

# plot everything
pl.scatter(x, np.real(y_noise), c='b',alpha=.2, label="re(f) + noise")
pl.scatter(x, np.imag(y_noise), c='r',alpha=.2, label="im(f) + noise")

pl.plot(x, np.real(y), c='b', lw=1., label="re(f)")
pl.plot(x, np.imag(y), c='r', lw=1., label="im(f)")

pl.plot(x, np.real(f(x,*p_opt)) ,'--', c='b', lw=2.5, label="fitted re(f)")
pl.plot(x, np.imag(f(x,*p_opt)) ,'--', c='r', lw=2.5, label="fitted im(f)")

pl.xlabel("x")
pl.ylabel("f(x)")

pl.legend(loc="best")
pl.xlim([-xlim*1.01,xlim*1.01])

pl.show()

笔记

似乎最小化器可能对初始值有点敏感 - 因此我将我的第一个猜测 (p_init) 放在离最佳值不远的地方。如果您不得不与此作斗争，除了全局优化循环之外，您还可以使用相同的最小化程序，例如盆地跳跃或brute。