You can try using the Custom theano Op to do numerical integration. Write down your model in PyMC3 as usual, and then get the model.logpt (which is a tensor of the output log(p(\theta \ | \, x)p(x \ | \,a)p(a))), and use the custom theano numerical integration to integrate re x, and then maximize this new function.
Not sure how straightforward it is though…