# Hierarchical Modelling

**URL:** <https://openturns.discourse.group/t/hierarchical-modelling/374>\
**Category:** Methodology\
**Tags:** calibration\
**Created:** [July 17, 2024, 11:30am UTC](https://openturns.discourse.group/t/hierarchical-modelling/374 "2024-07-17T11:30:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![guro86](https://avatars.discourse-cdn.com/v4/letter/g/47e85d/32.png) [@guro86](https://openturns.discourse.group/u/guro86)\
**Post date:** [July 17, 2024, 11:30am UTC](https://openturns.discourse.group/t/hierarchical-modelling/374/1 "2024-07-17T11:30:14Z")

</div>

Hello all,

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:

 ![Screenshot from 2024-07-17 13-23-29](https://global.discourse-cdn.com/free1/uploads/openturns/original/1X/d5b927c7158b263e0de15c120964085d3d9ad56b.png)

And the calibration is performed and published here:  
[https://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-530505](https://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-530505)

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.

Thanks,

---

<div class="post-metadata">

**Author:** ![josephmure](https://yyz2.discourse-cdn.com/free1/user_avatar/openturns.discourse.group/josephmure/32/114_2.png) [@josephmure](https://openturns.discourse.group/u/josephmure)\
**Post date:** [August 1, 2024, 6:46am UTC](https://openturns.discourse.group/t/hierarchical-modelling/374/2 "2024-08-01T06:46:12Z")

</div>

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.

```python
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:

```python
log_pdf_10_2 = ot.Function(LogPDF(10, 2))

```