如何对 SciPy 曲线拟合施加约束？答案

【问题标题】：How do I put a constraint on SciPy curve fit?如何对 SciPy 曲线拟合施加约束？
【发布时间】：2013-05-08 14:46:20
【问题描述】：

我正在尝试使用自定义概率密度函数来拟合一些实验值的分布。显然，结果函数的积分应该始终等于 1，但是简单 scipy.optimize.curve_fit(function, dataBincenters, dataCounts) 的结果永远不会满足这个条件。解决此问题的最佳方法是什么？

【问题讨论】：

标签： python optimization scipy curve-fitting

【解决方案1】：

您可以定义自己的残差函数，包括一个惩罚参数，如下面的代码中详述的那样，其中预先知道沿区间的积分必须是2.。如果你在没有惩罚的情况下进行测试，你会发现你得到的是传统的curve_fit：

import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit, minimize, leastsq
from scipy.integrate import quad
from scipy import pi, sin
x = scipy.linspace(0, pi, 100)
y = scipy.sin(x) + (0. + scipy.rand(len(x))*0.4)
def func1(x, a0, a1, a2, a3):
    return a0 + a1*x + a2*x**2 + a3*x**3

# here you include the penalization factor
def residuals(p,x,y):
    integral = quad( func1, 0, pi, args=(p[0],p[1],p[2],p[3]))[0]
    penalization = abs(2.-integral)*10000
    return y - func1(x, p[0],p[1],p[2],p[3]) - penalization

popt1, pcov1 = curve_fit( func1, x, y )
popt2, pcov2 = leastsq(func=residuals, x0=(1.,1.,1.,1.), args=(x,y))
y_fit1 = func1(x, *popt1)
y_fit2 = func1(x, *popt2)
plt.scatter(x,y, marker='.')
plt.plot(x,y_fit1, color='g', label='curve_fit')
plt.plot(x,y_fit2, color='y', label='constrained')
plt.legend(); plt.xlim(-0.1,3.5); plt.ylim(0,1.4)
print 'Exact   integral:',quad(sin ,0,pi)[0]
print 'Approx integral1:',quad(func1,0,pi,args=(popt1[0],popt1[1],
                                                popt1[2],popt1[3]))[0]
print 'Approx integral2:',quad(func1,0,pi,args=(popt2[0],popt2[1],
                                                popt2[2],popt2[3]))[0]
plt.show()

#Exact   integral: 2.0
#Approx integral1: 2.60068579748
#Approx integral2: 2.00001911981

其他相关问题：

SciPy LeastSq Goodness of Fit Estimator

【讨论】：

@Axon 哪个警告？如果您可以将代码粘贴到网络上的某个位置，那就太好了...
我尝试了这种方法，但在这种情况下，我收到了这个警告：dpaste.org/NfLwy，并且得到的拟合曲线甚至与分布几乎不一样。 scipy.optimize.curve_fit with penalization
@Axon 这是一个集成错误。我在这里检查，但你可以尝试另一个惩罚因素10000 看看会发生什么。您还可以将初始guess=(1.,1.,1.,1.) 更改为另一次尝试
@Axon This answer 或 this answer 可能会让您了解如何拟合分布函数
@altroware 没有特殊原因，但由于 curve_fit 是 leastsq 周围的 Python 包装器，我更喜欢使用后者......但如果有一个新的答案 curve_fit 就好了； )

【解决方案2】：

这是一个几乎相同的 sn-p，它只使用curve_fit。

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt
import scipy.integrate as integr


x = np.linspace(0, np.pi, 100)
y = np.sin(x) + (0. + np.random.rand(len(x))*0.4)

def Func(x, a0, a1, a2, a3):
    return a0 + a1*x + a2*x**2 + a3*x**3

# modified function definition with Penalization
def FuncPen(x, a0, a1, a2, a3):
    integral = integr.quad( Func, 0, np.pi, args=(a0,a1,a2,a3))[0]
    penalization = abs(2.-integral)*10000
    return a0 + a1*x + a2*x**2 + a3*x**3 + penalization


popt1, pcov1 = opt.curve_fit( Func, x, y )
popt2, pcov2 = opt.curve_fit( FuncPen, x, y )

y_fit1 = Func(x, *popt1)
y_fit2 = Func(x, *popt2)

plt.scatter(x,y, marker='.')
plt.plot(x,y_fit2, color='y', label='constrained')
plt.plot(x,y_fit1, color='g', label='curve_fit')
plt.legend(); plt.xlim(-0.1,3.5); plt.ylim(0,1.4)
print 'Exact   integral:',integr.quad(np.sin ,0,np.pi)[0]
print 'Approx integral1:',integr.quad(Func,0,np.pi,args=(popt1[0],popt1[1],
                                                popt1[2],popt1[3]))[0]
print 'Approx integral2:',integr.quad(Func,0,np.pi,args=(popt2[0],popt2[1],
                                                popt2[2],popt2[3]))[0]
plt.show()

#Exact   integral: 2.0
#Approx integral1: 2.66485028754
#Approx integral2: 2.00002116217

【讨论】：

【解决方案3】：

按照上面的示例是添加任何约束的更一般的方法：

from scipy.optimize import minimize
from scipy.integrate import quad
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, np.pi, 100)
y = np.sin(x) + (0. + np.random.rand(len(x))*0.4)

def func_to_fit(x, params):
    return params[0] + params[1] * x + params[2] * x ** 2 + params[3] * x ** 3

def constr_fun(params):
    intgrl, _ = quad(func_to_fit, 0, np.pi, args=(params,))
    return intgrl - 2

def func_to_minimise(params, x, y):
    y_pred = func_to_fit(x, params)
    return np.sum((y_pred - y) ** 2)

# Do the parameter fitting
#without constraints
res1 = minimize(func_to_minimise, x0=np.random.rand(4), args=(x, y))
params1 = res1.x
# with constraints
cons = {'type': 'eq', 'fun': constr_fun}
res2 = minimize(func_to_minimise, x0=np.random.rand(4), args=(x, y), constraints=cons)
params2 = res2.x

y_fit1 = func_to_fit(x, params1)
y_fit2 = func_to_fit(x, params2)

plt.scatter(x,y, marker='.')
plt.plot(x, y_fit2, color='y', label='constrained')
plt.plot(x, y_fit1, color='g', label='curve_fit')
plt.legend(); plt.xlim(-0.1,3.5); plt.ylim(0,1.4)
plt.show()
print(f"Constrant violation: {constr_fun(params1)}")

违反约束：-2.9179325622408214e-10

【讨论】：

【解决方案4】：

如果您能够提前标准化您的概率拟合函数，那么您可以使用此信息来限制您的拟合。一个非常简单的例子是将高斯拟合到数据中。如果要拟合以下三参数 (A, mu, sigma) 高斯，那么它通常是非归一化的：

但是，如果改为对 A 强制执行标准化条件：

那么高斯只有两个参数，会自动归一化。

【讨论】：

【解决方案5】：

您可以确保通过数值积分对拟合的概率分布进行归一化。例如，假设您有数据 x 和 y，并且您已经为概率分布定义了带有参数 a 和 b 的 unnormalised_function(x, a, b)，该概率分布在区间 x1 到 @987654327 上定义@（可能是无限的）：

from scipy.optimize import curve_fit
from scipy.integrate import quad

# Define a numerically normalised function
def normalised_function(x, a, b):
    normalisation, _ = quad(lambda x: unnormalised_function(x, a, b), x1, x2)
    return unnormalised_function(x, a, b)/normalisation

# Do the parameter fitting
fitted_parameters, _ = curve_fit(normalised_function, x, y)

【讨论】：