# A hierarchical Bayesian model in pymc3

Suppose we have the following model:

$$X$$ unobserved

$$Y$$ such that $$Y|X \sim \mathcal{N}(X,\sigma^2)$$, observed

$$Z$$ such that $$Z|X \sim \mathcal{B}(1,X)$$, observed

and suppose, given observed data $$(y_i,z_i)_{i}$$, we want to estimate $$\sigma^2$$

To do this in a Bayesian fashion, we can put some priors on $$X$$ and $$\sigma^2$$ and run a MCMC. Consider the following pymc3 implementation (uniform prior on $$X$$ and half normal on $$\sigma^2$$)

with pm.Model() as my_model:

x = pm.Uniform('x', lower=0, upper=1)
my_sd = pm.HalfNormal('my_sd', sd=1)
y = pm.Normal('y', mu=x, sd=my_sd,observed=observed_y)
z = pm.Bernoulli('z',p=x,observed=observed_z)

with my_model:
step = pm.Metropolis()
trace = pm.sample(100000, step=step)


If I remove the line z = pm.Bernoulli('z',p=x,observed=observed_z) I get the same results. Is this expected? One could expect that providing more data that are dependent on $$X$$ would make the difference ...

You can see that \begin{align} p(\sigma \vert Y) &\propto p(\sigma)p(Y\vert\sigma) %\\&= %p(\sigma) \int_{0}^1\frac{1}{\sigma}\phi\left(\frac{x - y}{\sigma}\right)dx\\ %&= p(\sigma)\left[\Phi\left(\frac{y}{\sigma}\right) - \Phi\left(\frac{y-1}{\sigma}\right)\right], \end{align}
\begin{align*} p(\sigma \vert Y, Z) &\propto p(\sigma)p(Y, Z\vert \sigma) \\ &= p(\sigma)p(Y\vert\sigma)p(Z\vert\sigma). \end{align*}
But, the distribution of $$Z$$ does not depend on $$\sigma$$, so \begin{align*} p(\sigma \vert Y, Z) &\propto p(\sigma\vert Y), \end{align*} Therefore, $$Z$$ carries no additional information about $$\sigma.$$