# Implementing Gibbs sampler in R from posterior distribution

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$.

• You have the three full conditionals, where is the difficulty for you? (Please add self-study as a tag.) Mar 26 '15 at 6:30
• I guess the difficulty for me is knowing how to generate this in R. I am new to programming in R. I wanted to work on learning R and simulating this data to better understand it at the same time. Yes, I hope to supplement this all with self-study, but learning a new programming language (at least for me) is a slow process and that isn't always linear (I use multiple sources, including reading code written by others, and then dissecting that code to learn what functions they used etc.) Mar 26 '15 at 9:36
• You should first learn R then, since this question has to do with R and not with Gibbs sampling. Mar 26 '15 at 10:01

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$:

• Thanks, @Xi'an. The only things I am not sure about after looking at the code is 1) Is it using Bayesian confidence intervals of 95%, 2) Is it using the data generated in the first part of this set? (stats.stackexchange.com/questions/143468/…). In any case, thanks again. Mar 26 '15 at 10:51
• did you run the R code yourself? then you should find an answer to your questions... Mar 26 '15 at 10:56