Implementing Gibbs sampler in R from posterior distribution

Question

I am referencing a follow-up idea from something I posted earlier (Zero-inflated Poisson and Gibbs sampling, proofs and sampling).

I want to implement the Gibbs sampler, by generating a large (dependent) sample from the posterior distribution and use that to construct 95% Bayesian confidence intervals for $p$ and $\lambda$ using the data I generated in the first question on this page (Zero-inflated Poisson and Gibbs sampling, proofs and sampling).

Basically, I want to know how to do this in R, so that I can play around with different values of $a$ and $b$.

Answer 1

Given the three full conditionals \begin{align*} \lambda|p,\mathbf{r},\mathbf{x}&\sim Gamma\left(a+ \sum_{i}x_i, b+ \sum_{i}r_i\right)\\ p|\lambda,\mathbf{r},\mathbf{x}&\sim Beta\left(1+ \sum_{i}r_i, n+1 - \sum_{i}r_i\right)\\ r_i|\lambda,p,\mathbf{x}&\sim Bernoulli\left(\frac{pe^{- \lambda}}{pe^{- \lambda}+(1-p)\mathbb{I}{\{x_i=0}\}}\right) \end{align*} a basic Gibbs sampler would go round and round through the simulation of those three entities:

a=b=2
T=10^4
lamb=pe=rep(.5,T)
for (t in 2:T){
  r=(x==0)*(runif(n)<1/(1+(1-pe[t-1])/(pe[t-1]*exp(-lamb[t-1]))))+(x>0)
  lamb[t]=rgamma(1,a+sum(x),b+sum(r))
  pe[t]=rbeta(1,1+sum(r),n-sum(r)+1)}

leading to the following outcome (based on data simulated with $p=0.3$ and $\lambda=2$: