I’m trying to create a custom likelihood function similar to the one described here: P(d|m) \propto exp \frac{-||d-f(m)||^2}{\sigma^2} where again d is observed data and m are parameters.

In my case however, f involves a numerical integration, something like f(\lambda, m) \propto \int{\eta(\lambda) h(\lambda, m) d\lambda}, where \lambda is the numerical domain over which I need to evaluate f. The function h(\lambda, m) is analytically defined, but \eta(\lambda) is some empirically-measured function defined over discrete values of \lambda. Am I wrong in thinking this rules out using the popular solution found here? If that won’t work, is there any way to do this integration?

I am new to the world of Bayesian data analysis, PyMC3, and Theano, but I think if Theano had a tt.trapz() function like numpy’s np.trapz(), I might be able to make this work, using a strategy like this one.

Any advice is greatly appreciated!