MCMC diagnostics
Explore the difference between MC and MCMC with a sample example.
A discrete variable
δ{1,2,3}\delta\sim\{1,2,3\}
and a continuous variable
θI ⁣R\theta\in\mathrm{I\!R}
.
1
# discrete variable
2
delta = c(.45, .10, .45)
3
4
# continuous variable
5
mu = c(-3, 0, 3)
6
sigma2 = c(1/3, 1/3, 1/3)
7
8
# exact marginal density of theta
9
ext_margin_den <- function(x)
10
{
11
dnorm(x, mu[1], sqrt(sigma2[1])) * delta[1] +
12
dnorm(x, mu[2], sqrt(sigma2[2])) * delta[2] +
13
dnorm(x, mu[3], sqrt(sigma2[3])) * delta[3]
14
}
15
16
theta = seq(-6, 6, length.out = 1000)
17
ptheta = ext_margin_den(theta)
18
plot(theta, ptheta, type = "l")
Copied!
The marginal density of
θ\theta
would be
Firstly, generate 1000 Monte Carlo
θ\theta
-samples.
1
sample_delta <- function(x)
2
{
3
denx = c(delta[1] * dnorm(x, mu[1], sqrt(sigma2[1])),
4
delta[2] * dnorm(x, mu[2], sqrt(sigma2[2])),
5
delta[3] * dnorm(x, mu[3], sqrt(sigma2[3])))
6
condenx = denx / sum(denx)
7
sample(1:3, 1, prob = condenx)
8
}
9
# gibbs sampler
10
gibbs <- function(m, delta = 3)
11
{
12
DELTA = delta
13
THETA = NULL
14
for (i in 1:m)
15
{
16
# sample theta from full conditional distribution
17
theta = sample_theta(delta)
18
THETA = c(THETA, theta)
19
delta = sample_delta(theta)
20
DELTA = c(DELTA, delta)
21
}
22
return(THETA)
23
}
24
res = gibbs(1000)
25
plot(res)
26
hist(res, breaks = 20)
Copied!
For 1000 MCMC samples, we have
Actually, re-run the above code, we can get much different figure. Let's try 10000 MCMC samples,
It turns out to be much stable when you re-run the above code.

sample autocorrelation

Use R-function acf. If a Markov chain with high autocorrelation, then it will move around the parameter space slowly, taking a long time to achieve the correct balance among the different regions of the parameter space.

effective sample size

Use R command effectiveSize in the coda package, which can be interpreted as the number of independent Monte Carlo samples necessary to give the same precision as the MCMC samples.
Last modified 3yr ago