import matplotlib.pyplot as plt
import numpy as np
from tqdm.auto import tqdm

7. Metropolis: Normal-Cauchy*

Adapted from Unit 5: norcaumet.m.

For this example:

\[\begin{split} \begin{align} X \mid \theta &\sim N(\theta, 1) \\ \theta &\sim \text{Cauchy}(0, 1) \end{align} \end{split}\]

Original method

The professor chooses a \(N(x, 1)\) proposal (where \(x=2\) is our sole datapoint). This is not symmetric around the current state of the chain, so it doesn’t satisfy \(q(\theta^\prime| \theta) = q(\theta | \theta^\prime)\) and we must use Metropolis-Hastings. It does simplify the acceptance ratio quite a bit. Here’s how we get the ratio:

\[\begin{split} \begin{align} \pi(\theta\mid X) &\propto f(x|\theta)\pi(\theta) \\ &\propto \frac{e^{-\frac{(x - \theta)^2}{2}}}{1+\theta^2} \end{align} \end{split}\]

\[\begin{split} \begin{align} \rho &= \frac{\pi(\theta')q(\theta\mid\theta')}{\pi(\theta)q(\theta'\mid\theta)} \\ &= \frac{\frac{e^{-\frac{(x - \theta')^2}{2}}}{1+{\theta'}^2}e^{-\frac{(\theta - x)^2}{2}}}{\frac{e^{-\frac{(x - \theta)^2}{2}}}{1+\theta^2}e^{-\frac{(\theta' - x)^2}{2}}} &= \frac{1 + \theta^2}{1 + {\theta'}^2} \end{align} \end{split}\]

Skipping some intermediate steps for now, but if you’d like to see them let me know and I’ll take some time to add them.

Warning

I usually pre-generate arrays of random numbers (see theta_prop and unif variables in the below cell) because it’s computationally faster. However, that only works when those numbers don’t depend on the previous step in the sampling loop.

rng = np.random.default_rng(1)

n = 1000000  # observations
burn = 500
theta = 1  # init
thetas = np.zeros(n)
x = 2  # observed

accepted = np.zeros(n)

# generating necessary randoms as arrays is faster
theta_prop = rng.standard_normal(n) + x
unif = rng.uniform(size=n)

for i in tqdm(range(n)):
    r = (1 + theta**2) / (1 + theta_prop[i] ** 2)
    rho = min(r, 1)
    if unif[i] < rho:
        theta = theta_prop[i]
        accepted[i] = 1
    thetas[i] = theta

thetas = thetas[burn:]

print(f"{np.mean(thetas)=:.3f}")
print(f"{np.var(thetas)=:.3f}")
print(f"{np.sum(accepted)/n=:.3f}")

fig, ax = plt.subplots()
ax.grid(True)
plt.hist(thetas, color="lightblue", density=True, bins=50)
plt.show()

np.mean(thetas)=1.283
np.var(thetas)=0.863
np.sum(accepted)/n=0.588

Random-walk proposal

How would we do this with classic Metropolis? We would use a random-walk proposal, centering the proposal at the previous state of the chain at each step.

Notice I can still pregenerate proposal values, because it’s easy to adjust the pregenerated \(N(0, 1)\) values by adding theta, which centers the distribution at the previous state of the chain.

n = 1000000
burn = 500
theta = 1
thetas = np.zeros(n)
x = 2

accepted = np.zeros(n)

theta_prop = rng.standard_normal(n)
unif = rng.uniform(size=n)


def f(x, theta):
    return np.exp(-0.5 * (x - theta) ** 2) / (1 + theta**2)


for i in tqdm(range(n)):
    theta_prop_current = theta_prop[i] + theta

    r = f(x, theta_prop_current) / f(x, theta)
    rho = min(r, 1)

    if unif[i] < rho:
        theta = theta_prop_current
        accepted[i] = 1

    thetas[i] = theta

thetas = thetas[burn:]

print(f"{np.mean(thetas)=:.3f}")
print(f"{np.var(thetas)=:.3f}")
print(f"{np.sum(accepted)/n=:.3f}")

fig, ax = plt.subplots()
ax.grid(True)
plt.hist(thetas, color="lightblue", density=True, bins=50)
plt.show()

np.mean(thetas)=1.283
np.var(thetas)=0.865
np.sum(accepted)/n=0.686

The end results are the same, with a higher acceptance ratio!

%load_ext watermark
%watermark -n -u -v -iv

Last updated: Sat Oct 14 2023

Python implementation: CPython
Python version       : 3.11.5
IPython version      : 8.15.0

numpy     : 1.25.2
matplotlib: 3.7.2