Hey @MNGuenther, I also depend on many external functions (sometimes quite complicated like ode solvers etc.) to produce my modeled data. Currently the solution I have been using is to follow the prescription in this post: https://docs.pymc.io/notebooks/blackbox_external_likelihood.html, using a custom likelihood function, since it is in most cases non-trivial to write a theano version of the external functions in my case. I do not make use of the cython aspect of the code currently but following the same principle has worked well for me, though I am sure it is likely not the most efficient method since many of the computations happen outside of theano so likely the entire model is not very optimized (?)
Either way it is a solution, but I am very interested to hear if there is a more appropriate solution.