未曾知全貌，我亦知何如（蒙特卡洛算法）

Target I = ∫_0^1 exp(x) dx = e - 1 = 1.718281828459045
Plain      : Ihat=1.72215757, SE=3.466e-03, 95% CI=(np.float64(1.7153641635016394), np.float64(1.7289509773653648))
Antithetic : Ihat=1.71772080, SE=6.245e-04, 95% CI=(np.float64(1.7164966886828101), np.float64(1.7189449125866085))
ControlVar : Ihat=1.71855180, SE=4.446e-04, 95% CI=(np.float64(1.7176804373321641), np.float64(1.7194231571928984)), b≈1.6884

import numpy as np
import matplotlib.pyplot as plt

# =========================
# Monte Carlo: mathematical modeling + simulation demo (fast version)
# =========================
# Target: I = ∫_0^1 exp(x) dx = e - 1 = E[exp(U)], U ~ Unif(0,1)
# Show:
#   - Estimator + 95% CI
#   - Convergence ~ 1/sqrt(N)
#   - Variance reduction: antithetic + control variate
#   - CLT intuition via standardized errors
# =========================

TRUE_I = np.e - 1.0

def mc_plain(N, rng):
    u = rng.random(N)
    y = np.exp(u)
    Ihat = y.mean()
    se = y.std(ddof=1) / np.sqrt(N)
    return Ihat, se

def mc_antithetic(N_total, rng):
    # Use N_total samples by pairing => N_pairs = N_total/2
    N_pairs = N_total // 2
    u = rng.random(N_pairs)
    z = 0.5 * (np.exp(u) + np.exp(1.0 - u))
    Ihat = z.mean()
    se = z.std(ddof=1) / np.sqrt(N_pairs)
    return Ihat, se

def mc_control(N, rng):
    u = rng.random(N)
    x = np.exp(u)
    c = u
    c0 = 0.5
    cov = np.cov(x, c, ddof=1)[0, 1]
    varc = np.var(c, ddof=1)
    b = cov / varc
    z = x - b * (c - c0)
    Ihat = z.mean()
    se = z.std(ddof=1) / np.sqrt(N)
    return Ihat, se, b

def ci95(Ihat, se):
    return (Ihat - 1.96*se, Ihat + 1.96*se)

# ---- Single run summary
rng = np.random.default_rng(0)
N = 20000
I1, se1 = mc_plain(N, rng)
I2, se2 = mc_antithetic(N, rng)
I3, se3, b3 = mc_control(N, rng)

print("Target I = ∫_0^1 exp(x) dx = e - 1 =", TRUE_I)
print(f"Plain      : Ihat={I1:.8f}, SE={se1:.3e}, 95% CI={ci95(I1,se1)}")
print(f"Antithetic : Ihat={I2:.8f}, SE={se2:.3e}, 95% CI={ci95(I2,se2)}")
print(f"ControlVar : Ihat={I3:.8f}, SE={se3:.3e}, 95% CI={ci95(I3,se3)}, b≈{b3:.4f}")

# ---- Convergence study (reduced)
Ns = np.unique(np.logspace(2, 5, 16).astype(int))  # 1e2..1e5, 16 points
n_rep = 60  # repetitions per N

def rms_error(method, N, seed_base=123):
    errs = np.empty(n_rep)
    for k in range(n_rep):
        rg = np.random.default_rng(seed_base + 10000*k + N)
        if method == "plain":
            Ihat, _ = mc_plain(N, rg)
        elif method == "antithetic":
            Ihat, _ = mc_antithetic(N, rg)
        elif method == "control":
            Ihat, _, _ = mc_control(N, rg)
        else:
            raise ValueError
        errs[k] = Ihat - TRUE_I
    return np.sqrt(np.mean(errs**2))

rms_plain = np.array([rms_error("plain", int(n)) for n in Ns])
rms_anti  = np.array([rms_error("antithetic", int(n)) for n in Ns])
rms_ctrl  = np.array([rms_error("control", int(n)) for n in Ns])

ref = rms_plain[0] * np.sqrt(Ns[0] / Ns)

plt.figure()
plt.loglog(Ns, rms_plain, marker="o", label="Plain MC")
plt.loglog(Ns, rms_anti,  marker="o", label="Antithetic")
plt.loglog(Ns, rms_ctrl,  marker="o", label="Control variate")
plt.loglog(Ns, ref, linestyle="--", label="~ 1/sqrt(N)")
plt.xlabel("N samples")
plt.ylabel("RMS error")
plt.title("Monte Carlo convergence (RMS error vs N)")
plt.grid(True, which="both")
plt.legend()
plt.show()

# ---- CLT intuition: standardized errors (reduced)
N_clt = 2000
reps = 1500
Ihats = np.empty(reps)
ses = np.empty(reps)
for k in range(reps):
    rg = np.random.default_rng(7 + k)
    Ihat, se = mc_plain(N_clt, rg)
    Ihats[k] = Ihat
    ses[k] = se
z = (Ihats - TRUE_I) / ses

plt.figure()
plt.hist(Ihats, bins=50)
plt.xlabel("Ihat")
plt.ylabel("Count")
plt.title(f"Sampling distribution of Ihat (plain), N={N_clt}, reps={reps}")
plt.grid(True)
plt.show()

plt.figure()
plt.hist(z, bins=50)
plt.xlabel("z = (Ihat - I)/SE")
plt.ylabel("Count")
plt.title("Standardized errors (approximately N(0,1) by CLT)")
plt.grid(True)
plt.show()