时间戳

首页 > 解决方案 > 对 Python 中不定积分定义的函数使用curve_fit

python - 对 Python 中不定积分定义的函数使用curve_fit

问题描述

我正在尝试编写代码以将具有 5 个参数的 2 条曲线拟合到真实数据。它们显示在这里：

在此处输入图像描述

第一条曲线仅取决于 a、b 和 gamma。所以我决定对这 3 个（有效）使用一次curve_fit，然后在第二条曲线上再次使用它来调整最后两个 alpha 和 k_0。

问题是第二个是由这个不定积分定义的，我无法正确编码。我试图将 x 视为一个符号并使用 sym.integrate 进行积分，然后与 quad 正常积分。都没有奏效。在第二种情况下，我得到“ValueError：具有多个元素的数组的真值不明确。在“mortes”函数中使用 a.any() 或 a.all()。

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy.integrate as integrate
import numpy as np
import sympy as sym


#Experimental x and y data points
#Dados de SP
xData = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34])

ycasos = np.array([2, 13, 65, 459, 1406, 4466, 8419, 13894, 20004, 31174, 44411, 61183, 80558, 107142, 140549, 172875, 215793, 265581, 312530, 366890, 412027, 479481, 552318, 621731, 697530, 749244, 801422, 853085, 890690, 931673, 970888, 1003429, 1034816, 1062634, 1089255])

ymortes = np.array([0, 0, 15, 84, 260, 560, 991, 1667, 2586, 3608, 4688, 6045, 7532, 9058, 10581, 12494, 14263, 15996, 17702, 19647, 21517, 23236, 25016, 26780, 28392, 29944, 31313, 32567, 33927, 35063, 36136, 37223, 37992, 38726, 39311])

#Dados do Brasil    
#xData = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45])  

#ycasos = np.array([2,9,121,1128,3912,10298,20818,36739,58973,96559,155939,233142, 347398, 498440, 672846, 850514, 1067579, 1313667, 1577004, 1839850, 2074860, 2394513, 2707877, 3012412, 3317096, 3582362, 3846153, 4123000, 4315687, 4528240, 4717991, 4906833, 5082637, 5224362, 5380635, 5535605, 5653561, 5848959, 6052786, 6290272, 6577177, 6880127, 7213155, 7465806, 7716405, 8013708])

#ymortes = np.array([])


u0 = ycasos[0]
v0 = ymortes[0]

#u(t)
def casos(x,a,b,gama):
    return u0 * (a ** (1/gama)) * np.exp(a*x) *((a + b * (u0 ** gama) * (np.exp(a*gama*x)-1)) ** (-1/gama))
    

#Plot experimental data points
plt.plot(xData, ycasos, 'bo', label='reais')
 
# Initial guess for the parameters
#initialGuess = [3.0,1.5,0.05]    
 
#Primeiro fit
copt, ccov = curve_fit(casos, xData, ycasos,bounds=(0, [1., 1., np.inf]),maxfev=100000)
a_opt = copt[0]
b_opt = copt[1]
gama_opt = copt[2]
print('Primeira etapa \n')
print('Parametros encontrados: a=%.9f, b=%.9f,gama=%.9f \n' % tuple(copt))


def integrand(t,alpha):
    return np.exp((a_opt - alpha)*t) *((a_opt + b_opt * (u0 ** gama_opt) * (np.exp(a_opt*gama_opt*t)-1)) ** (-1/gama_opt)) 


def mortes(x,k0,alpha):
    return u0 * (a_opt ** (1/gama_opt)) * k0 * integrate.quad(integrand, 0, x, args=(alpha)) + v0

#Segundo fit
mopt, mcov = curve_fit(mortes, xData, ymortes, bounds=(0, [np.inf, a_opt]), maxfev=100000)

print('Segunda etapa \n')
print('Parametros encontrados: k0=%.9f, alpha=%.9f \n' % tuple(mopt))


#x values for the fitted function
xFit = np.arange(0.0, float(len(xData)), 0.01)
 
#Plot the fitted function
plt.plot(xFit, casos(xFit, *copt), 'r', label='estimados')
 
plt.xlabel('t')
plt.ylabel('casos')
plt.legend()
plt.show()

标签： python

解决方案

积分的上限（integrate.quad）必须是浮点数，而不是 x 数组（mortes() 的参数）：

这样它应该可以工作：

def mortes(x,k0,alpha):
    integralRes = []
    for upBound in x: 
        integralRes.append(integrate.quad(integrand, 0, upBound, args=(alpha))[0])
        
    return u0 * (a_opt ** (1/gama_opt)) * k0 * np.array(integralRes) + v0

ps 代码风格的优雅版本非常受欢迎（比如允许将数组传递到Integrate.quad 的上限和下限）。

推荐阅读