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?