Python에서 무한 적분으로 정의 된 함수에 curve_fit 사용

실제 데이터에 5 개의 매개 변수를 사용하여 2 개의 곡선을 맞추는 코드를 만들려고합니다. https://i.stack.imgur.com/Yrydc.png 여기에 표시됩니다.

첫 번째 곡선은 a, b 및 감마에만 의존합니다. 그래서 나는 curve_fit 을이 3 (작동하는)에 한 번 사용하고 두 번째 곡선에서 다시 사용하여 마지막 두 알파와 k_0을 조정하기로 결정했습니다.

문제는 두 번째가이 무한 적분으로 정의되어 제대로 코딩 할 수 없다는 것입니다. 나는 x를 기호로 취급하고 sym.integrate를 사용하여 통합하고 normaly를 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()