Truth value of an array with more than one element is ambiguous when using scipy optimize

Question

I have been looking though similar questions, but I can't understand what is going on. My code is the following, I want to minimize a function:

def hazard(Tx,a,b,c):
    Abac=lambda x:(a/b)*(e**(b*x)-1)*e**(-c*(x**2))
    A=(a/b)*(e**(b*Tx)-1)*e**(-c*(Tx**2))
    i=integrate.quad(Abac,0,Tx)
    H=A/(1-i[0])
    return(H)



ages_population=[11.57,10.94,10.11,9.87,10.05,10.51,9.98,8.39,6.79,4.47,2.75] #male
ages_cases_2018=[40,73,123,214,381,447,542,586,567,455,373]#male
ages=[30,35,40,45,50,55,60,65,70,75,80]
incidence=[ages_cases_2018[i]/(ages_population[i]*0.2*0.094*1e6) for i in range(len(ages_population))]
one,two=optimize.curve_fit(hazard,ages,incidence,p0=[1.78e-8,0.204,0.994e-3])

I get the following error:

File "C:\Users\nicol\.spyder-py3\temp.py", line 44, in <module>
    one,two=optimize.curve_fit(hazard,ages,incidence,p0=[1.78e-8,0.204,0.994e-3])

  File "C:\Users\nicol\AppData\Local\Programs\Spyder\pkgs\scipy\optimize\minpack.py", line 789, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)

  File "C:\Users\nicol\AppData\Local\Programs\Spyder\pkgs\scipy\optimize\minpack.py", line 410, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)

  File "C:\Users\nicol\AppData\Local\Programs\Spyder\pkgs\scipy\optimize\minpack.py", line 24, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))

  File "C:\Users\nicol\AppData\Local\Programs\Spyder\pkgs\scipy\optimize\minpack.py", line 485, in func_wrapped
    return func(xdata, *params) - ydata

  File "C:\Users\nicol\.spyder-py3\temp.py", line 23, in hazard
    i=integrate.quad(Abac,0,Tx)

  File "C:\Users\nicol\AppData\Local\Programs\Spyder\pkgs\scipy\integrate\quadpack.py", line 348, in quad
    flip, a, b = b < a, min(a, b), max(a, b)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I don't understand what is wrong, and it's driving me crazy. I know that its related to the "integrate" part of the hazard function because whenever I do something else, it works. But that is the part that is somehow malfunctioning, and I don't know what to do

Answer 1

Your function with scalar inputs:

In [8]: hazard(1,1,1,1)
Out[8]: 1.1243074314863974

giving an array for Tx produces your error:

In [9]: hazard(np.array([1,2,3]),1,1,1)
Traceback (most recent call last):
  Input In [9] in <cell line: 1>
    hazard(np.array([1,2,3]),1,1,1)
  Input In [2] in hazard
    i=integrate.quad(Abac,0,Tx)
  File /usr/local/lib/python3.8/dist-packages/scipy/integrate/_quadpack_py.py:348 in quad
    flip, a, b = b < a, min(a, b), max(a, b)
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

If I add:

def hazard(Tx,a,b,c):
    print('Tx',Tx)
    ...

And run your optimize:

In [13]: one,two=optimize.curve_fit(hazard,ages,incidence,p0=[1.78e-8,0.204,0
    ...: .994e-3])
Tx [30. 35. 40. 45. 50. 55. 60. 65. 70. 75. 80.]
Traceback (most recent call last):
...

curve_fit is passing an array as the Tx. That's the xdata, ages list.

ydata is incidences, a matching list. That's constructed with a list comprehension - it could probably be done with numpy whole-array expression, but it's only done once, so that's not a big deal.

hazard needs to return an array (or list) that matches Tx in size.

Do a straightforward list iteration in hazard:

def fhazard(Tx,a,b,c):
    H = []
    Abac=lambda x:(a/b)*(e**(b*x)-1)*e**(-c*(x**2))
    for t in Tx:
        A=(a/b)*(e**(b*t)-1)*e**(-c*(t**2))
        i=integrate.quad(Abac,0,t)
        h=A/(1-i[0])
        H.append(h)
    return(H)

In [25]: fhazard([1,2,3],1,1,1)
Out[25]: [1.1243074314863974, 0.614543471762487, 0.015056268137507168]

In [27]: one,two=optimize.curve_fit(fhazard,ages,incidence,p0=[1.78e-8,0.204,
    ...: 0.994e-3])
In [28]: one
Out[28]: array([3.00698561e-06, 1.01042613e-01, 4.31148110e-04])
In [29]: two
Out[29]: 
array([[ 3.70806221e-12, -5.25556908e-08, -3.85724391e-10],
       [-5.25556908e-08,  7.62947960e-04,  5.67783901e-06],
       [-3.85724391e-10,  5.67783901e-06,  4.25991424e-08]])

Testing the fit:

In [33]: plt.plot(ages,incidence, ages, fhazard(ages,*one))