Using pm.DensityDist and customized likelihood with a black-box function

Yes. There are a few thread on that, see eg Custom theano Op to do numerical integration