Generating samples from Gibbs sampling

Question

I am quite new to sampling. I am doing Gibbs sampling for a Bayesian network. I am aware about the algorithm for the Gibbs sampling but there's one thing I am not able to understand.

For example let's assume you have 3 variables $x_1, x_2, x_3$ means you need to generate sample of the form $(x_1,x_2,x_3)$. And suppose that all three parameters can take 3 values 0,1 or 2. Now, let's take initial value of the sample as (2,1,1) and compute the $\mathbb P(x_1=0 | x_2=1,x_3=1)$, $\mathbb P(x_1=1 |x_2=1,x_3=1)$ and $\mathbb P(x_1=2 |x_2=1,x_3=1)$.

Now, I am selecting the largest probability among these 3 and assigning that value to $x_1$ and using it for the next sample. I want to ask if this method for generating model is correct or I am doing something wrong.

Answer 1

In your simple example, once you've computed $\mathbb P(x_1=0 | x_2=1,x_3=1)=p_{0,11}$, $\mathbb P(x_1=1 |x_2=1,x_3=1)=p_{1,11}$ and $\mathbb P(x_1=2 |x_2=1,x_3=1)=p_{2,11}$, you draw $x_1=0$ with probability $p_{0,11}$, $x_1=1$ with probability $p_{1,11}$, and $x_1=2$ with probability $p_{2,11}$. You can operationalize this as

1. Draw $U\sim U[0,1]$
2. If $U \le p_{0,11}$, assign $x_1=0$, and you are done.
3. If $p_{0,11} < U \le p_{0,11} + p_{1,11}$, assign $x_1=1$, and you are done.
4. If $ p_{0,11} + p_{1,11}<U$, assign $x_1=2$.

It is easy to see that thus generated $x_1$ has the required distribution.

If you drew the modal value, you will quickly get stuck in the local mode of the joint distribution, which is not what you want.

Answer 2

I have posted a link to a previous great answer which is more of an intuitive nature. I am going to focus on more technical details and to provide a simple example implemented in R.

Gibbs sampling is a Markov Chain Monte Carlo (MCMC) method that is used to simulate from a random vector $(X_1,...,X_n)$. Suppose that you want to obtain a simulated sample of size $N$ from this vector and denote the $j$th simulation as $(x_1^j,...,x_n^j)$. Assume also that you know the conditional distributions $X_j\vert X_1,{j-1},X_{j+1},...,X_n$ for $j=1,...,n$ and that you can simulate from them.

The Gibbs sampler is defined as the following iterative algorithm

1. Start at the initial value $(x_1^0,...,x_n^0)$.
2. For $k=1,...,N$ obtain a simulation of each conditional $X_j^{k}\vert x_1^{k},...,x_{j-1}^k,x_{j+1}^{k-1},...,x_{n}^{k-1}$. For this step it is essential to be able to simulate from each univariate conditional for the model parameters of interest.

With this algorithm you will obtain a sample $\{(x_1^j,...,x_n^j)\}$ of size $N$ from the joint distribution of $(X_1,...,X_n)$. Note that there is no need to evaluate any probability.

In order to ensure that the sample is approximately independent and from the target distribution you have to ignore the first $M$ samples, for a large $M$, this is called burn-in, and also to subsample every $m$ iterations, this is called thinning.

In your particular case, you have to be able to simulate from the conditional distributions of $X_1\vert X_2,X_3$; $X_2\vert X_1,X_3$; and $X_3\vert X_1,X_2$.

A toy example

Suppose that you want to simulate from $(X_1,X_2)$ distributed as a bivariate normal with mean $(0,0)$ and covariance matrix $\begin{bmatrix}1 & \rho\\ \rho & 1\end{bmatrix}$ with $\rho=0.5$. It is known that the conditional distributions are given by

$$X_1\vert X_2 \sim \mbox{Normal}(\rho X_2,1-\rho^2),\\ X_2\vert X_1 \sim \mbox{Normal}(\rho X_1,1-\rho^2).$$

Note that these distributions are easy to simulate and then we can easily implement a Gibbs sampler. The following R code shows how to do this.

N = 26000 # Total number of iterations
rho = 0.5

simx = simy = vector()
simx[1] = simy[1] = 0 # Initialise vectors

# Gibbs sampler
for(i in 2:N){
simx[i] = rnorm(1,rho*simy[i-1],sqrt(1-rho^2))
simy[i] = rnorm(1,rho*simx[i],sqrt(1-rho^2))
}

# burnin and thinning
burnin = 1000
thinning = 25

# Gibbs sample
sim = cbind(simx[seq(from=burnin,to=N,by=thinning)],simy[seq(from=burnin,to=N,by=thinning)])

# A diagnosis tool to assess the convergence of the sampler
plot(ts(sim[,1]))
plot(ts(sim[,2]))

# Plots of the marginals and the joint distribution

hist(sim[,1])
hist(sim[,2])

plot(sim,col=1:1000)
plot(sim,type="l")

As you can see, there is no need for evaluating of the conditional densities or distributions. The only thing we need to do is to simulate from the corresponding conditionals. In your case you just need to simulate from the aforementioned conditionals and to conduct an analogous iterative algorithm.