# Samples from Marginal Posterior Distribution in Pymc3

Let us consider the following Hierarchical Bayesian model:

• $$mu \sim\ Beta(1, 1)$$
• $$k \sim\ Exponential(1)$$
• $$a = k*mu$$
• $$b = (1-mu) * k$$
• $$theta \sim\ Beta(a, b)$$
• $$y \sim\ Bern(theta)$$

The above example is a simplification of the example in figure 2.19 in the book Bayesian Analysis with Python by Osvaldo Martin.

I use the pymc3 to construct the model as follows:

import numpy as np
import pymc3 as pm
import arviz as az

data = np.random.binomial(1, 0.3, 10000)

with pm.Model() as model:
mu = pm.Beta('mu', 1., 1.)
k = pm.Exponential('k', 1)
a = pm.Deterministic('a', mu*k)
b = pm.Deterministic('b', (1.0-mu)*k)
theta = pm.Beta('theta', alpha=a, beta=b)
y = pm.Bernoulli('y', p=theta, observed=data)
trace = pm.sample(1000)


I have the following two questions:

• The $$trace$$ variable contains samples from the posterior distribution for the three random variables $$mu$$, $$k$$ and $$theta$$. Does the samples correspond to each of the marginal posterior distribution i.e. $$p(mu|data)$$, $$p(k|data)$$, $$p(theta|data)$$ or they correspond to the joint posterior i.e. $$p(theta,mu,k|data)$$ ?
• If the samples correspond to the joint posterior, how can I obtain samples from each of the marginal posterior distribution i.e. $$p(mu|data)$$, $$p(k|data)$$, $$p(theta|data)$$?

It contains samples from the joint posterior $$p(\mu, k, \theta| data)$$.
To obtain samples from the marginals, you can use those samples from the joint distribution: E.g., to obtain a sample from $$p(\mu|data)$$ you simply chose one of the samples from the joint distribution $$p(\mu, k, \theta| data)$$ and then you just pick from that sample the $$\mu$$-value (and throw the other two values away).