How do I put a constraint on SciPy curve fit?

ctzwtxfj  于 2022-11-10  发布在  其他
关注(0)|答案(5)|浏览(117)

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

ecr0jaav

ecr0jaav1#

您可以定义自己的残差函数,包括惩罚参数,如下面的代码所示,其中预先知道沿区间的积分必须为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

其他相关问题:

fcipmucu

fcipmucu2#

下面是一个几乎相同的代码片段,它只使用了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

cczfrluj

cczfrluj3#

下面的示例是添加任何约束的更通用的方法:

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

wlwcrazw

wlwcrazw4#

如果你能预先对概率拟合函数进行归一化处理,那么你就可以使用这个信息来约束拟合。一个非常简单的例子就是对数据进行高斯拟合。如果要拟合以下三参数(A,mu,sigma)高斯函数,那么它通常是未归一化的:

然而,如果改为对A:

则高斯仅为两个参数并且被自动归一化。

2ul0zpep

2ul0zpep5#

您可以确保拟合的概率分布是通过数值积分进行归一化的。例如,假设您有数据xy,并且您已经为概率分布定义了一个带有参数abunnormalised_function(x, a, b),该概率分布定义在x1x2的区间上(该区间可以是无穷大):

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)

相关问题