Gibbs

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

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

Then he argues that

$m$

copies of $\mathbf y_{mis}$

in each iteration is not really necessary. And briefly summary the DA algorithm:Another example:

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.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.

Last modified 3yr ago