Using curve_fit to a function defined by indefinite integral in Python

I’m trying to make a code to fit 2 curves with 5 parameters to real data. They are shown here:

The first curve only depends on a,b and gamma. So I decided to use curve_fit once to these 3 (which works) and then use it again on the second curve to adjust the last two alpha and k_0.

Problem is that the second is defined by this indefinite integral and i can’t code it properly. I have tried to treat x as a symbol and integrate using sym.integrate and just integrate normally with quad. Neither worked. In the second case, I get “ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()” in “mortes” function.

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()

Answer

The upper bound of an integral (integrate.quad) has to be a float, not an array as your x (argument of mortes()):

In this way it should work:

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

p.s. Elegant editions of the code style are more than welcomed (like allowing passing an array to upper and lower bounds of integrate.quad ).

Advertisement

Answer