3
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ If you only have one parameter Gibbs sampling doesn't make sense, in this case you'd just be sampling directly from the posterior. Try looking at a problem with multiple parameters. $\endgroup$ – aleshing Mar 2 '18 at 1:23
  • 2
    $\begingroup$ You mention MacKay, are you referencing a paper or book? Could you add this to the question? $\endgroup$ – Jon Mar 2 '18 at 6:13
  • $\begingroup$ information theory, inference and learning algorithms $\endgroup$ – wirrbel Mar 2 '18 at 7:14
4
$\begingroup$

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

enter image description here

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I was hoping to solve this without the analytical posterior distribution. Basically performing sampling to obtain $\Pr(k\mid\theta)$. $\endgroup$ – wirrbel Mar 2 '18 at 7:12
  • $\begingroup$ @Xi'an thank you, I had misunderstood it's purpose then. $\endgroup$ – wirrbel Mar 2 '18 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.