fisher_mm_variance.py

from functools import partial
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

@partial(np.vectorize, otypes=["float64"] * 2)
def mm_stds(rng, n_draws, alpha, beta, k=500):
    draws = np.log(rng.gamma(shape=beta, scale=alpha, size=(k, n_draws)))
    grad = beta - np.exp(draws) / alpha
    draws_var = draws.var(-1)
    grad_var = grad.var(-1)
    std_based = np.log(draws_var)
    fisher_based = np.log(np.sqrt(draws_var / grad_var))
    return std_based.std(), fisher_based.std()

rng = np.random.default_rng(42)
alpha = 1
beta = 10

n_draws = np.logspace(0.5, 3, 30, dtype=int)
stds, fishers = mm_stds(rng, n_draws, alpha=alpha, beta=beta)

plt.plot(n_draws, stds, label="stds")
plt.plot(n_draws, fishers, label="fisher")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Number of draws")
plt.ylabel("Std of mass matrix estimate")
plt.title(f"log-gamma(alpha={alpha}, beta={beta})")
plt.legend();