I have a hierarchical probabilistic model that I calibrated successfully using my own implementation of a mh-within-Gibbs sampler. The model looks like this:
I now want to implement the same application in openTurns since it includes tuning proposal covariances, and the implementation will be more general. However, I struggle to understand how one could implement a hierarchical model. I.e., where we have mean values and standard deviation and then another layer of distributions that depend on those parameters, those parameters, in turn, get used into a surrogate that finally is conditioned on observations through a normal likelihood.
I would be grateful if someone could give me an example of a hierarchical implementation or other ideas, tips, or tricks.
Hi @guro86 we talked about your use case yesterday but I would like to add one thing: if you really want to always use the smallest possible implementation of the conditional log-pdf for each Metropolis-Hastings sampler, you can define a kind of “template” function with the class OpenTURNSPythonFunction.
The idea is to implement a child class which will then be wrapped in a regular Function object.
Let us assume you have implemented a method x_index such that x_index(i,j)gets you the index of x_{i,j} in the state (i.e. the list with length 99 = 3 \times 2 + 3\times31). Let us also assume you have implemented a similar mu_index(j) to get the index of \mu_j and a similar sigma_index(j) method to get the index of \sigma_j.
import openturns as ot
class LogPDF(ot.OpenTURNSPythonFunction):
def __init__(self, arg_number, obs_number):
"""Setup the input and output dimensions, and also internal variables."""
# arg_number is j = 1,2 or 3
# obs_number is i = 1..31
super().__init__(99, 1)
self._arg_number = arg_number
self._obs_number = obs_number
def _exec(self, state):
"""_exec is what will actually run when the function is called"""
muj = state[mu_index(self._arg_number)]
sigmaj = state[sigma_index(self._arg_number)]
xij = state[x_index(self._obs_number, self._arg_number)]
xi = [state[x_index(self._obs_number, j) for j in range(1, 4)]
logpdf = # sum of the logpdf of xij conditioned on muj and sigmaj and of the log-likelihood of xi=(xi1, xi2, xi3)
return [logpdf]
Then, to build for example the conditional log-pdf of, e.g., the x_{10,2} sampler, you can simply write:
