0

I have been trying to apli this code to pymc3, but I got ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all(). I got the error on evaluating the integral, here is the code. 

import pymc3 as pm
import numpy as np
import scipy.integrate as integrate
data= np.loadtxt('sn_union2.txt',comments="#")
x=data[:,0]
y=data[:,1]
yerr=data[:,2]
c=3e5
mul=[]
def Or(h):
    Neff= 3.04
    Orad= 2.469*10**(-5)*h**(-2)*(1+0.2271* Neff)
    return Orad
def H(z,Om,Orc,h):
    Orad= Or(h)
    H0= 100*h
    E=np.sqrt((1-(np.sqrt(Orc)+np.sqrt(Orc+Om))**2)*(1+z)**2+(np.sqrt(Orc)+np.sqrt(Orc+Om*(1+z)**3))**2)
    return E
def Hinv(z,Om,Orc,h):
    return 1.0/H(z,Om,Orc,h)         
def dl(z,Om,Orc,h):
    int=integrate.quad(Hinv,0,z,args=(Om,Orc,h))[0]
    return int
def mu(z,Om,Orc,h, mu0):
    H0=100.0*h
    smu=5.0*np.log10((c/H0)*(1.0+z)*dl(z,Om,Orc,h)) + mu0
    return smu
with pm.Model() as model:
    Orc=pm.Uniform('Orc', lower=0, upper=0.25, transform=None)
    Om=pm.Uniform('Om', lower=0, upper=0.5, transform=None)
    mu0= pm.Uniform('mu0', lower=0, upper=100,   transform=None)
    h=pm.Normal('h',mu=0.7324, sigma=0.0174)
    mu=mu(x,Om,Orc,h,mu0)
    likelihood=pm.Normal('Y', mu=mu,sd=yerr, observed=y)
with model:
    step = pm.Metropolis ()
    trace = pm.sample (2000, tune=100, init=None, step=step, cores=8)
x=trace.get_values('Orc',burn=100, combine=True)
y=trace.get_values('Om',burn=100, combine=True)
z=trace.get_values('h',burn=100, combine=True)
mu0=trace.get_values('mu0',burn=100, combine=True)
t=list(zip(x,y,z,mu0))
np.savetxt('SNIaDGP.txt',t)
JohanC
  • 71,591
  • 8
  • 33
  • 66
Monica Castillo
  • 1
  • 1
    Normally this happens when using the "or" or "and" keywords between numpy arrays. I don't see this here. Can you include the full error message and traceback in your question body? – ijustlovemath Apr 13 '20 at 22:14
  • Here's the solution to the typical problem: https://stackoverflow.com/questions/10062954/valueerror-the-truth-value-of-an-array-with-more-than-one-element-is-ambiguous – ijustlovemath Apr 13 '20 at 22:14
  • 2
    A more prohibitive issue is that throwing a numerical integration over a random variable doesn't really work like this code seems to expect it to. PyMC3 random variables are not numerical values (e.g., numpy ndarrays) but instead nodes in a Theano computational graph, i.e., non-trivial data structures. If you have a complex operation in the sampling, you'll likely need to [write a custom Theano Op](https://docs.pymc.io/Advanced_usage_of_Theano_in_PyMC3.html#writing-custom-theano-ops). And [the PyMC Discourse](https://discourse.pymc.io) is likely a better support forum if you go this route. – merv Apr 14 '20 at 02:40

0 Answers0