Using quad from scipy.integrate with pymc RVs as parameter

Questions
1

I’m trying to write a custom random variable, in which I need to do numerical integration with pymc RVs as parameter. An example is stated as below:

#define some pymc RVs
a1=pm.Normal(‘a1’,mu=0,sd=1)
a2=pm.Nromal(‘a2’,mu=0,sd=1)

#the function to be integrated
def f(x,p):
retrun tt.power(x,p[0]*p[1])

scipy.integrate.quad(f,0,1,args=[a1,a2])

However, quad would return an error since a1,a2 are not real number, and I can’t find something like a1.value() to get the current value of the random variable.

Does anyones knows how to do this integration?

Thanks a lot!

1 Like
2

I don’t think you can use the integration from Scipy directly. You should try to find a closed-form solution to the integration (or an approximation) and rewrite it in theano.

You can have a look at PR #2688, where @domenzain is implementing logcdf for distribution. There are quite a lot of integration involved.

3

Thank you! The actual integration is quite complicated, I’ll try pymc2 for my problem.

Would it be possile to add the .value() attribute to the random varibles in future version? It’s quite helpful at least in my case.

1 Like
4

You can have a look at this post, which has a theano integrateOut op which the use scope is also in PyMC3 https://stackoverflow.com/questions/42678490/custom-theano-op-to-do-numerical-integration/43154053

2 Likes
5

I’ve tried this integrateOut op but it doesn’t really work in my case, although it seems to be helpful in many other cases.
For the example in this integrateOut op:

intZ = integrateOut(z,t,0,1)(x,y)
functionZ = theano.function([x,y],intZ)

I can’t treated intZ, a theano.op, as a pymc3 RV(I’m not sure about this, please tell me if there is some way to do it) and define a new RV like:

M=pm.Normal('M',mu=intZ,sd=intZ*0.1,observed=0.25)

If I want to calculate the logp manually, I still need to know the value of x and y(in order to use the functionZ), which are RVs in my case.

6

You are right, you cannot treat a theano.op as a pymc3 RV. You need to evaluate it. For example, is it work doing?:

M = pm.Normal('M',mu=intZ(),sd=1,,observed=0.25)
7

Sorry, I’ve made a mistake. intZ here is not a theano.op. It’s a tensor variable since it has already been evaluated by (x,y).

In fact, it won’t give an error in the step defining M, but would return errors while sampling the model, if x and y are RVs.

8

pinging @aseyboldt for a more expert take :wink: