Columbia University
IEOR4703 – Monte Carlo Simulation (Hirsa)
Assignment 8 – Due on Tuesday May 2nd, 2023 by midnight
Problem 1 (Hamiltonian Monte-Carlo): As you recall Gibbs Sampler failed to sample from the following distribution:
f(x1, x2) =
{ 1
2 : 0 < x1 < 1 & 0 < x2 < 1
1
2 : −1 < x1 < 0 & −1 < x2 < 0
We can write p(x) ∝ exp(0)1{x∈A} where A = ([0, 1]× [0, 1])∪ ([−1, 0]× [−1, 0]). For Hamiltonian Monte Carlo, introduce the momentum
variable y and we assume that it is normal, so that p(y) ∝ exp(− 12 (y21 + y22)). The joint density is then
p(x,y) =
1
Z
exp
(
0− 1
2
(y21 + y
2
2)
)
1{x∈A},
where Z is a constant. Form here, write Hamiltonian equations, use leapfrog discretization and the rest should be clear.
Problem 2 (Hamiltonian vs. Metropolis-Hastings): Consider a 10-dimensional multivariate normal distribution N (θ;µ,Σ) its
pdf is given by
pi(θ) ∝ 1√|Σ|e− 12 (θ−µ)>Σ−1(θ−µ)
where θ = (θ1, . . . θ10). Assume µ = (0, . . . , 0) and
Λ =

1 . . . 0 0
...
. . .
...
...
0 . . . 1 ρ
0 . . . ρ 1
 σ =

0.1 0 . . . 0 0
0 0.2 . . . 0 0
...
...
. . .
...
...
0 0 . . . 0.9 0
0 0 . . . . . . 1

where
Σ = σ>Λ σ
• In Metropolis-Hastings, assume the following distribution for proposal distribution, q(θ|θ′),
θ′ = θ + ηZ Z ∼ N (0, I)
Find Metropolis ratio and use Metropolis-Hasting algorithm to sample from the 10-dimensional multivariate normal
• For Hamiltonian Monte Carlo, using the following kinetic energy1
K(p) = p>p/2
to find Hamiltonian equations and sample from the 10-dimensional multivariate normal.
Plot & Compare traceplots and running averages and conclude on each method’s performance. Use few different values for ρ, e.g. 0.1,
0.5, 0.99, -0.1, -0.5, -0.99.
1where p = (p1, . . . , pd)