Integrate function of random variables in pymc3

Question

I am trying to construct a model in pymc3 which requires me to integrate a function of random variables. The basic idea using actual numbers is this:

from scipy.integrate import quad

def PowerLaw(x,N0,alpha):
    """ A PowerLaw distribution"""
    return N0 * x**-alpha

print quad(PowerLaw,0.1,10,args=(1e-4,2)) #(0.0009900000000000002, 3.008094474110274e-12)

I can also do this in theano:

from theano import function,tensor as tt

xt = tt.dscalar('x')
N0 = tt.dscalar('N0')
alpha = tt.dscalar('alpha')
y = PowerLaw(xt,N0,alpha)
func = function([xt,N0,alpha],y)
print quad(func,0.1,10,args=(1e-4,2)) #Same answer as before

Here is an example of what I want to do:

with pm.Model() as myModel:
    N0 = pm.Uniform("N0",1e-5,1e-1)
    alpha = pm.Uniform("alpha",1,5)
    yval = quad(PowerLaw,0.1,10,args=(N0,alpha))

But of course when I try this I get a TypeError because N0 and alpha are not floats. Of course, in this simple case, I know the analytical solution of the integral; my actual model requires more complicated integrals where I do not know the closed form. Is there any way to accomplish this in pymc3?

Integrate function of random variables in pymc3

0 Answers0