带约束和固定点的scipy曲线拟合

deyfvvtc  于 2023-08-05  发布在  其他
关注(0)|答案(1)|浏览(120)

我试图使用SciPy的optimize.curve_fit将函数拟合到一些分散的数据,但我需要拟合曲线下的面积与基于分散的数据计算的面积相同,并且曲线通过数据的起始点和终点。为了做到这一点,我在惩罚公式中使用由散射数据定义的面积(积分),如this answer,同时用参数sigma加权拟合,如建议here
不幸的是,当包含积分约束时,我不能让我的拟合通过初始点和终点。如果忽略积分约束,拟合工作正常并通过点。积分约束和积分约束不可能同时满足吗?我在Windows 10上使用Python 3.7.10。

import scipy
import numpy as np
import matplotlib.pyplot as plt

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

def Func(x,a,b,c):
    return a*x**2 + b*x + c

# modified function definition with penalization
def FuncPen(x,a,b,c):
    integral = scipy.integrate.quad(Func, x[0], x[-1], args=(a,b,c))[0]
    penalization = abs(np.trapz(y,x)-integral)*10000
    return a*x**2 + b*x + c + penalization

sigma = np.ones(len(x))
sigma[[0, -1]] = 0.0001 # first and last points

popt1, _ = scipy.optimize.curve_fit(Func, x, y, sigma=sigma)
popt2, _ = scipy.optimize.curve_fit(FuncPen, x, y, sigma=sigma)

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

fig, ax = plt.subplots(1)
ax.scatter(x,y)
ax.plot(x,y_fit1, color='g', alpha=0.75, label='curve_fit')
ax.plot(x,y_fit2, color='b', alpha=0.75, label='constrained')
plt.legend()

字符串


的数据

fcy6dtqo

fcy6dtqo1#

非常感谢埃尔温·卡尔维拉根对这个问题的开放性评论。我在这里发布我的解决方案:

import scipy
import numpy as np
import matplotlib.pyplot as plt

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

def Func(x,a,b,c):
    return a*x**2 + b*x + c

# modified function definition with penalization
def FuncPen(x,a,b,c):
    integral = scipy.integrate.quad(Func, x[0], x[-1], args=(a,b,c))[0]
    penalization = abs(np.trapz(y,x)-integral)*10000
    return a*x**2 + b*x + c + penalization

# Writing as a general constraint problem
def FuncNew(x,params):
    return params[2]*x**2 + params[1]*x + params[0]

def ConstraintIntegral(params):
    integral = scipy.integrate.quad(FuncNew, x[0], x[-1], args=(params,))[0]
    return integral- np.trapz(y,x)

def ConstraintBegin(params):
    return y[0] - FuncNew(x[0],params)

def ConstraintEnd(params):
    return y[-1] - FuncNew(x[-1],params)

def Objective(params,x,y):
    y_pred = FuncNew(x,params)
    return np.sum((y_pred - y) ** 2) # least squares

cons = [{'type':'eq', 'fun': ConstraintIntegral},
        {'type':'eq', 'fun': ConstraintBegin},
        {'type':'eq', 'fun': ConstraintEnd}]

sigma = np.ones(len(x))
sigma[[0, -1]] = 0.0001 # first and last points

popt1, _ = scipy.optimize.curve_fit(Func, x, y, sigma=sigma)
popt2, _ = scipy.optimize.curve_fit(FuncPen, x, y, sigma=sigma)

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

new = scipy.optimize.minimize(Objective, x0=popt1, args=(x,y), constraints=cons)
popt3 = new.x
y_fit3 = FuncNew(x,popt3)

fig, ax = plt.subplots(1)
ax.scatter(x,y)
ax.plot(x,y_fit1, color='g', alpha=0.75, label='curve_fit')
ax.plot(x,y_fit2, color='b', alpha=0.75, label='constrained')
ax.plot(x,y_fit3, color='r', alpha=0.75, label='generally constrained')
plt.legend()

字符串


的数据

相关问题