Gibbs

Search…
Gibbs
​Robert and Casella (2013) gives the following algorithm,
while Liu, Jun S. (2008) introduces two types of Gibbs sampling strategy.
Bivariate Gibbs sampler
It is easy to implement this sampler:
1
## Julia program for Bivariate Gibbs sampler
2
## author: weiya <[email protected]>
3
## date: 2018-08-22
4
​
5
function bigibbs(T, rho)
6
    x = ones(T+1)
7
    y = ones(T+1)
8
    for t = 1:T
9
        x[t+1] = randn() * sqrt(1-rho^2) + rho*y[t]
10
        y[t+1] = randn() * sqrt(1-rho^2) + rho*x[t+1]
11
    end
12
    return x, y
13
end
14
​
15
## example
16
bigibbs(100, 0.5)
Copied!
Completion Gibbs Sampler
Example:
We can use the following Julia program to implement this algorithm.
1
## Julia program for Truncated normal distribution
2
## author: weiya <[email protected]>
3
## date: 2018-08-22
4
​
5
# Truncated normal distribution
6
function rtrunormal(T, mu, sigma, mu_down)
7
    x = ones(T)
8
    z = ones(T+1)
9
    # set initial value of z
10
    z[1] = rand()
11
    if mu < mu_down
12
        z[1] = z[1] * exp(-0.5 * (mu - mu_down)^2 / sigma^2)
13
    end
14
    for t = 1:T
15
        x[t] = rand() * (mu - mu_down + sqrt(-2*sigma^2*log(z[t]))) + mu_down
16
        z[t+1] = rand() * exp(-(x[t] - mu)^2/(2*sigma^2))
17
    end
18
    return(x)
19
end
20
​
21
## example
22
rtrunormal(1000, 1.0, 1.0, 1.2)
Copied!
Slice Sampler
​Robert and Casella (2013) introduces the following slice sampler algorithm,
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.
1
## Julia program for Slice sampler
2
## author: weiya <[email protected]>
3
## date: 2018-08-22
4
​
5
function rnorm_slice(T)
6
    x = ones(T+1)
7
    w = ones(T+1)
8
    for t = 1:T
9
        w[t+1] = rand() * exp(-1.0 * x[t]^2/2)
10
        x[t+1] = rand() * 2 * sqrt(-2*log(w[t+1])) - sqrt(-2*log(w[t+1]))
11
    end
12
    return x[2:end]
13
end
14
​
15
## example
16
rnorm_slice(100)
Copied!
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 mmm
 copies of ymis\mathbf y_{mis}ymis​
 in each iteration is not really necessary. And briefly summary the DA algorithm:
It seems to agree with the algorithm presented by Robert and Casella (2013).
Another example:
It seems that we do not need to derive the explicit form of g(x,z)g(x, z)g(x,z)
, if we can directly obtain the conditional distribution. We can use the following Julia program to sample.
1
## Julia program for Grouped Multinomial Data (Ex. 7.2.3)
2
## author: weiya <[email protected]>
3
## date: 2018-08-26
4
​
5
# call gamma function
6
#using SpecialFunctions
7
​
8
​
9
# sample from Dirichlet distributions
10
using Distributions
11
​
12
function gmulti(T, x, a, b, alpha1 = 0.5, alpha2 = 0.5, alpha3 = 0.5)
13
    z = ones(T+1, size(x, 1)-1) # initial z satisfy `z <= x`
14
    mu = ones(T+1)
15
    eta = ones(T+1)
16
    for t = 1:T
17
        # sample from g_1(theta | y)
18
        dir = Dirichlet([z[t, 1] + z[t, 2] + alpha1, z[t, 3] + z[t, 4] + alpha2, x[5] + alpha3])
19
        sample = rand(dir, 1)
20
        mu[t+1] = sample[1]
21
        eta[t+1] = sample[2]
22
        # sample from g_2(z | x, theta)
23
        for i = 1:2
24
            bi = Binomial(x[i], a[i]*mu[t+1]/(a[i]*mu[t+1]+b[i]))
25
            z[t+1, i] = rand(bi, 1)[1]
26
        end 
27
        for i = 3:4
28
            bi = Binomial(x[i], a[i]*eta[t+1]/(a[i]*eta[t+1]+b[i]))
29
            z[t+1, i] = rand(bi, 1)[1]
30
        end
31
    end
32
    return mu, eta
33
end
34
​
35
# example
36
​
37
## data
38
a = [0.06, 0.14, 0.11, 0.09];
39
b = [0.17, 0.24, 0.19, 0.20];
40
x = [9, 15, 12, 7, 8]; 
41
​
42
gmulti(100, x, a, b)
Copied!
Reversible Data Augmentation
Reversible Gibbs Sampler
Random Sweep Gibbs Sampler
Random Gibbs Sampler
Hybrid Gibbs Samplers
Metropolization of the Gibbs Sampler
Let us illuatrate this algorithm with the following example.