# Beta-Binomial Gibbs Sampler

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.

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.

• Thanks for the help! I understand how to sample $h(\theta \mid x,n)$ using the rbeta() function, but I am still unsure of $h(n \mid \theta, x)$ - why would I need to use the sample() function? Won't rbinom() with size $n$ and probability $\theta$ be enough?
– user310180
Feb 14, 2022 at 5:52
• No, rbinom won't do because $h(n| \theta, x)$ isn't binomial (look at its pmf).
– fool
Feb 14, 2022 at 15:18