import numpy as np
import matplotlib.pyplot as plt

def rejection_sampler(target_pdf, proposal_sampler, proposal_pdf, M, n_samples, batch=10000):
    np.random.seed(42)  # 设置随机种子以确保结果可重复
    samples = []
    total_draws = 0
    while len(samples) < n_samples:
        x = proposal_sampler(batch)
        u = np.random.rand(batch)
        # 防止 proposal_pdf 为 0 导致除零
        denom = M * proposal_pdf(x)
        # 若 denom 为 0，则 accept_prob 视为 0（因为 proposal 在该点没有质量）
        with np.errstate(divide='ignore', invalid='ignore'):
            accept_prob = np.where(denom>0, target_pdf(x) / denom, 0.0)
        # 把数值误差之外的 accept_prob 截断到 [0,1]
        accept_prob = np.clip(accept_prob, 0.0, 1.0)
        accept = u < accept_prob
        accepted = x[accept]
        samples.extend(accepted.tolist())
        total_draws += batch
    samples = np.array(samples[:n_samples])
    acceptance_empirical = n_samples / total_draws
    return samples, acceptance_empirical, total_draws

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta, norm

np.random.seed(2025)

def plot_rejection_geometry_beta(n_points=5000):
    a, b = 2, 5
    xs = np.linspace(0,1,400)
    f_vals = beta.pdf(xs, a, b)
    q_vals = np.ones_like(xs)
    M = f_vals.max()
    X = np.random.rand(n_points)
    Y = np.random.rand(n_points) * M
    accept_mask = Y <= beta.pdf(X, a, b)
    accepted = np.sum(accept_mask)
    empirical_accept = accepted / n_points
    theoretical_accept = 1.0 / M
    
    plt.figure(figsize=(8,5))
    plt.scatter(X[~accept_mask], Y[~accept_mask], s=6, alpha=0.6, label='rejected')
    plt.scatter(X[accept_mask], Y[accept_mask], s=6, alpha=0.6, label='accepted')
    plt.plot(xs, f_vals, linewidth=2, label='target f(x)=Beta(2,5)')
    plt.plot(xs, M*q_vals, linewidth=1.5, linestyle='--', label='M q(x) (envelope)')
    plt.text(0, -0.8, 'Beta(2,5) 的目标密度 f(x)（实线）和包络 Mq(x)（虚线），\nq 为 Uniform(0,1)）以及在区域 R={0≤y≤Mq(x)} 内的随机点。\n橙色点表示被接受（落在 y≤f(x)），蓝色点表示被拒绝。\n图上也显示了经验接受率与理论接受率。', fontsize=12, color='red')
    plt.title(f'Beta(2,5) rejection geometry: empirical accept={empirical_accept:.3f}, theoretical={theoretical_accept:.3f}')
    plt.xlabel('x')
    plt.ylabel('y (height)')
    plt.ylim(-0.02, M*1.05)
    plt.legend()
    plt.show()

def plot_rejection_geometry_normal(n_points=8000, sigma=2.0):
    xs = np.linspace(-4,4,600)
    f_vals = norm.pdf(xs, loc=0, scale=1.0)
    q_vals = norm.pdf(xs, loc=0, scale=sigma)
    ratio = f_vals / q_vals
    M = ratio.max()
    X = np.random.normal(loc=0.0, scale=sigma, size=n_points)
    qX = norm.pdf(X, loc=0.0, scale=sigma)
    Y = np.random.rand(n_points) * (M * qX)
    accept_mask = Y <= norm.pdf(X, loc=0.0, scale=1.0)
    accepted = np.sum(accept_mask)
    empirical_accept = accepted / n_points
    theoretical_accept = 1.0 / M
    
    mask_plot = (X >= -4) & (X <= 4)
    plt.figure(figsize=(8,5))
    plt.scatter(X[mask_plot & ~accept_mask], Y[mask_plot & ~accept_mask], s=6, alpha=0.5, label='rejected')
    plt.scatter(X[mask_plot & accept_mask], Y[mask_plot & accept_mask], s=6, alpha=0.6, label='accepted')
    plt.plot(xs, f_vals, linewidth=2, label='target f(x)=N(0,1)')
    plt.plot(xs, M*q_vals, linewidth=1.5, linestyle='--', label=f'M q(x), M≈{M:.3f}')
    plt.text(-4, -0.11, 'Normal(0,1) 目标与 Normal(0,2²) 提议的对应可视化（仅显示 x∈[−4,4]），\n同样以点的上/下位置区分接受与拒绝，并标出理论/经验接受率。', fontsize=12, color='red')
    plt.title(f'Normal(0,1) vs Normal(0,{sigma}^2) proposal: empirical accept={empirical_accept:.3f}, theoretical={theoretical_accept:.3f}')
    plt.xlabel('x')
    plt.ylabel('y (height)')
    plt.ylim(-0.01, (M * q_vals).max()*1.05)
    plt.xlim(-4,4)
    plt.legend()
    plt.show()

plot_rejection_geometry_beta(n_points=6000)
plot_rejection_geometry_normal(n_points=12000, sigma=2.0)

def summary():
    a,b = 2,5
    xs = np.linspace(0,1,10001)
    M_beta = beta.pdf(xs,a,b).max()
    sigma=2.0
    xs2 = np.linspace(-10,10,20001)
    M_normal = (norm.pdf(xs2,0,1) / norm.pdf(xs2,0,sigma)).max()
    return {"beta": {"M":float(M_beta), "theoretical_accept": float(1.0/M_beta)},
            "normal": {"sigma":sigma, "M":float(M_normal), "theoretical_accept": float(1.0/M_normal)}}

summary()

# ---------- Beta(2,5) target, Uniform(0,1) proposal ----------
def beta25_pdf(x):
    x = np.asarray(x)
    pdf = np.zeros_like(x, dtype=float)
    mask = (x>=0) & (x<=1)
    pdf[mask] = 30.0 * x[mask] * (1.0 - x[mask])**4
    return pdf

def uniform01_sampler(n):
    return np.random.rand(n)

def uniform01_pdf(x):
    x = np.asarray(x)
    pdf = np.zeros_like(x, dtype=float)
    mask = (x>=0) & (x<=1)
    pdf[mask] = 1.0
    return pdf

# 用细网格数值求最大值（严格做法通常需要解析，但网格足够细即可）
xs = np.linspace(0,1,10001)
p_vals = beta25_pdf(xs)
p_max = p_vals.max()
M_beta = p_max  # q(x)=1, 因此 M = max_x p(x)
theoretical_accept_beta = 1.0 / M_beta

samples_beta, acc_emp_beta, draws_beta = rejection_sampler(beta25_pdf, uniform01_sampler, uniform01_pdf, M_beta, 5000)

plt.figure(figsize=(8,4))
plt.hist(samples_beta, bins=50, density=True)
xx = np.linspace(0,1,300)
plt.plot(xx, beta25_pdf(xx))
plt.title(f'Beta(2,5) target with Uniform(0,1) proposal\nM={M_beta:.6f}, theoretical accept={theoretical_accept_beta:.6f}, empirical accept~{acc_emp_beta:.6f}')
plt.xlabel('x')
plt.ylabel('density')
plt.show()

print('Beta example: M={:.6f}, theoretical_accept={:.6f}, empirical_accept~{:.6f}, draws_used={:d}'.format(M_beta, theoretical_accept_beta, acc_emp_beta, draws_beta))

# ---------- 标准正态 target N(0,1), proposal N(0,σ^2) ----------
sigma = 2.0

def normal_pdf(x, mu=0.0, sigma=1.0):
    x = np.asarray(x)
    coef = 1.0 / np.sqrt(2.0 * np.pi * sigma**2)
    return coef * np.exp(-0.5 * ((x - mu) / sigma)**2)

def normal_sampler(n, mu=0.0, sigma=1.0):
    return np.random.normal(loc=mu, scale=sigma, size=n)

def normal_proposal_pdf(x):
    return normal_pdf(x, mu=0.0, sigma=sigma)

# 数值求 M（理论上若同均值且 sigma>1, M = sigma）
xs = np.linspace(-10,10,20001)
ratio = normal_pdf(xs, mu=0.0, sigma=1.0) / normal_proposal_pdf(xs)
M_normal = ratio.max()
theoretical_accept_normal = 1.0 / M_normal

# 运行采样（注意 target_pdf 的参数要正确）
samples_norm, acc_emp_norm, draws_norm = rejection_sampler(lambda x: normal_pdf(x, mu=0.0, sigma=1.0),
                                                           lambda n: normal_sampler(n, mu=0.0, sigma=sigma),
                                                           normal_proposal_pdf,
                                                           M_normal,
                                                           5000)

plt.figure(figsize=(8,4))
plt.hist(samples_norm, bins=50, density=True)
xx = np.linspace(-5,5,400)
plt.plot(xx, normal_pdf(xx, mu=0.0, sigma=1.0))
plt.title(f'Normal(0,1) target with Normal(0,{sigma}^2) proposal\nM={M_normal:.6f}, theoretical accept={theoretical_accept_normal:.6f}, empirical accept~{acc_emp_norm:.6f}')
plt.xlabel('x')
plt.ylabel('density')
plt.show()

print('Normal example: sigma={:.2f}, M={:.6f}, theoretical_accept={:.6f}, empirical_accept~{:.6f}, draws_used={:d}'.format(sigma, M_normal, theoretical_accept_normal, acc_emp_norm, draws_norm))

# 重新执行（内核已重置）
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(2025)

def rejection_sampling_pi(N=250000):
    x = np.random.uniform(-1, 1, size=N)
    y = np.random.uniform(-1, 1, size=N)
    inside = (x*x + y*y) <= 1.0
    k = inside.sum()
    p_hat = k / N
    pi_hat = 4.0 * p_hat
    z = 1.959963984540054  # 95%
    n = N
    denom = 1.0 + z**2 / n
    center = (p_hat + (z*z)/(2*n)) / denom
    halfwidth = z * np.sqrt((p_hat*(1-p_hat) + (z*z)/(4*n)) / n) / denom
    p_low, p_high = center - halfwidth, center + halfwidth
    pi_low, pi_high = 4.0 * p_low, 4.0 * p_high
    return {"N": N, "k_inside": int(k), "p_hat": p_hat, "pi_hat": pi_hat,
            "pi_CI95": (pi_low, pi_high), "x": x, "y": y, "inside": inside}

res = rejection_sampling_pi(N=30000)

print(f"N={res['N']}, inside={res['k_inside']} -> p_hat={res['p_hat']:.6f}")
print(f"π̂ = {res['pi_hat']:.6f}")
print(f"95% Wilson CI for π: [{res['pi_CI95'][0]:.6f}, {res['pi_CI95'][1]:.6f}]")

m = 6000
xv = res["x"][:m]
yv = res["y"][:m]
insv = res["inside"][:m]

theta = np.linspace(0, 2*np.pi, 600)
cx = np.cos(theta)
cy = np.sin(theta)

plt.figure(figsize=(6,6))
plt.scatter(xv[~insv], yv[~insv], s=8, alpha=0.6, label='rejected (outside circle)')
plt.scatter(xv[insv], yv[insv], s=8, alpha=0.6, label='accepted (inside circle)')
plt.plot(cx, cy, linewidth=2, label='unit circle boundary')
plt.text(-1, -1.2, '正方形内的提议点（蓝=拒绝，橙=接受）和单位圆边界', fontsize=12, color='red')
plt.title('Rejection Sampling for π: points in square, accept if x²+y² ≤ 1')
plt.gca().set_aspect('equal', 'box')
plt.xlim(-1,1); plt.ylim(-1,1)
plt.legend()
plt.show()

inside_prefix = res["inside"].astype(float)
running_p = np.cumsum(inside_prefix) / np.arange(1, res["N"]+1)
running_pi = 4.0 * running_p

plt.figure(figsize=(9,4))
plt.plot(running_pi)
plt.axhline(np.pi, linestyle='--', linewidth=1.5, label='π (true)')
plt.text(-1, 1.5, 'π 的运行估计随样本量收敛到 π', fontsize=12, color='red')
plt.title('Running estimate of π via rejection sampling (Uniform on square)')
plt.xlabel('N (number of proposals)')
plt.ylabel('estimate of π')
plt.legend()
plt.show()

{"pi_hat": float(res["pi_hat"]), "pi_CI95": [float(res["pi_CI95"][0]), float(res["pi_CI95"][1])], "N": int(res["N"]), "accepted_ratio": float(res["k_inside"]/res["N"])}

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

np.random.seed(2025)

def importance_estimate(samples, target_pdf, proposal_pdf, h):
    qx = proposal_pdf(samples)
    px = target_pdf(samples)
    w = px / qx
    wh = w * h(samples)
    is_hat = wh.mean() # 计算 IS 估计
    sn_hat = wh.sum() / w.sum() #  自归一化估计
    return is_hat, sn_hat, w

N1 = 5000
sigma_prop = 2.0
x_prop = np.random.normal(0.0, sigma_prop, size=N1)
target_pdf = lambda x: norm.pdf(x, 0.0, 1.0)
proposal_pdf = lambda x: norm.pdf(x, 0.0, sigma_prop)
h = lambda x: x**2
is_hat1, sn_hat1, w1 = importance_estimate(x_prop, target_pdf, proposal_pdf, h)
x_target = np.random.normal(0.0, 1.0, size=N1)
ref_est = (x_target**2).mean()
w_norm = w1 / w1.sum()
ESS1 = 1.0 / np.sum(w_norm**2)

print("=== 示例：估计 E_{N(0,1)}[X^2] ===")
print(f"目标真实值（理论） = 1.0")
print(f"\tIS (未归一) 估计 (1/n Σ w h) = {is_hat1:.6f}")
print(f"\tIS（自归一） 估计 = {sn_hat1:.6f}")
print(f"\t直接在 target 上采样的估计 (参照) = {ref_est:.6f}")
print("-> 检查两种 IS 估计（未归一与自归一）是否都接近真实值 1.0。")
print(f"样本数 N = {N1}, ESS (approx) = {ESS1:.1f} （接近 N 则说明 q 与 π 匹配较好，权重并不太集中）")
print("")

plt.figure(figsize=(8,6.5))
plt.hist(w1, bins=80)
plt.text(0.5, -120, '重要性权重 w = π(x)/q(x) 的直方图，\n注意权重分布的集中程度：\n- 是否有少数样本主导估计？\n- 权重分布是否集中在有限范围？\n- 是否出现极端爆炸？', fontsize=12, color='red')
plt.title("Importance weights histogram")
plt.xlabel("weight w = π(x)/q(x)")
plt.ylabel("count")
plt.show()

running_sn = np.cumsum(w1 * h(x_prop)) / np.cumsum(w1)
plt.figure(figsize=(9,4))
plt.plot(running_sn)
plt.axhline(1.0, linestyle='--', linewidth=1.5, label='true E[X^2]=1')
plt.text(100, 1.2, 'E[X^2] 的运行自归一化 IS 估计随样本量收敛到 1', fontsize=12, color='red')
plt.title("Running self-normalized IS estimate of E[X^2]")
plt.xlabel("n (number of proposals used)")
plt.ylabel("estimate")
plt.legend()
plt.show()

p_true = 1.0 - norm.cdf(3.0)
print("=== 示例：估计尾概率 p = P(Z>3) ===")
print(f"真实 p = {p_true:.6e}")

def single_run_compare(N, mu_shift):
    xs_naive = np.random.normal(0.0, 1.0, size=N)
    p_naive = np.mean(xs_naive > 3.0)
    xs_is = np.random.normal(loc=mu_shift, scale=1.0, size=N)
    qx = norm.pdf(xs_is, loc=mu_shift, scale=1.0)
    px = norm.pdf(xs_is, loc=0.0, scale=1.0)
    w = px / qx
    indicators = (xs_is > 3.0).astype(float)
    p_is = np.mean(w * indicators)
    ess = (w.sum()**2) / (np.sum(w**2) + 1e-300)
    return p_naive, p_is, ess, w, xs_is

N = 2000
reps = 300
mus = [3.0, 2.5]
results = {}
for mu in mus:
    p_naive_list = []
    p_is_list = []
    ess_list = []
    for r in range(reps):
        p_naive, p_is, ess, w, xs_is = single_run_compare(N, mu_shift=mu)
        p_naive_list.append(p_naive)
        p_is_list.append(p_is)
        ess_list.append(ess)
    p_naive_arr = np.array(p_naive_list)
    p_is_arr = np.array(p_is_list)
    results[mu] = {
        "p_naive_mean": p_naive_arr.mean(),
        "p_naive_std": p_naive_arr.std(ddof=1),
        "p_is_mean": p_is_arr.mean(),
        "p_is_std": p_is_arr.std(ddof=1),
        "ESS_mean": np.mean(ess_list),
        "p_naive_arr": p_naive_arr,
        "p_is_arr": p_is_arr
    }

print("比较结果（使用 N=2000, reps=300）：")
for mu in mus:
    r = results[mu]
    print(f"\nproposal N({mu},1):")
    print(f"  naive: mean={r['p_naive_mean']:.3e}, std={r['p_naive_std']:.3e} (naive Monte Carlo（直接从 N(0,1) 采样）估计的标准差)")
    print(f"  IS   : mean={r['p_is_mean']:.3e}, std={r['p_is_std']:.3e} (IS 得到的标准差)")
    print(f"  mean ESS (IS) ~ {r['ESS_mean']:.1f} (out of N={N})")
    if r['p_is_std']>0:
        print(f"  variance reduction factor (naive_var / IS_var) ≈ {(r['p_naive_std']**2)/(r['p_is_std']**2):.3f} (观察是否 IS 显著减少了方差，若有大幅度降低，则说明 IS 在估计稀有事件（尾概率）时非常有效。)")
    else:
        print("  IS variance is (near) zero")

# single-run running plots for visual comparison
xs_naive = np.random.normal(0.0, 1.0, size=N)
running_naive = np.cumsum((xs_naive > 3.0).astype(float)) / np.arange(1, N+1)
xs_is = np.random.normal(loc=3.0, scale=1.0, size=N)
qx = norm.pdf(xs_is, loc=3.0, scale=1.0); px = norm.pdf(xs_is, loc=0.0, scale=1.0)
w = px/qx; inds = (xs_is > 3.0).astype(float)
running_is = np.cumsum(w*inds) / np.arange(1, N+1)

plt.figure(figsize=(9,4))
plt.plot(running_naive, label='naive running estimate')
plt.plot(running_is, label='IS running estimate (mu=3)')
plt.axhline(p_true, linestyle='--', linewidth=1.5, label='true p')
plt.title('Running estimates: naive vs IS (single realization)')
plt.xlabel('n')
plt.ylabel('estimate of P(Z>3)')
plt.legend()
plt.show()

plt.figure(figsize=(8,4))
plt.hist(w, bins=80)
plt.text(50, 1020, '权重直方图说明权重并非完全集中到一个点，但有一定变异；\n若权重过于尖锐，会造成估计数值不稳定。', fontsize=12, color='red')
plt.title('Importance weights histogram (IS for tail prob, mu=3)')
plt.xlabel('w = π(x)/q(x)')
plt.ylabel('count')
plt.show()

# 重要性采样估计积分：比较不同 proposal 的效果，体会“好 proposal 分布 = 更低方差”。
# 积分函数
f = lambda x: 1/(1+x**2)

# 提议分布 1：均匀(0,1)
N = 10000
x = np.random.rand(N)  # g(x)=1
weights = f(x) / 1
I_est_uniform = np.mean(weights)

# 提议分布 2：Beta(0.5,0.5) → 偏向 0 和 1，贴合 f(x) 的形状
from scipy.stats import beta
x = beta.rvs(0.5,0.5,size=N)
g = beta.pdf(x,0.5,0.5)
weights = f(x)/g
I_est_beta = np.mean(weights)

print(f"真实值: {np.pi/4:.5f}")
print(f"均匀 proposal 估计: {I_est_uniform:.5f}")
print(f"Beta proposal 估计: {I_est_beta:.5f}")

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from math import exp

np.random.seed(2025)

mu1_true = np.exp(0.5)

def mc_naive_expX(N):
    X = np.random.normal(0.0, 1.0, size=N)
    g = np.exp(X)
    return g.mean(), X, g

def mc_antithetic_expX(N):
    M = (N // 2)
    Z = np.random.normal(0.0, 1.0, size=M)
    Y = (np.exp(Z) + np.exp(-Z)) / 2.0
    return Y.mean(), Z, Y

def mc_controlvar_expX(N):
    X = np.random.normal(0.0, 1.0, size=N)
    h = np.exp(X)
    Y = X**2 - 1.0
    Y_mean = Y.mean()
    #print(f"Control variates: Y (=X^2-1) mean = {Y_mean:.6f}")
    h_mean = h.mean()
    cov_hY = np.mean((h - h_mean)*(Y - Y_mean))
    var_Y = np.var(Y, ddof=0)
    beta_hat = cov_hY / (var_Y + 1e-12)
    est = np.mean(h - beta_hat * Y) # here, E[Y] = 0, so we can use h directly
    return est, X, h, Y, beta_hat

N_run = 4000
val_naive, Xn, gn = mc_naive_expX(N_run)
running_naive = np.cumsum(gn) / np.arange(1, N_run+1)

val_anti, Za, Ya = mc_antithetic_expX(N_run)
running_anti = np.cumsum(Ya) / np.arange(1, Ya.size+1)

val_cv, Xc, hc, Yc, beta_hat = mc_controlvar_expX(N_run)
cv_terms = hc - beta_hat * Yc
running_cv = np.cumsum(cv_terms) / np.arange(1, N_run+1)

plt.figure(figsize=(9,4))
plt.plot(running_naive, label='Naive running mean')
plt.axhline(mu1_true, linestyle='--', linewidth=1.5, label='True E[e^X]')
plt.title('Running estimate (Naive) for E[e^X], X~N(0,1)')
plt.xlabel('n')
plt.ylabel('estimate')
plt.legend()
plt.show()

plt.figure(figsize=(9,4))
plt.plot(running_anti, label='Antithetic running mean')
plt.axhline(mu1_true, linestyle='--', linewidth=1.5, label='True E[e^X]')
plt.title('Running estimate (Antithetic) for E[e^X], X~N(0,1)')
plt.xlabel('pairs')
plt.ylabel('estimate')
plt.legend()
plt.show()

plt.figure(figsize=(9,4))
plt.plot(running_cv, label=f'Control-Variates running mean (beta≈{beta_hat:.3f})')
plt.axhline(mu1_true, linestyle='--', linewidth=1.5, label='True E[e^X]')
plt.title('Running estimate (Control Variates) for E[e^X], X~N(0,1)')
plt.xlabel('n')
plt.ylabel('estimate')
plt.legend()
plt.show()

def replicate_stats(mu_true, N, reps=300):
    naive_vals = []
    anti_vals = []
    cv_vals = []
    betas = []
    for _ in range(reps):
        v1, *_ = mc_naive_expX(N)
        v2, *_ = mc_antithetic_expX(N)
        v3, *rest = mc_controlvar_expX(N)
        beta_hat_r = rest[-1]
        naive_vals.append(v1)
        anti_vals.append(v2)
        cv_vals.append(v3)
        betas.append(beta_hat_r)
    naive_vals = np.array(naive_vals)
    anti_vals  = np.array(anti_vals)
    cv_vals    = np.array(cv_vals)
    betas      = np.array(betas)
    sd_naive = np.std(naive_vals, ddof=1)
    sd_anti  = np.std(anti_vals,  ddof=1)
    sd_cv    = np.std(cv_vals,    ddof=1)
    vrf_anti = (sd_naive**2) / (sd_anti**2 + 1e-12)
    vrf_cv   = (sd_naive**2) / (sd_cv**2   + 1e-12)
    bias_naive = naive_vals.mean() - mu_true
    bias_anti  = anti_vals.mean()  - mu_true
    bias_cv    = cv_vals.mean()    - mu_true
    return {
        "naive_mean": naive_vals.mean(), "naive_sd": sd_naive, "naive_bias": bias_naive,
        "anti_mean": anti_vals.mean(),   "anti_sd": sd_anti,   "anti_bias": bias_anti, "VRF_anti": vrf_anti,
        "cv_mean": cv_vals.mean(),       "cv_sd": sd_cv,       "cv_bias": bias_cv,     "VRF_cv": vrf_cv,
        "beta_hat_mean": betas.mean(),   "beta_hat_sd": betas.std(ddof=1)
    }

stats_A = replicate_stats(mu1_true, N=2000, reps=300)

dfA = pd.DataFrame({
    "method": ["Naive", "Antithetic", "Control Variates"],
    "mean_estimate": [stats_A["naive_mean"], stats_A["anti_mean"], stats_A["cv_mean"]],
    "std_over_reps": [stats_A["naive_sd"], stats_A["anti_sd"], stats_A["cv_sd"]],
    "bias": [stats_A["naive_bias"], stats_A["anti_bias"], stats_A["cv_bias"]],
    "VRF_vs_Naive": [1.0, stats_A["VRF_anti"], stats_A["VRF_cv"]]
})
display(dfA)

print("=== 示例: E[e^X], X~N(0,1) ===")
print(f"True value = {mu1_true:.6f}")
print(f"Naive   : mean={stats_A['naive_mean']:.6f}, sd={stats_A['naive_sd']:.6f}")
print(f"Anti    : mean={stats_A['anti_mean']:.6f}, sd={stats_A['anti_sd']:.6f}, VRF≈{stats_A['VRF_anti']:.2f}")
print(f"Control : mean={stats_A['cv_mean']:.6f}, sd={stats_A['cv_sd']:.6f}, VRF≈{stats_A['VRF_cv']:.2f}")
print(f"(beta_hat across reps: mean≈{stats_A['beta_hat_mean']:.3f}, sd≈{stats_A['beta_hat_sd']:.3f})")
print("")

# 普通 MC 和分层采样在区间 [0,1] 的采样点分布
import numpy as np
import matplotlib.pyplot as plt

N = 10  # 分层数量

# 普通蒙特卡洛采样
mc_samples = np.random.rand(N)

# 分层采样
strata = np.linspace(0, 1, N+1)
u = np.random.rand(N)
strat_samples = strata[:-1] + u * (strata[1:] - strata[:-1])

# 绘图
plt.figure(figsize=(8, 4))

# 普通 MC
plt.subplot(1, 2, 1)
plt.scatter(mc_samples, np.zeros(N), color="blue", s=80, label="MC samples")
plt.hlines(0, 0, 1, colors="gray", linestyles="dashed")
plt.text(0.5, -0.05, '普通MC的点分布\n可能随机聚在某些地方。', fontsize=12, color='red', ha='center')
plt.title("Plain Monte Carlo (N=10)")
plt.xlabel("x in [0,1]")
plt.yticks([])
plt.legend()

# 分层采样
plt.subplot(1, 2, 2)
plt.scatter(strat_samples, np.zeros(N), color="green", s=80, label="Stratified samples")
for i in range(N+1):
    plt.axvline(strata[i], color="gray", linestyle="dashed", linewidth=0.8)  # 画出分层边界
plt.hlines(0, 0, 1, colors="black", linestyles="dashed")
plt.text(0.5, -0.05, '分层采样的点在每个小区间都有一个，\n均匀覆盖整个 [0,1]', fontsize=12, color='red', ha='center')
plt.title("Stratified Sampling (N=10)")
plt.xlabel("x in [0,1]")
plt.yticks([])
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 积分目标函数
f = lambda x: np.exp(x)
true_val = np.e - 1

def mc_integral(N=1000):
    x = np.random.rand(N)
    return f(x).mean()

def stratified_integral(N=1000):
    # 将区间[0,1]均匀分成N个子区间。例如，当N=4时，strata=[0.0,0.25,0.5,0.75,1.0]， 即分成 [0,0.25),[0.25,0.5),[0.5,0.75),[0.75,1.0)
    strata = np.linspace(0, 1, N+1)
    # 每个区间采样一个点，即生成 N 个在 [0,1] 区间的随机数。
    u = np.random.rand(N)
    # 把这些随机数“缩放到”各自的小区间
    # strata[:-1] 是所有小区间的 左端点，例如，[0,0.25,0.5,0.75]。
    # strata[1:] 是所有小区间的 右端点，例如，[0.25,0.5,0.75,1.0]。
    # u * (strata[1:] - strata[:-1]) 是每个小区间的长度乘以随机数 u，得到在每个小区间内的随机点。
    # 最终得到的 x 是在 [0,1] 区间内均匀分布的 N 个点。
    x = strata[:-1] + u * (strata[1:] - strata[:-1])
    return f(x).mean() # 计算积分估计

# 多次实验比较方差
R = 500
N = 1000
mc_vals = np.array([mc_integral(N) for _ in range(R)])
strat_vals = np.array([stratified_integral(N) for _ in range(R)])

print(f"True value = {true_val:.6f}")
print(f"MC: mean={mc_vals.mean():.6f}, std={mc_vals.std():.6f}")
print(f"Stratified: mean={strat_vals.mean():.6f}, std={strat_vals.std():.6f}")

# 可视化分布
plt.figure(figsize=(8,5))
plt.hist(mc_vals, bins=30, alpha=0.6, label="Plain MC")
plt.hist(strat_vals, bins=30, alpha=0.6, label="Stratified Sampling")
plt.axvline(true_val, color="red", linestyle="--", label="True value")
plt.title("MC vs Stratified Sampling (Integral of exp(x) on [0,1])")
plt.xlabel("Estimate value")
plt.ylabel("Frequency")
plt.legend()
plt.show()