sampling behind bayesian hierarchical models

Question

I'm unsure how sampling is done in Bayesian Hierarchical modelling, i'm reading a book on how to use it in PyMC3 but it doesn't explain the math and i'd like to understand it.

Suppose i want to estimate the parameters $\theta_{1},\theta_{2}$ and $\theta_{3}$ and i have samples $y_{1},y_{2},y_{3}$ corresponding to the model involving parameter $\theta_{1}$, the samples $x_{1},x_{2},x_{3},x_{4},x_{5}$ corresponding to the model involving parameter $\theta_{2}$ and samples $z_{1},z_{2}$ corresponding to the model involving parameter $\theta_{3}$. In addition suppose $\theta_{i} \sim Beta(\alpha,\beta)$ with $\alpha \sim HalfCauchy(\beta_{\alpha})$ and $\beta \sim HalfCauchy(\beta_{\beta})$. The samples for each $\theta_{i}$ are obtained bernoulli i.e $bern(\theta_{i})$.

How do we obtain the posteriors for each $\theta_{i}$? Do we feed the data $y_{1},\ldots,y_{3}$, resp($x_{1},\ldots,x_{5}$ and $z_{1},z_{2}$) in separately or together as one? I would have guessed i would feed the corresponding data seperately to obtain each estimate $\theta_{i}$ but it appears in the code it is done all at once. Can someone explain how the math works?

Here is the scenario that was detailed in the book along with the code, not i've changed the samples amounts in each of the cases to not be even as they are in the example.

"To illustrate the main concepts of hierarchical models, we are going to use a toy model of the water quality example we discussed at the beginning of this section, and we are going to use synthetic data. Imagine we have collected water samples from three different regions of the same city and we have measured the lead content of water; samples with lead concentration above recommendations from the World Health Organization (WHO) are marked with zero and samples with the following values are marked with one. This is just a pedagogic example; in a more realistic example, we would have a continuous measurement of lead concentration and probably many more groups. Nevertheless, for our current purposes, this example is good enough to uncover the details of hierarchical models."

Data Generation

N_samples =  [30, 30, 30]
G_samples =  [18, 18, 18]
group_idx = np.repeat(np.arange(len(N_samples)), N_samples)
data = []  
for i in range(0, len(N_samples)):    
    data.extend(np.repeat([1, 0], [G_samples[i], N_samples[i]-G_samples[i]]))

Bayesian Modelling

with pm.Model() as model_h:    
    alpha = pm.HalfCauchy('alpha', beta=10)    
    beta = pm.HalfCauchy('beta', beta=10)
    theta = pm.Beta('theta', alpha, beta, shape=len(N_samples))        
    y = pm.Bernoulli('y', p=theta[group_idx], observed=data)

    trace_j = pm.sample(2000)
chain_h = trace_h[200:]
pm.traceplot(chain_h)

I'd really love to know what is going on under the hood of the model how are the posteriors calculated? How is the data fed into the model? If for example to obtain the posterior for $\theta_{1}$ i should only feed in $y_{1},y_{2},y_{3}$ why does it look like all the data is aggregated and fed into the model-model_h?

In addition another question i have is how does the PyMC3 code distinguish which portions of the data belong to which $\theta_{i}$ since we feed in a long list of data in the variable under the variable name data which is all observations for all $\theta$? The guess that i have is that it could be due to the group_idx, which gives an index for which $\theta_{i}$ the data observation belongs to. Is this correct?

Answer 1

Here is the basic structure of a hierarchical model. In order to simplify the exposition, I'm going to modify the notation a bit.

Let there be $n$ groups (or units), $Y = (Y_1, \ldots, Y_n)$, where $Y_i = (y_{i1}, \ldots, y_{iT_i})$. Let the sampling distribution be given by \begin{equation} p(Y|\theta) = \prod_{i=1}^n p(Y_i|\theta_i), \end{equation} where $\theta = (\theta_1,\ldots,\theta_n)$ \begin{equation} p(Y_i|\theta_i) = \prod_{t=1}^{T_i} p(y_{it}|\theta_i) . \end{equation} Let the prior for $\theta$ be given by \begin{equation} p(\theta|\phi) = \prod_{i=1}^n p(\theta_i|\phi) , \end{equation} and let the prior for $\phi$ be given by $p(\phi)$.

In the example, $\phi = (\alpha,\beta)$ and $p(\phi)$ involves two half-Cauchy distributions. In addition, $p(\theta_i|\phi)$ is a Beta distribution and $p(y_{it}|\theta_i)$ is a Bernoulli distribution.

The joint posterior for all the unknowns given all the data is characterized by \begin{equation}\tag{A} p(\theta,\phi|Y) \propto p(Y,\theta,\phi) = p(Y|\theta)\,p(\theta|\phi)\,p(\phi). \end{equation} The joint posterior distribution can be characterized by the following collection of full conditional distributions: \begin{align}\tag{B} p(\theta|Y,\phi) &\propto p(Y|\theta)\,p(\theta|\phi) = \prod_{i=1}^n p(Y_i|\theta_i)\,p(\theta_i|\phi) \\ \tag{C} p(\phi|Y,\theta) &\propto p(\theta|\phi)\,p(\phi) = \bigg(\prod_{i=1}^n p(\theta_i|\phi)\bigg) p(\phi) . \end{align} The right-hand sides display important simplifications. To begin with, each of (B) and (C) omits unnecessary factors: $p(\phi)$ is omitted from (B) while $p(Y|\theta)$ is omitted from (C). The right-hand side of (B) shows that each $\theta_i$ occurs in two factors by itself and no other. Consequently the posteriors for the $\theta$ parameters (conditional on $\phi$) are independent of each other. (This follows from the "product rule.") The right-hand side of (C) shows that the conditional posterior for the hyperparameter does not involve the data $Y$. Instead, the parameters $\theta$ play the role of the "observations."