Search code examples
pythoncurve-fittingintegral

Fitting an equation with an integral, where upper limit is a variable

I'm trying to fit an equation that has an integral. The upper limit of the integral is a variable, and as such I'm having trouble using quad as it only accepts floats as limits. I've tried to get around this using a for loop but I still keep getting the same error of 'ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()'.

The equation I'm trying to fit is:

Fit Equation

I know the constants M and c, and have data for m(z) and z. I'm trying to fit ΩM and ΩΛ. This is my code: 

import numpy as np
import scipy.integrate as integrate
from scipy.optimize import curve_fit

z=[0.03, 0.05, 0.026, 0.075, 0.026, 0.014, 0.101, 0.02, 0.036, 0.045, 0.043, 0.018, 0.079, 
0.088, 0.063, 0.071, 0.052, 0.05]
m=[16.26, 17.63, 16.08, 18.43, 16.28, 14.47, 19.16, 15.18, 16.66, 17.61, 17.19, 15.61, 
18.27, 19.28, 18.24, 18.33, 17.54, 17.69]
c=299792.458 #speed of light
M=-18.316469239 

def integrand(Z,OM,OV):
    return np.power((((1+Z)**2)*(1+OM*Z)-Z*(2+Z)*OV),-0.5)

def curve(z,OM,OV):
    for i in z:
        I=integrate.quad(integrand,0,i,args=(OM,OV))[0]
    return M+5*np.log10(c*(1+z))*I)                       

popts, pcov=curve_fit(curve,z,m)

Thanks in advance and hope that includes everything!

Solution

  • In addition to what was commented you only take the last result of the integration loop into account. The corrected code should be like this:

    from matplotlib import pylab as plt

import numpy as np
import scipy.integrate as integrate
from scipy.optimize import curve_fit

z=np.asarray([0.03, 0.05, 0.026, 0.075, 0.026, 0.014, 0.101, 0.02, 0.036, 0.045, 0.043, 0.018, 0.079, 
0.088, 0.063, 0.071, 0.052, 0.05])
m=np.asarray([16.26, 17.63, 16.08, 18.43, 16.28, 14.47, 19.16, 15.18, 16.66, 17.61, 17.19, 15.61, 
18.27, 19.28, 18.24, 18.33, 17.54, 17.69])

# sorting of the data is not necessary but nicer when plotting
idx = np.argsort(z)
z = np.take(z, idx)
m = np.take(m, idx)

c=299792458 #speed of light
M=-18.316469239 

def integrand(Z,OM,OV):
    return np.power((((1+Z)**2)*(1+OM*Z)-Z*(2+Z)*OV),-0.5)

def curve(z,OM,OV):
    I = [integrate.quad(integrand,0,i,args=(OM,OV))[0] for i in z]
    return M+5*np.log10(c*(1+z)*I)

popts, pcov=curve_fit(curve,z,m, p0=(1,1))

plt.plot(z,m,'o')
plt.plot(z,curve(z,*popts))
plt.show()

    And indeed it gives a nice match:

    data with fit