control_variates_simulation.py

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(42)

def g(x):
    return np.log(1 + x) - np.log(x)

def f(x):
    return x

Ef = 0.5
n_samples = 100000
n_repeats = 10000

def monte_carlo(n):
    x = np.random.uniform(0, 1, n)
    return np.mean(g(x))

def control_variate(n):
    x = np.random.uniform(0, 1, n)
    gx = g(x)
    fx = f(x)
    cov = np.cov(gx, fx, ddof=1)[0, 1]
    var_f = np.var(fx, ddof=1)
    c = -cov / var_f
    return np.mean(gx + c * (fx - Ef))

mc_est = [monte_carlo(n_samples) for _ in range(n_repeats)]
cv_est = [control_variate(n_samples) for _ in range(n_repeats)]

std_mc = np.std(mc_est)
std_cv = np.std(cv_est)