Using Pymc3+Theano+astropy for Bayesian inference with integral expressions

This is my first time using Pymc3 or Theano, so I apologize if this question is straightforward. I'm interested in using Bayesian inference to see how effective the (non) observation of something can be at determining the value of model parameters.

I am interested in the probability of a model parameter being equal to a certain value given N observations of a phenomenon. I'd like to do this with Bayesian statistics. So, I would take a Poisson distribution and multiply that by my prior for the model parameter in question.

The problem is that I don't have a rate, I have a complicated expression for the differential rate of the phenomenon with respect to redshift and luminosity. These expressions in turn use integral equations to compute things such as the proper distance between the observer and source. Some of these can be handled by astropy, but I eventually have to numerically integrate if I want a rate in number per second that I can use.

I gather the correct way to do this would be handled by a Theano Op, but I'm having trouble understanding how to correctly structure my code.

I guess this boils down to the question: how would I structure code that computes:

$$p(\mu|N)=p(N|q)p(q)$$

with

$$p(N|q)$$ a Poisson distribution with parameter $$R\times observing~time$$ given that

$$R=\int_0^zdz'\int dL\frac{d^2R}{dz'dL}(q, z, ...)$$

using Theano.

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Nov 4 '21 at 20:56