while Liu, Jun S. (2008) introduces two types of Gibbs sampling strategy.
Bivariate Gibbs sampler
It is easy to implement this sampler:
## Julia program for Bivariate Gibbs sampler## author: weiya <szcfweiya@gmail.com>## date: 2018-08-22functionbigibbs(T, rho) x =ones(T+1) y =ones(T+1)for t =1:T x[t+1] =randn() *sqrt(1-rho^2) + rho*y[t] y[t+1] =randn() *sqrt(1-rho^2) + rho*x[t+1]endreturn x, yend## examplebigibbs(100,0.5)
Completion Gibbs Sampler
Example:
We can use the following Julia program to implement this algorithm.
## Julia program for Truncated normal distribution## author: weiya <szcfweiya@gmail.com>## date: 2018-08-22# Truncated normal distributionfunctionrtrunormal(T, mu, sigma, mu_down) x =ones(T) z =ones(T+1)# set initial value of z z[1] =rand()if mu < mu_down z[1] = z[1] *exp(-0.5* (mu - mu_down)^2/ sigma^2)endfor t =1:T x[t] =rand() * (mu - mu_down +sqrt(-2*sigma^2*log(z[t]))) + mu_down z[t+1] =rand() *exp(-(x[t] - mu)^2/(2*sigma^2))endreturn(x)end## examplertrunormal(1000,1.0,1.0,1.2)
and Liu, Jun S. (2008) also presents the slice sampler with slightly different expression:
In my opinion, we can illustrate this algorithm with one dimensioanl case. Suppose we want to sample from normal distribution (or uniform distribution), we can sample uniformly from the region encolsed by the coordinate axis and the density function, that is a bell shape (or a square).
Consider the normal distribution as an instance.
It is also easy to write the following Julia program.
## Julia program for Slice sampler## author: weiya <szcfweiya@gmail.com>## date: 2018-08-22functionrnorm_slice(T) x =ones(T+1) w =ones(T+1)for t =1:T w[t+1] =rand() *exp(-1.0* x[t]^2/2) x[t+1] =rand() *2*sqrt(-2*log(w[t+1])) -sqrt(-2*log(w[t+1]))endreturn x[2:end]end## examplernorm_slice(100)
Data Augmentation
A special case of Completion Gibbs Sampler.
Let's illustrate the scheme with grouped counting data.
And we can obtain the following algorithm,
But it seems to be not obvious to derive the above algorithm, so I wrote some more details
Liu, Jun S. (2008) also presents the DA algorithm which based on Bayesian missing data problem.
Then he argues that m copies of ymis in each iteration is not really necessary. And briefly summary the DA algorithm:
It seems that we do not need to derive the explicit form of g(x,z), if we can directly obtain the conditional distribution. We can use the following Julia program to sample.