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