Hierarchical modeling with a black box likelihood function

Questions
3

I think you’re looking for DensityDist. As long as you have an implementation of your likelihood using theano you should be good to go.

def loglik_x_theano(theta, t, x):
  ...

def loglik_y_theano(theta, t, b, y):
  ...

with pm.Model() as model:
  theta_hp = pm.Bound(pm.Normal, lower=0.0, upper=3.1415)('theta_hp', 0., 1.)
  b_hp = pm.Normal(0, 10.)
  for k, (xvec, yvec, tvec) in enumerate(experimental_data):
    theta_k = pm.Bound(pm.Normal, lower=0.0, upper=3.1415)('theta_%d' % k, theta_hp, 1.)
    b_k = pm.Normal('b_%d' % k, b_hp, 1.)
    def lik_x_k(x_, t_):
      return loglik_x_theano(theta_k, t_, x_)
    def lik_y_k(y_, t_, b_):
      return loglik_y_theano(theta_k, t_, b_, y_)
    x_k = pm.DensityDist('x_%d' % k, loglik_x_theano, observed={'t_': tvec, 'x_': xvec})
    y_k = pm.DensityDist('y_%d' % k, loglik_y_theano, observed={'t_': tvec, 'b_': b, 'y_': yvec}

If your likelihood is truly black-box (you can call a non-Theano function; but cannot write it as a sequence of Theano operations) see Using pm.DensityDist and customized likelihood with a black-box function

2 Likes