Assigned 2023-11-14, due 2023-11-21 on Gradescope
There is no written part to this assignment. This assignment consists of three coding parts (three separate notebooks). Please submit the notebooks with the output included on Gradescope.
Download this notebook, and follow the instructions. If you want to see what the results should look like, you can find it here.
Download this notebook, and follow the instructions. If you want to see what the results should look like, you can find it here. I used 5 books from Project Guttenberg as a data sample. Please download a few of your own (one long book should suffice; download it as a plain text file).
There is no template for this part; it is a modification of code you wrote in homework 8.
Recall the distribution from Assignment 8: Let $\mu_1 = (1, 0)$, $\mu_2 = (-1, 0)$, $\kappa_1 = 1/100$, $\kappa_2 = 1/10$. Let $G_i$ be the PDF of a $2$-dimensional Gaussian with mean $\mu_i$ and diagonal covariance matrix $\kappa_i I$, and define \begin{equation}\label{e:stationary} p = 0.9 G_1 + 0.1 G_2 \,, \quad U = -\ln p\,. \end{equation}
The code you created in Assignment 8 simulated the Langevin Monte Carlo algorithm to sample from this distribution. On this homework, we will modify that code to simulate the Metropolis Adjusted Langevin Monte Carlo Algorithm (MALA).
Write a function to implement MALA. That is, time discretize the (two dimensional) SDE $dX = -\grad U(X) \, dt + \sqrt{2} \, dW_t$ (where $W$ is a two dimensional Brownian motion). But instead of updating $X_n$ using Euler Maruyama, we view the point Euler Maruyama gives us as a proposal and accept/reject it by Metropolis. (Look in the class notes for a detailed description and exact formula.) Also compute the average number of points which are to the left of the $y$-axis. Here’s roughly what the function should look like, and what it should return.
def mala( X0, n_iters, dt=.001 ):
# X0 is the initial distribution. Should be a 2 x 1000 array.
# n_iters is the number of iterations, and dt is the time step.
X = X0.copy()
on_left = []
for _ in trange( n_iters ):
# Find the new X by MALA
X = ...
on_left.append( average( X[0, :] < 0 ) )
# Return final distribution, and the on_left array
return X, on_left
You already generated $1,000$ samples from the distribution $p$ in a $2 \times 1000$ array $Y$ in homework 8.
Use this as your initial data for MALA with a time step $dt = 0.001$, n_iters = int(10/dt) and generate the same plot you did in assignment 8.
(We now know that $Y$ is the stationary distribution – so the final distribution should look the same and the number on the left should be roughly $10\%$.
This is mainly to test if your code works.)
One advantage of MALA is that it still gives reasonable results with a larger timestep.
Try the previous part with $dt = 0.1$, n_iters = int(10/dt) and generate the same plots.
(With LMC you get garbage; with MALA they should look OK.)
If you’re running locally, and have the computational power, try running MALA for a really long time (and say uniform initial data) to get close to the stationary distribution. (Don’t turn this in with your homework.) The example given here is a worst case scenario of LMC / MALA, and both algorithms perform similarly, and will take an exceptionally long time to converge to the stationary distribution. However, in this case I can take a larger time step with MALA, and so I will need less iterations to get close to the stationary distribution.