Question 17. Hierarchical Poisson model
Suppose J factories manufacture components. Factory j produces n_j observations of the number of defective components per batch.
Assume
Y_{ij}\mid\theta_j\stackrel{\mathrm{ind}}{\sim}\operatorname{Poisson}(\theta_j), \qquad i=1,\ldots,n_j,\quad j=1,\ldots,J.
Here, \theta_j>0 represents the expected number of defects per batch in factory j.
Rather than modelling the factory rates independently, we assume
\theta_j\mid\lambda\overset{\text{iid}}{\sim}\operatorname{Gamma}(a,\lambda),
where a>0 is known and \lambda>0 is an unknown common hyperparameter.
Finally, assign the hyperprior
\lambda\sim\operatorname{Gamma}(c,d),
where c,d>0 are known.
All Gamma distributions use the shape–rate parameterisation:
f(x\mid a,b)=\frac{b^{a}}{\Gamma(a)}x^{a-1}e^{-bx}, \qquad x>0.
(a) Write down the joint density p({\boldsymbol{Y}},{\boldsymbol{\theta}},\lambda) up to proportionality in the unknown parameters.
(b) Derive the full conditional distribution p(\theta_j\mid\lambda,{\boldsymbol{Y}}) for each factory j. Identify its distribution and parameters.
(c) Derive the full conditional distribution p(\lambda\mid{\boldsymbol{\theta}},{\boldsymbol{Y}}). Identify its distribution and parameters.
(d) Using your results from (b) and (c), describe a Gibbs sampling algorithm targeting p({\boldsymbol{\theta}},\lambda\mid{\boldsymbol{Y}}).
(e) Explain how the hierarchical model allows observations from one factory to influence inference about another factory’s defect rate.
Definitions. Y_{ij} is the count from batch i in factory j, \theta_j is that factory’s mean count, and \lambda is a shared unknown gamma rate. The known shape is a, while c,d define the shape–rate gamma hyperprior on \lambda. Write s_j=\sum_iY_{ij} and \bar Y_j=s_j/n_j. The conditional hierarchy is \lambda\longrightarrow{\boldsymbol{\theta}}\longrightarrow{\boldsymbol{Y}}.
(a) Let s_j=\sum_{i=1}^{n_j}Y_{ij}. Since {\boldsymbol{Y}}\perp\lambda\mid{\boldsymbol{\theta}},
p({\boldsymbol{Y}},{\boldsymbol{\theta}},\lambda) = p({\boldsymbol{Y}}\mid{\boldsymbol{\theta}})\,p({\boldsymbol{\theta}}\mid\lambda)\,p(\lambda).
For \theta_j>0 and \lambda>0,
\begin{aligned} p({\boldsymbol{Y}}\mid{\boldsymbol{\theta}}) &= \Biggl(\prod_{j=1}^{J}\prod_{i=1}^{n_j}\frac{1}{Y_{ij}!}\Biggr) \prod_{j=1}^{J}e^{-n_j\theta_j}\theta_j^{s_j},\\ p({\boldsymbol{\theta}}\mid\lambda) &= \Biggl(\frac{\lambda^{a}}{\Gamma(a)}\Biggr)^{J} \prod_{j=1}^{J}\theta_j^{a-1}e^{-\lambda\theta_j},\\ p(\lambda) &= \frac{d^{c}}{\Gamma(c)}\lambda^{c-1}e^{-d\lambda}. \end{aligned}
The factors 1/Y_{ij}!, \Gamma(a)^{-J} and d^{c}/\Gamma(c) do not depend on ({\boldsymbol{\theta}},\lambda). Dropping them leaves
\boxed{ p({\boldsymbol{Y}},{\boldsymbol{\theta}},\lambda) \propto \lambda^{aJ+c-1}e^{-d\lambda} \prod_{j=1}^{J} \theta_j^{s_j+a-1}e^{-(n_j+\lambda)\theta_j}. }
(b) In the kernel from (a), the factors containing \theta_j are
\theta_j^{s_j+a-1}e^{-(n_j+\lambda)\theta_j} = \theta_j^{(s_j+a)-1}e^{-(n_j+\lambda)\theta_j}, \qquad\theta_j>0.
No other rate appears. Here s_j+a>0 because a>0 and s_j\ge 0, and n_j+\lambda>0 because \lambda>0. This is a shape–rate gamma kernel, so
\boxed{ \theta_j\mid\lambda,{\boldsymbol{Y}} \sim \operatorname{Gamma}\!\left(a+\sum_{i=1}^{n_j}Y_{ij},\ n_j+\lambda\right). }
(c) In the kernel from (a), the factors containing \lambda are
\lambda^{aJ+c-1} \exp\Bigl(-\lambda\sum_{j=1}^{J}\theta_j\Bigr)e^{-d\lambda} = \lambda^{(aJ+c)-1} \exp\Bigl(-\Bigl(d+\sum_{j=1}^{J}\theta_j\Bigr)\lambda\Bigr),
for \lambda>0. Given {\boldsymbol{\theta}}, the observations do not appear. The shape aJ+c>0, and the rate d+\sum_j\theta_j>0 because d>0. Thus
\boxed{ \lambda\mid{\boldsymbol{\theta}},{\boldsymbol{Y}} \sim \operatorname{Gamma}\!\left(c+aJ,\ d+\sum_{j=1}^{J}\theta_j\right). }
(d) By (b) and (c), alternate a block of factory-rate draws with a shared-rate draw. Recall s_j=\sum_iY_{ij}.
All gamma distributions in the algorithm use rates. The J factory draws are conditionally independent given the previous \lambda, so their order does not matter. The final draw uses all newly updated rates. These exact full-conditional updates preserve p({\boldsymbol{\theta}},\lambda\mid{\boldsymbol{Y}}); assess convergence and mixing before using the retained output.
(e) By (b), \theta_j depends on the other factories only through \lambda. Counts from factory k update \theta_k, and by (c) every rate enters the conditional distribution of the shared \lambda. That \lambda is the second parameter of \theta_j. This dependence through the common hyperparameter is borrowing strength: observations from one factory inform the defect rate at another.
At a fixed \lambda, the conditional posterior mean makes the shrinkage explicit:
\mathbb{E}(\theta_j\mid\lambda,{\boldsymbol{Y}}) =\frac{s_j+a}{n_j+\lambda} =\frac{n_j}{n_j+\lambda}\bar Y_j +\frac{\lambda}{n_j+\lambda}\frac a\lambda.
This averages factory j’s sample mean and the common prior mean a/\lambda. For a given \lambda, smaller factories receive more pooling. The marginal posterior mean averages this expression over p(\lambda\mid{\boldsymbol{Y}}), which is informed by all factories. Fixing \lambda externally removes this route for other factories’ data to influence \theta_j. Integrating over a shared uncertain \lambda also generally makes the factory rates dependent in the joint posterior.