Gibbs sampling and Bayesian inference

Question

I wanted to get a more in-depth understanding of sampling algorithms, and so I thought I could start with the very simple example of a binomial or bernoulli likelihood with a beta prior (since it is solvable in closed form, I thought it might be suitable).

Assume we have trials where have draw $N$ times, and count $k$ successes. $N$ is a deterministic parameter predefined in the experimental setting. We model $k \mid \theta, N \sim \text{Binomial}(\theta, N)$ and model $\theta \sim \text{Beta}(\alpha_0, \beta_0)$. We would like to obtain $\Pr(\theta \mid k, N)$ by sampling, given several observations of $k$, without relying on the closed-form solution of a Beta-posterior.

Looking at MacKay's Information theory, inference and learning algorithms, we would start with an initial $x=(k_0, \theta_0)$, and use the $\text{Binomial}(k\mid \theta_0, N)$ to draw a new $k_1$. This part I can understand, because it is drawing a random sample from a binomial distribution. Then we would use $\Pr(\theta\mid k_1, N)$ to draw $\theta_1$. This part I do not understand anymore, because at this point, we don't have a posterior yet, don't we? So without a posterior distribution $\Pr(\theta\mid k, N)$, how can we sample a new $\theta$?. And furthermore, **shouldn't the algorithm also sample from the observed $k$ values? Or am I missing a step?

Answer 1

If the distribution to simulate is the posterior $$\pi(\theta|k_0) \propto {N \choose k_0} \theta^{k_0+\alpha_0 - 1}(1-\theta)^{N-k_0+\beta_0 - 1}\propto \theta^{k_0+\alpha_0 - 1}(1-\theta)^{N-k_0+\beta_0 - 1},$$a Gibbs sampler can be constructed by introducing auxiliary variables and making the above density function an integral over these auxiliary variables: $$\pi(\theta|k_0) = \int \pi(\theta,u|k_0)\text{d}u$$ A standard completion algorithm is the slice sampler, where terms in the density are replaced by integrals of indicator functions: \begin{align*}\theta^{k_0+\alpha_0-1}&=\int_0^\infty \mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\text{d}u_1\\ (1-\theta)^{N-k_0+\beta_0 - 1}&=\int_0^\infty \mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\text{d}u_2\\\end{align*} which makes the joint distribution equal to $$\pi(\theta,u|k_0)\propto\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\qquad u=(u_1,u_2)$$ The associated Gibbs sampler is then \begin{align*}U_1|\theta,k_0,U_2 &\sim \mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}} \\ U_2|\theta,k_0,U_2 &\sim \mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\\ \theta|k_0,U_1,U_2 &\sim\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\\\end{align*} Meaning \begin{align*}U_1|\theta,k_0,U_2 &\sim \mathcal{U}(0,\theta^{N-k_0+\beta_0 - 1}) \\ U_2|\theta,k_0,U_2 &\sim \mathcal{U}(0,(1-\theta)^{N-k_0+\beta_0 - 1}) \\ \theta|k_0,U_1,U_2 &\sim \mathcal{U}(U_1^{1/{k_0+\alpha_0-1}},1-U_2^{1/{N-k_0+\beta_0 - 1}})\\\end{align*} since $$\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}=\mathbb{I}_{u_1^{1/{k_0+\alpha_0-1}}\le\theta}\mathbb{I}_{u_2^{1/{N-k_0+\beta_0 - 1}}\le(1-\theta)}$$

An R code implementing this slice sampler is

  sliz=function(T=1e3){
   the0=k0/N #starting value
   thez=rep(0,T)
   for (t in 1:T){
     u1=runif(1)^{1/{k0+a0-1}}*the0
     u2=runif(1)^{1/{N-k0+b0-1}}*(1-the0)
     thez[t]=the0=runif(1,u1,1-u2)}
   return(thez)}

resulting in a neat fit to the intended target

Answer 2

The joint distribution $$ f(\theta, k) \propto {N \choose k} \theta^{k+\alpha_0 - 1}(1-\theta)^{N-k+\beta_0 - 1} \tag{1}. $$ That means $$ f(\theta \mid k) = \text{Beta}(k + \alpha_0, N-k+\beta_0) \tag{2}. $$ Gibbs sampling in this case would involve iterations that alternate by drawing from $f(\theta \mid k)$ and $f(k \mid \theta)$ (which is given). This would give you draws from a Markov chain with a stationary distribution equal to the joint distribution above. They would look something like: $$ (\theta_0, k_0),(\theta_0, k_1) ,(\theta_1, k_1) ,(\theta_1, k_2) ,\ldots $$

You are correct that it's strange to do Gibbs sampling in this situation because you are targeting a joint distribution, and usually in practice you target the posterior. The posterior is completely known (2), and the joint is also known ((1) is a Beta-Binomial). However, this is a common example used for for illustrative purposes.