9.2 期望值與變異數的計算#

import numpy as np, random, math

# --- Welford 單趟演算法 vs 天真做法 ---
def welford(data):
    n = mean = m2 = 0
    for x in data:
        n += 1
        d = x - mean
        mean += d / n
        m2 += d * (x - mean)
    return mean, m2 / (n - 1)

big = 1e9 + np.concatenate([np.zeros(10),  # 巨大均值 + 微小波動
                            np.random.default_rng(0).normal(0, 1, 10)])
mu1, var1 = welford(big)
mu2 = big.mean(); var2 = big.var(ddof=1)
print(f"Welford 變異數 = {var1:.6f}")
print(f"兩步法變異數   = {var2:.6f}（真值 ≈ 1.05，兩步法有抵銷誤差）")

# --- 反變換法：模擬幾何分布（丟銅板到第一次正面） ---
def geom(u, p=0.5):
    return math.ceil(math.log(1 - u) / math.log(1 - p))

rng = np.random.default_rng(1)
u = rng.uniform(0, 1, 200_000)
xs = [geom(v) for v in u]
print(f"幾何分布模擬平均 {np.mean(xs):.3f}（理論 1/p = 2）")
print(f"模擬變異數 {np.var(xs):.3f}（理論 (1-p)/p² = 2）")

# --- Box-Muller：產生常態亂數並驗證 CLT ---
def box_muller(u1, u2):
    r = np.sqrt(-2 * np.log(u1))
    return r * np.cos(2 * np.pi * u2)

u1, u2 = rng.uniform(0, 1, 100_000), rng.uniform(0, 1, 100_000)
z = box_muller(u1, u2)
print(f"Box-Muller 樣本: 平均 {z.mean():.4f}, 標準差 {z.std():.4f}, "
      f"偏度 {float(((z - z.mean())**3).mean() / z.std()**3):.4f}")

# CLT 驗證：擲骰子 30 次取平均，重複 20 萬次
rolls = rng.integers(1, 7, size=(200_000, 30))
sample_means = rolls.mean(axis=1)
mu_theory, sd_theory = 3.5, math.sqrt(35/12) / math.sqrt(30)
print(f"樣本平均 ≈ {sample_means.mean():.4f}（理論 3.5），"
      f"標準誤 ≈ {sample_means.std():.4f}（理論 {sd_theory:.4f}）")

# --- 收集郵票問題 E = n·Hn ---
def coupon_collect(n, trials=20_000):
    total = 0
    for _ in range(trials):
        seen, count = set(), 0
        while len(seen) < n:
            seen.add(rng.integers(0, n)); count += 1
        total += count
    return total / trials

for n in [5, 10]:
    print(f"集滿 {n} 種：模擬 {coupon_collect(n):.1f} 次，"
          f"理論 n·Hn = {n * sum(1/k for k in range(1, n+1)):.1f}")