Beta-Binomial Gibbs Sampler

Question

I am self-studying Bayesian statistics from the book Computational Bayesian Statistics by Turkman et al, but I am stuck on Problem 6.3 from the book:

Suppose we want to consider a Binomial (unknown $\theta, n) \land$ Beta model, in particular, a binomial sampling model $x \mid \theta, n \sim \text{Bi}(n, \theta)$. Assume that $\theta \sim \text{Be}(a_0,b_0)$ and $h(n) \propto \frac{1}{n^2}$. Find the conditional posterior distributions $h(\theta \mid n, x)$ and $h(n \mid \theta, x)$.

Next, describe and implement a Gibbs sampling algorithm to generate $(n_m, \theta_m) \sim h(n, \theta \mid x)$. Plot the joint posterior $h(n, \theta \mid x)$, and plot on top of the same figure the simulated posterior draws $(n_m , \theta_m), m = 1, ..., 50$ (connected by line segments showing the moves). Use $x = 50$, and $(a_0,b_0) = (1,4)$, a grid on $0.01 \leq \theta \leq 0.99$ and $x \leq n \leq 500$. Use the R function lgamma(n+1) to evaluate $\ln{(n!)}$.

Finally, implement Metropolis-Hastings posterior simulation. Add the simulated posterior draws on top of the plot from above.

I was able to find the conditional posterior distributions $h(\theta \mid n, x)$ and $h(n \mid \theta, x)$ to be $$h(\theta \mid n, x) = \text{Be}(a_0 + x, b_0 + n - x)$$ and $$h(n, \theta \mid x) \propto \frac{1}{n^2} \theta^{a_0 + x - 1}(1 - \theta)^{n - x + b_0 - 1} \frac{n!}{(n - x)!},$$ but I'm not sure how to implement the Gibbs sampler in this case. I'm new to R, so I'd appreciate any help in this manner! Thanks in advance.

Answer 1

You have the conditional posterior $$h(n | \theta, x) \propto n^{-2} (1-\theta)^{n} \frac{n!}{(n-x)!}.$$ To sample from $h(n, \theta|x)$ using Gibbs sampling, you take turns sampling from $h(\theta|x,n)$ and $h(n |\theta, x)$. So your problem boils down to sampling from $h(n|\theta,x)$.

A first approach is an inexact method, but offers a taste of the essence of solutions. The approach is as follows. Compute explicitly $h(n | \theta, x$) for $n=1,2,\dots,M$ where $M$ is some very large integer like $100000$, then you can sample using R's sample() function. This will introduce a truncation bias, as you're sampling from a distribution which is unsupported for values larger than $M$.

The better approach is to use inverse transform sampling$^\dagger$. The idea is to sample a value $u$ uniformly from $[0,1]$, then find the inverse CDF to find a value $n$ such that $\text{CDF}(n | \theta, x) = u$, and take $n$ to be your realization from $h(n | \theta, x)$. $$\text{CDF}(n|\theta,x) = \sum_{i=1}^n h(i|\theta,x).$$

$^\dagger$ you may find the discussions in this question on inverse sampling helpful.