python - 如何对 SciPy 曲线拟合施加约束？

Question

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

Answer 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 拟合优度估计器

Answer 2

这是一个几乎相同的片段，它只使用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

Answer 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

Answer 4

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

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

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

Answer 5

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

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)