# QUESTION 1
# Slice sampling
rm(list = ls())
set.seed(5444)
# Given data
x1 <- 0.92
x2 <- 0.72
# x1 + x2
a <- x1 + x2
# Target:
# f(theta | x1, x2) proportional to
# theta^2 * exp(-(x1+x2)*theta), theta > 0
#
# The normalized target is Gamma(shape = 3, rate = 1.64)
# Number of retained samples
n.sample <- 1000
# Burn-in
burnin <- 200
# Total number of iterations
nMC <- n.sample + burnin
# Storage
theta.sample <- rep(NA, nMC)
# Initial value
theta <- 1
# Slice sampler
#
# Factor the target as
#
# g1(theta) = theta^2
# g2(theta) = exp(-a*theta)
#
# Introduce
# u1 | theta ~ Uniform(0, theta^2)
# u2 | theta ~ Uniform(0, exp(-a*theta))
#
# Conditional on u1 and u2:
#
# theta > sqrt(u1)
#
# and
#
# theta < -log(u2)/a
for (i in 1:nMC) {
# Step 1: sample auxiliary variables
u1 <- runif(
1,
min = 0,
max = theta^2
)
u2 <- runif(
1,
min = 0,
max = exp(-a * theta)
)
# Step 2: determine slice interval
lower <- sqrt(u1)
upper <- -log(u2) / a
# Step 3: sample theta from the slice
theta <- runif(
1,
min = lower,
max = upper
)
theta.sample[i] <- theta
}
# Removing burn-in
theta.keep <- theta.sample[(burnin + 1):nMC]
# histogram
hist(
theta.keep,
breaks = 35,
probability = TRUE,
col = "grey",
border = "white",
main = "Slice Sampling",
xlab = expression(theta)
)
# theoretical density
#
# theta | x1,x2 ~ Gamma(shape=3, rate=1.64)
theta.grid <- seq(
0,
max(theta.keep) * 1.05,
length.out = 500
)
f.theta <- dgamma(
theta.grid,
shape = 3,
rate = a
)
lines(
theta.grid,
f.theta,
lwd = 2
)
# numerical summaries
mean(theta.keep)
## [1] 1.723104
var(theta.keep)
## [1] 1.068989
range(theta.keep)
## [1] 0.1134067 6.4857268
# QUESTION 2
# Parallel Tempering
rm(list = ls())
set.seed(5444)
# Mixture distribution parameters
mu1 <- c(0, 0)
mu2 <- c(15, 0)
mu3 <- c(9, 1.5)
# Variances
# The covariance matrices are:
# 0.10*I, 0.08*I, and 0.09*I
sig1 <- 0.10
sig2 <- 0.08
sig3 <- 0.09
# Mixture weights
w1 <- 0.3
w2 <- 0.3
w3 <- 0.4
# Bivariate normal density
#
# If X ~ N(mu, sig*I),
#
# f(x) = 1/(2*pi*sig) *
# exp{-||x-mu||^2/(2*sig)}
#
# Here sig is the variance, not the standard deviation.
dmvnorm.iso <- function(x, mu, sig) {
f <- 1 / (2 * pi * sig) *
exp(
-sum((x - mu)^2) / (2 * sig)
)
return(f)
}
# Target mixture density
#
# f(x1,x2)
# =
# 0.3 N(mu1,0.10I)
# + 0.3 N(mu2,0.08I)
# + 0.4 N(mu3,0.09I)
poster <- function(x) {
f1 <- w1 *
dmvnorm.iso(
x,
mu = mu1,
sig = sig1
)
f2 <- w2 *
dmvnorm.iso(
x,
mu = mu2,
sig = sig2
)
f3 <- w3 *
dmvnorm.iso(
x,
mu = mu3,
sig = sig3
)
f <- f1 + f2 + f3
return(f)
}
# Log target density
log.poster <- function(x) {
f <- poster(x)
# Prevent taking log(0)
f <- max(
f,
.Machine$double.xmin
)
return(
log(f)
)
}
# PARALLEL TEMPERING FUNCTION
parallel.tempering <- function(
Tset,
x0,
alpha0,
nMC,
delta
) {
# Number of temperature chains
I <- length(Tset)
# Dimension of x
d <- ncol(x0)
# Initial states
x.current <- x0
# Current log target densities
logf.current <- rep(
NA,
I
)
for (k in 1:I) {
logf.current[k] <-
log.poster(
x.current[k, ]
)
}
# Storage for all chains
sample.array <- array(
NA,
dim = c(
I,
d,
nMC
)
)
# Acceptance records
MH.accept <- rep(
0,
I
)
MH.total <- rep(
0,
I
)
swap.accept <- rep(
0,
I - 1
)
swap.total <- rep(
0,
I - 1
)
# Main parallel tempering iterations
for (m in 1:nMC) {
# Draw u to determine whether to conduct
# a parallel step or a swapping step
u <- runif(1)
# STEP 1: PARALLEL MCMC UPDATE
#
# With probability alpha0, update every chain.
if (u < alpha0) {
for (k in 1:I) {
MH.total[k] <-
MH.total[k] + 1
# Random-walk proposal
x.star <-
x.current[k, ] +
rnorm(
d,
mean = 0,
sd = delta
)
# Log density at proposed value
logf.star <-
log.poster(
x.star
)
# Tempered target:
#
# pi_k(x) proportional to pi(x)^(1/T_k)
#
# Therefore,
#
# log alpha =
# [log pi(x*) - log pi(x)] / T_k
log.alpha <-
(
logf.star -
logf.current[k]
) /
Tset[k]
# Accept / reject
if (
log(runif(1)) <
min(
0,
log.alpha
)
) {
x.current[k, ] <-
x.star
logf.current[k] <-
logf.star
MH.accept[k] <-
MH.accept[k] + 1
}
}
# STEP 2: SWAPPING STEP
#
# With probability 1-alpha0, choose neighboring chains
# and propose a swap.
} else {
# Choose one neighboring pair
j <- sample(
1:(I - 1),
size = 1
)
swap.total[j] <-
swap.total[j] + 1
# Swap acceptance probability
#
# log alpha_swap
#
# =
#
# (1/T_j - 1/T_{j+1})
#
# *
#
# [log pi(x_{j+1}) - log pi(x_j)]
log.alpha.swap <-
(
1 / Tset[j] -
1 / Tset[j + 1]
) *
(
logf.current[j + 1] -
logf.current[j]
)
# Accept / reject swap
if (
log(runif(1)) <
min(
0,
log.alpha.swap
)
) {
# Swap current states
temp.x <-
x.current[j, ]
x.current[j, ] <-
x.current[j + 1, ]
x.current[j + 1, ] <-
temp.x
# Swap corresponding log densities
temp.f <-
logf.current[j]
logf.current[j] <-
logf.current[j + 1]
logf.current[j + 1] <-
temp.f
swap.accept[j] <-
swap.accept[j] + 1
}
}
# Store the state of all chains
sample.array[, , m] <-
x.current
}
# Acceptance rates
rMH <-
MH.accept /
MH.total
rSW <-
swap.accept /
swap.total
# Return results
return(
list(
s.sample = sample.array,
rMH = rMH,
rSW = rSW
)
)
}
# RUN PARALLEL TEMPERING
# Temperatures
Tset <- c(
100,
3,
1,
0.5
)
I <- length(Tset)
# Initial states
x0 <- cbind(
runif(
I,
min = -5,
max = 15
),
runif(
I,
min = -10,
max = 10
)
)
# Probability of conducting parallel step
alpha0 <- 0.5
# Number of MCMC iterations
nMC <- 80000
# Burn-in
burnin <- 10000
# Random-walk proposal standard deviation
delta <- 0.5
# Run algorithm
Mod <- parallel.tempering(
Tset = Tset,
x0 = x0,
alpha0 = alpha0,
nMC = nMC,
delta = delta
)
# Extract samples
s.sample <-
Mod$s.sample
# Print MH acceptance rates
Mod$rMH
## [1] 0.9252465 0.5642925 0.3568997 0.2369238
# Print swap acceptance rates
Mod$rSW
## [1] 0.06902219 0.49809003 0.65867408
# FUNCTION TO PLOT TEMPERED DENSITY + SAMPLE POINTS
plot.temperature <- function(
temp,
samples,
main.title
) {
# Grid
x1.grid <- seq(
-5,
20,
length.out = 160
)
x2.grid <- seq(
-8,
8,
length.out = 130
)
den <- matrix(
NA,
nrow = length(x1.grid),
ncol = length(x2.grid)
)
# Evaluate tempered density on grid
for (i in 1:length(x1.grid)) {
for (j in 1:length(x2.grid)) {
xtemp <- c(
x1.grid[i],
x2.grid[j]
)
# Tempered target:
#
# pi_T(x) proportional to pi(x)^(1/T)
den[i, j] <-
exp(
log.poster(xtemp) /
temp
)
}
}
# Contour plot
contour(
x1.grid,
x2.grid,
den,
nlevels = 12,
drawlabels = FALSE,
xlab = expression(x[1]),
ylab = expression(x[2]),
main = main.title,
cex.main = 0.9
)
# Plot a thinned version of sample points
use <- seq(
1,
nrow(samples),
by = 20
)
points(
samples[use, 1],
samples[use, 2],
pch = 20,
cex = 0.35
)
}
# PLOT TWO DIFFERENT TEMPERATURES
#
# T = 100: hot chain
# T = 1: original target distribution
par(
mfrow = c(1, 2),
mar = c(5, 4, 3, 1)
)
# Temperature T = 100
id.hot <- which(
Tset == 100
)
sample.hot <- t(
s.sample[
id.hot,
,
(burnin + 1):nMC
]
)
plot.temperature(
temp = 100,
samples = sample.hot,
main.title = "T = 100"
)
# Temperature T = 1
id.target <- which(
Tset == 1
)
sample.target <- t(
s.sample[
id.target,
,
(burnin + 1):nMC
]
)
plot.temperature(
temp = 1,
samples = sample.target,
main.title = "T = 1"
)
# Return plot layout to normal
par(
mfrow = c(1, 1)
)
# QUESTION 3
# Stochastic volatility model
rm(list = ls())
set.seed(100)
# Parameters
T <- 100
phi <- 0.9
sig2 <- 1
sigma <- sqrt(sig2)
gamma <- 1
# Storage
x <- rep(NA, T)
y <- rep(NA, T)
# QUESTION 3(b)
# Simulate latent states and observations
# Initial state
x[1] <- rnorm(
1,
mean = 0,
sd = sqrt(sig2 / (1 - phi^2))
)
# Latent process
for (t in 2:T) {
x[t] <-
phi * x[t - 1] +
rnorm(
1,
mean = 0,
sd = sigma
)
}
# Observation process
# w_t ~ N(0, 0.25), so sd(w_t) = 0.5
w <- rnorm(
T,
mean = 0,
sd = 0.5
)
y <- exp(gamma + x) * w
# FIGURE 3
# Latent state and observed process
par(
mfrow = c(2, 1),
mar = c(4, 4, 3, 1)
)
# Latent state
plot(
1:T,
x,
type = "l",
lwd = 2,
xlab = "Time",
ylab = expression(x[t]),
main = "Latent State Process"
)
# Observed process
plot(
1:T,
y,
type = "l",
lwd = 1,
xlab = "Time",
ylab = expression(y[t]),
main = "Observed Process"
)
points(
1:T,
y,
pch = 16,
cex = 0.4
)
par(mfrow = c(1, 1))
# QUESTION 3(c)
# Sequential Importance Resampling
N <- 1000
# Particle storage
Xs <- matrix(
NA,
nrow = N,
ncol = T
)
# Posterior summaries
post.mean <- rep(NA, T)
post.lower <- rep(NA, T)
post.upper <- rep(NA, T)
# SIR algorithm
for (t in 1:T) {
# Step 1: Propagate
if (t == 1) {
Xs[, t] <- rnorm(
N,
mean = 0,
sd = sqrt(sig2 / (1 - phi^2))
)
} else {
Xs[, t] <- rnorm(
N,
mean = phi * Xs[, t - 1],
sd = sigma
)
}
# Step 2: Importance weights
# y_t | x_t ~ N(0, 0.25*exp(2*(gamma+x_t)))
# Conditional standard deviation:
# 0.5*exp(gamma+x_t)
obs.sd <-
0.5 *
exp(
gamma + Xs[, t]
)
logweight <- dnorm(
y[t],
mean = 0,
sd = obs.sd,
log = TRUE
)
# Stable normalization
logweight <-
logweight -
max(logweight)
weight <- exp(logweight)
weight <-
weight /
sum(weight)
# Step 3: Resample
id <- sample(
1:N,
size = N,
replace = TRUE,
prob = weight
)
Xs[, t] <-
Xs[
id,
t
]
# Step 4: Posterior summaries
post.mean[t] <-
mean(
Xs[, t]
)
post.lower[t] <-
quantile(
Xs[, t],
probs = 0.025
)
post.upper[t] <-
quantile(
Xs[, t],
probs = 0.975
)
}
# Numerical summaries
range(x)
## [1] -4.142705 3.130338
range(y)
## [1] -66.44380 13.63763
sqrt(
mean(
(post.mean - x)^2
)
)
## [1] 0.8027361
mean(
x >= post.lower &
x <= post.upper
)
## [1] 0.95
# QUESTION 3(d)
# FIGURE 4
# Filtering results
layout(
matrix(
c(1, 2, 3),
ncol = 1
),
heights = c(4, 1, 4)
)
# Panel 1: Latent state and filtering results
par(
mar = c(4, 4, 3, 1)
)
ylim.state <- range(
c(
x,
post.mean,
post.lower,
post.upper
)
)
plot(
1:T,
x,
type = "l",
lwd = 2,
lty = 1,
ylim = ylim.state,
xlab = "Time",
ylab = "State",
main = "Filtering Results for the Latent State"
)
# Posterior mean
lines(
1:T,
post.mean,
lwd = 2,
lty = 2
)
# Lower credible bound
lines(
1:T,
post.lower,
lwd = 1,
lty = 3
)
# Upper credible bound
lines(
1:T,
post.upper,
lwd = 1,
lty = 3
)
# Panel 2: Legend
par(
mar = c(0, 0, 0, 0)
)
plot.new()
legend(
"center",
legend = c(
"latent x",
"posterior mean",
"95% credible interval"
),
lty = c(
1,
2,
3
),
lwd = c(
2,
2,
1
),
horiz = TRUE,
bty = "n",
cex = 0.9
)
# Panel 3: Observed data
par(
mar = c(4, 4, 3, 1)
)
plot(
1:T,
y,
type = "l",
lwd = 1,
xlab = "Time",
ylab = "Observation",
main = "Observed Data"
)
points(
1:T,
y,
pch = 16,
cex = 0.35
)
# Reset layout
layout(matrix(1))
par(mfrow = c(1, 1))