I am trying to draw samples from a posterior distribution of random parameters of an experiment simulation. I have a model that takes the random parameters as an input, solves a PDE, and returns the value. Then, my likelihood is
p(Y | theta) = C * Exp(-0.5 * (g(theta) - Y)’ Sigma^-1 (g(theta) - Y))
theta can be either Gaussian or Uniform distributed, g is my model, and is the experiment observation.
C is the constant (2*pi)^d/2 |det(Sigma)|^-1/2
For simple toy model g, say
g(theta) = 2*theta -1
I could use pm.Normal like:
theta = pm.Uniform(‘theta’, lower=-3, upper=3)
Y_obs = pm.Normal(‘Y_obs’, mu=model(theta), sigma=sigma2_e**.5, observed=Y)
and Theano would just deal with the model.
However, when I try to use the real model, with the code for the solution of the PDE, it obviously crashes.
Is there a way to define my likelihood? I tried reading the pymc3 tutorials and examples and could not get any idea on how to solve my problem. Any help would be appreciated.