{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "BwUwrL3Vs6YA"
},
"source": [
"## Metropolis Hasting on a Mixture of Gaussians in 1d\n",
"$\n",
"\\newcommand{\\PP}{\\mathbb P}\n",
"\\newcommand{\\PE}{\\mathbb E}\n",
"\\newcommand{\\Xset}{\\mathsf{X}}\n",
"\\newcommand{\\nset}{\\mathbb{N}}\n",
"\\newcommand{\\invcdf}[1]{F_{#1}^{\\leftarrow}}\n",
"\\newcommand{\\rmd}{\\mathrm{d}}\n",
"\\newcommand{\\rme}{\\mathrm{e}}\n",
"$\n",
"This jupyter notebook is based on Randal Douc's google colab whose link can be found at the beginning of the course notes: https://wiki.randaldouc.xyz/lib/exe/fetch.php?media=world:polymcmc.pdf.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"id": "Yym8_tdfT3nu"
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import numpy.random as npr\n",
"import scipy.stats as stats\n",
"from scipy.stats import expon, geom, norm\n",
"from math import pi\n",
"from pylab import *\n",
"%matplotlib inline\n",
"\n",
"from IPython.display import HTML\n",
"from matplotlib import animation, rc"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "VUCIQn9ETl-R"
},
"source": [
"## Symmetric Random Walk Metropolis Hasting algorithm\n",
"\n",
"We now consider a target distribution which is the mixing of two gaussian distributions, one centered at $a$ and the other one centered at $-a$\n",
"$$\n",
"\\pi(x)=\\frac{1}{2}\\left(\\phi(x-a)+\\phi(x+a)\\right)=\\frac{1}{2} \\frac{\\rme^{-(x-a)^2}}{\\sqrt{2\\pi}}+\\frac{1}{2} \\frac{\\rme^{-(x+a)^2}}{\\sqrt{2\\pi}}\n",
"$$\n",
"where $\\phi$ is the density of the centered standard normal distribution.\n",
"\n",
"To target this distribution, we sample according to a Symmetric Random Walk Metropolis Hasting algorithm. When the chain is at the state $X_k$, we propose a candidate $Y_{k+1}$ according to $Y_{k+1}=X_k+Z_k$ where $Z_k\\sim {\\mathcal N}(0,1)$ and then we accept $X_{k+1}=Y_{k+1}$ with probability $\\alpha(X_k,Y_{k+1})$, where $\\alpha(x,y)=\\frac{\\pi(y)}{\\pi(x)} \\wedge 1$. Otherwise, $X_{k+1}=X_{k}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"colab": {
"background_save": true
},
"id": "xDS-IcTY0P7-",
"outputId": "94bbfa65-9bbe-40b0-c0c6-7ac503dbfbaf"
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"