I am self-learning Bayesian statistics with the book Computational Bayesian Statistics by Turkman et al. and I am currently stuck on Problem 6.5:
We consider a hierarchical event rate model with a Poisson sampling model $y_i \sim \text{Poi}(\lambda_i t_i)$, and a prior model $\lambda_i \sim \text{Ga}(\alpha, \beta)$ and $\beta \sim \text{Ga}(c,d)$, where $(\alpha, c, d)$ are fixed parameters.
- Find the conditional posterior distributions: $h(\lambda_i \mid \beta, \lambda_j, j \neq i,y)$ and $h(\beta \mid \lambda_1, …, \lambda_n,y)$.
- Using the data below, implement a Gibbs sampler to simulate from the posterior in the above model.
Table showing pump failures ($y_i$) over given exposure times ($t_i$), $i = 1,…,n = 10$:
| $y_i$ | $5$ | $1$ | $5$ | $14$ | $5$ | $19$ | $1$ | $1$ | $4$ | $22$ |
|---|---|---|---|---|---|---|---|---|---|---|
| $t_i$ | $94.32$ | $15.72$ | $62.88$ | $125.76$ | $5.24$ | $31.44$ | $1.048$ | $1.048$ | $2.096$ | $10.48$ |
- Alternatively, marginalize with respect to $\lambda_i$ to find the marginal likelihood $f(y_i \mid \beta)$ and $h(\beta \mid y)$.
I'm not sure where to start for 1 and 2, and I don't have much experience using R and implementing a Gibbs sampler, so I'd appreciate any hints/help! Thanks in advance.
Best Answer
The conditional denisty $p(\lambda_i|\beta,\lambda_j,j\neq i,y)$ is given by the following formula: $$p(\lambda_i|\beta,\lambda_j,j\neq i,y)=\frac{p(\lambda_1,\dots,\lambda_n,y,\beta)}{\int_0^{\infty}p(\lambda_1,\dots,\lambda_n,y,\beta)\mathrm{d}\lambda_i}$$ The joint density $p(\lambda_1,\dots,\lambda_n,y,\beta)$ apparent in the previous expression factors accordingly: $$\begin{eqnarray*}p(\lambda_1,\dots,\lambda_n,y,\beta)&=&p(y|\lambda_1,\dots,\lambda_n,\beta)p(\lambda_1,\dots,\lambda_n|\beta)p(\beta) \\ &=&e^{-(\lambda_1 t_1 + \dots + \lambda_n t_n)}\cdot \frac{(\lambda_1t_1)^{y_1}\times\dots\times(\lambda_n t_n)^{y_n}}{y_1!\times \dots \times y_n!}\cdot \left(\frac{\beta^{\alpha}}{\Gamma(\alpha)}\right)^n(\lambda_1 \times \dots \times \lambda_n)^{\alpha-1}e^{-\beta(\lambda_1 + \dots + \lambda_n)}\cdot \frac{d^c}{\Gamma(c)}\beta^{c-1}e^{-d\beta}\end{eqnarray*}$$ Plug this expression into our formula for $p(\lambda_i|\beta,\lambda_j,j\neq i,y)$ and simplify to get $$p(\lambda_i|\beta,\lambda_j,j\neq i,y)=\frac{\lambda_j^{\alpha+y_j-1}e^{-(\beta+t_j)\lambda_j}}{\int_0^{\infty}\lambda_j^{\alpha+y_j-1}e^{-(\beta+t_j)\lambda_j}\mathrm{d}\lambda_j}$$ You should now recognize that $$\lambda_i|\beta,\lambda_j,j\neq i,y\sim \text{Gamma}\left(\alpha+y_i,\beta+t_i\right)$$ Meanwhile, $$p(\beta|\lambda_1,\dots,\lambda_n,y)=\frac{p(\lambda_1,\dots,\lambda_n,y,\beta)}{\int_0^{\infty}p(\lambda_1,\dots,\lambda_n,y,\beta)\mathrm{d}\beta}$$ Plugging in our joint density and simplifying yields a simpler expression for this conditional distribution, namely $$p(\beta|\lambda_1,\dots,\lambda_n,y)=\frac{\beta^{n\alpha+c-1}e^{-(\lambda_1 + \dots + \lambda_n +d)\beta}}{\int_0^{\infty}\beta^{n\alpha+c-1}e^{-(\lambda_1 + \dots + \lambda_n +d)\beta}\mathrm{d}\beta}$$ It's now clear that $$\beta|\lambda_1,\dots,\lambda_n,y \sim \text{Gamma}\left(n\alpha +c, \lambda_1 + \dots + \lambda _n +d\right)$$ This is your answer for (1) so now let's take a look at (3). From the total law we gather that $$\begin{eqnarray*}p(\lambda_1,\dots,\lambda_n)&=& \int_0^{\infty}p(\lambda_1,\dots,\lambda_n|\beta)p(\beta)\mathrm{d}\beta \\ &=& \frac{(\lambda_1 \times \dots \times \lambda_n)^{\alpha-1}d^c}{\left(\Gamma(\alpha)\right)^n\Gamma(c)} \int_0^{\infty}\beta^{n\alpha +c -1}e^{-(\lambda_1 + \dots + \lambda_n +d)\beta}\mathrm{d}\beta \\&=& \frac{\Gamma\left(n\alpha+c\right)(\lambda_1 \times \dots \times \lambda_n)^{\alpha-1}d^c}{\left(\Gamma(\alpha)\right)^n\Gamma(c)(\lambda_1 + \dots + \lambda_n +d)^{n\alpha+c}}\end{eqnarray*}$$Notice from Baye's we have that $$\begin{eqnarray*}p(\lambda_1,\dots,\lambda_n|y)&=&\frac{p(y|\lambda_1,\dots,\lambda_n)p(\lambda_1,\dots,\lambda_n)}{\int_{(0,\infty)^n}p(y|\lambda_1,\dots,\lambda_n)p(\lambda_1,\dots,\lambda_n)\mathrm{d}\lambda_1\dots \mathrm{d}\lambda_n} \\ &=& C(y_1,\ldots,y_n)\cdot \frac{e^{-(\lambda_1t_1 + \dots + \lambda_n t_n)}\cdot (\lambda_1^{\alpha+y_1-1}\times \dots \times \lambda_n^{\alpha+y_n-1})}{(\lambda_1 + \dots + \lambda_n +d)^{n\alpha+c}}\end{eqnarray*}$$ where $$C(y_1,...,y_n)=\left(\int_{(0,\infty)^n}\frac{e^{-(\lambda_1t_1 + \dots + \lambda_n t_n)}\cdot (\lambda_1^{\alpha+y_1-1}\times \dots \times \lambda_n^{\alpha+y_n-1})}{(\lambda_1 + \dots + \lambda_n +d)^{n\alpha+c}}\mathrm{d}\lambda_1 \dots \mathrm{d}\lambda_n\right)^{-1}$$You can use total law to find $p(\beta|y)$: $$\begin{eqnarray*}p(\beta|y)&=& \int_{(0,\infty)^n}p(\beta|y,\lambda_1,\dots,\lambda_n)p(\lambda_1,\dots,\lambda_n|y)\mathrm{d}\lambda_1 \dots \mathrm{d}\lambda_n \\ &=& \frac{C(y_1,...,y_n)\Gamma(\alpha+y_1)\times \dots \times \Gamma(\alpha+y_n)}{\Gamma(n\alpha+c)}\cdot \frac{\beta^{n\alpha +c-1}e^{-d\beta}}{(\beta+t_1)^{\alpha+y_1}\times \dots \times (\beta+t_n)^{\alpha+y_n}} :\beta>0\end{eqnarray*}$$ Using the fact that $$\begin{eqnarray*}p(y_i|\beta) &=& \int_{(0,\infty)^n}p(y_i|\beta,\lambda_1,\dots,\lambda_n)p(\lambda_1,\dots,\lambda_n|\beta)\mathrm{d}\lambda_1 \dots \mathrm{d}\lambda_n \end{eqnarray*}$$ you will see $$p(y_i|\beta)=\frac{\Gamma(\alpha+y_i)}{y_{i}!\Gamma(\alpha)}\left(\frac{t_{i}}{t_{i}+\beta}\right)^{y_i}\left(\frac{\beta}{t_{i}+\beta}\right)^{\alpha}:y_i=0,1,2,\dots$$