큰수의 법칙¶

import numpy as np
import matplotlib.pyplot as plt

N_LIST = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
M = 10_000                    # 되풀이 횟수
EPS = 0.2                     # ε (σ 단위)

rng = np.random.default_rng(2026)
mu, sigma = 1.0, 1.0          # Exponential(1) 의 평균과 표준편차
eps = EPS * sigma

# 가장 큰 n 으로 한 번만 뽑고 앞쪽 n 개 열을 각 칸에 쓴다.
# 표본이 중첩되므로 칸끼리 다른 것은 오직 n 뿐이다.
X = rng.exponential(1.0, (M, max(N_LIST)))

lo, hi = np.percentile(X[:, :min(N_LIST)].mean(axis=1), [0.5, 99.5])
pad = 0.1 * (hi - lo)
bins = np.linspace(lo - pad, hi + pad, 60)

probs = []
fig, axes = plt.subplots(1, 10, figsize=(22, 4), constrained_layout=True)
for ax, n in zip(axes, N_LIST):
    xbar = X[:, :n].mean(axis=1)          # 크기 n 인 표본 M 개의 평균
    ax.hist(xbar, bins=bins, density=True, color="tab:blue", alpha=0.6,
            edgecolor="white", linewidth=0.3)
    ax.axvline(mu, color="red", lw=2, label=r"$\mu$")
    ax.axvspan(mu - eps, mu + eps, color="orange", alpha=0.18,
               label=r"$\mu\pm\varepsilon$")

    # 정리 2가 "0으로 간다"고 말하는 바로 그 확률을 세어 본다.
    p_out = np.mean(np.abs(xbar - mu) > eps)
    probs.append(p_out)
    ax.set_title(f"n = {n}\n"
                 rf"$\hat P(|\bar X_n-\mu|>\varepsilon)$ = {p_out:.3f}",
                 fontsize=11)
    ax.set_xlim(bins[0], bins[-1])
    ax.set_xlabel(r"$\bar X_n$")
axes[0].set_ylabel("density of $\\bar X_n$")
axes[0].legend(fontsize=9, loc="upper right")
plt.show()

# 정리 2가 0으로 간다고 말하는 그 확률만 따로 모아 본다.
for i in range(5):
    print(f"n = {N_LIST[i]:>2} : {probs[i]:.3f}      "
          f"n = {N_LIST[i + 5]:>2} : {probs[i + 5]:.3f}")

n =  5 : 0.661      n = 30 : 0.269
n = 10 : 0.522      n = 35 : 0.229
n = 15 : 0.434      n = 40 : 0.198
n = 20 : 0.370      n = 45 : 0.169
n = 25 : 0.312      n = 50 : 0.147

from scipy import stats

print(" n   정확(감마)  정규근사  체비쇼프   모의     차이/se")
for n, p in zip(N_LIST, probs):
    exact = (stats.gamma.cdf((mu - eps) * n, n)
             + stats.gamma.sf((mu + eps) * n, n))
    se = np.sqrt(exact * (1 - exact) / M)      # 모의값의 몬테카를로 오차
    print(f"{n:3d}    {exact:.4f}     {2 * stats.norm.sf(eps * np.sqrt(n)):.4f}"
          f"   {sigma**2 / (n * eps**2):7.3f}   {p:.3f}   {(p - exact) / se:+5.2f}")

 n   정확(감마)  정규근사  체비쇼프   모의     차이/se
  5    0.6562     0.6547     5.000   0.661   +1.03
 10    0.5258     0.5271     2.500   0.522   -0.77
 15    0.4361     0.4386     1.667   0.434   -0.45
 20    0.3680     0.3711     1.250   0.370   +0.41
 25    0.3140     0.3173     1.000   0.312   -0.39
 30    0.2700     0.2733     0.833   0.269   -0.18
 35    0.2335     0.2367     0.714   0.229   -1.19
 40    0.2029     0.2059     0.625   0.198   -1.11
 45    0.1769     0.1797     0.556   0.169   -2.10
 50    0.1547     0.1573     0.500   0.147   -2.17

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# 주사위를 1만 번 굴리고 "여태까지의 평균"을 매 시점 기록한다.
n_rolls = 10_000
rolls = np.random.randint(1, 7, size=n_rolls)

# cumsum(누적합)을 1, 2, 3, ... 로 나누면 각 시점까지의 평균이 된다.
# running_avg[k] = 처음 k+1번의 평균
# 반복문 없이 한 줄로 1만 개의 평균을 얻는 요령이다.
running_avg = np.cumsum(rolls) / np.arange(1, n_rolls + 1)

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(running_avg, alpha=0.8, label='Running Average')
# 참 기댓값 (1+2+3+4+5+6)/6 = 3.5. 큰수의 법칙은 곡선이 이 선에 붙는다고 말한다.
ax.axhline(y=3.5, color='r', linestyle='--', label='E[X] = 3.5')
ax.set_xlabel('Number of Rolls')
ax.set_ylabel('Sample Mean')
ax.set_title('Law of Large Numbers: Fair Die')
ax.legend()
ax.spines[['top', 'right']].set_visible(False)
plt.show()

# 그림에서 눈으로 재던 것을 수로 적는다.
sigma = np.sqrt(35 / 12)
for k in (10, 100, 1_000, 10_000):
    dev = running_avg[k - 1] - 3.5
    print(f"n = {k:>6,}: 누적평균 {running_avg[k-1]:.4f}  편차 {dev:+.4f}"
          f"  이론 sd {sigma / np.sqrt(k):.4f}  편차/sd {dev / (sigma / np.sqrt(k)):+.2f}")

# 누적합 쪽은 어떻게 되는가. S_n - 3.5n 은 줄지 않는다.
dev_sum = np.cumsum(rolls) - 3.5 * np.arange(1, n_rolls + 1)
k = int(np.abs(dev_sum).argmax())
print(f"|S_n - 3.5n| 의 최댓값 = {abs(dev_sum[k]):.0f}  (n = {k+1:,}),"
      f"  sigma*sqrt(n) = {sigma * np.sqrt(k + 1):.0f}")
out = np.where(np.abs(running_avg - 3.5) > 0.05)[0]
print(f"3.5 +- 0.05 띠를 마지막으로 벗어난 때 n = {out[-1] + 1:,}")

n =     10: 누적평균 3.8000  편차 +0.3000  이론 sd 0.5401  편차/sd +0.56
n =    100: 누적평균 3.6900  편차 +0.1900  이론 sd 0.1708  편차/sd +1.11
n =  1,000: 누적평균 3.4570  편차 -0.0430  이론 sd 0.0540  편차/sd -0.80
n = 10,000: 누적평균 3.4999  편차 -0.0001  이론 sd 0.0171  편차/sd -0.01
|S_n - 3.5n| 의 최댓값 = 134  (n = 6,730),  sigma*sqrt(n) = 140
3.5 +- 0.05 띠를 마지막으로 벗어난 때 n = 663

import numpy as np

rng = np.random.default_rng(0)
print("표본평균의 흩어짐 (각 n마다 200회 반복한 사분위범위)")
print(f"{'n':>9}{'파레토 a=1.5':>16}{'지수 (평균 3)':>16}{'배수':>8}")
for n in (10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
    pa = np.array([rng.pareto(1.5, n).mean() + 1 for _ in range(200)])
    ex = np.array([rng.exponential(3.0, n).mean() for _ in range(200)])
    iqr = lambda v: np.subtract(*np.percentile(v, [75, 25]))
    print(f"{n:>9}{iqr(pa):>16.5f}{iqr(ex):>16.5f}{iqr(pa)/iqr(ex):>8.1f}")

표본평균의 흩어짐 (각 n마다 200회 반복한 사분위범위)
        n       파레토 a=1.5       지수 (평균 3)      배수
     1000         0.38012         0.12907     2.9
    10000         0.17643         0.03738     4.7
   100000         0.09132         0.01276     7.2
  1000000         0.03894         0.00396     9.8

import numpy as np
from scipy import integrate

rng = np.random.default_rng(0)
c = integrate.quad(lambda t: np.exp(-t ** 2), 0, 1)[0]   # I_d = c^d
print(f"I_d = ∫_[0,1]^d exp(-|x|^2) dx = c^d,   c = {c:.6f}\n")

print(f"{'d':>4}{'참값':>12}{'몬테카를로 n=10^5':>20}{'상대오차':>11}{'격자 점 수':>14}")
for d in (1, 3, 10, 30, 100):
    U = rng.random((100_000, d))
    est = np.exp(-(U ** 2).sum(1)).mean()
    print(f"{d:>4}{c ** d:>12.3e}{est:>20.3e}{abs(est - c ** d) / c ** d:>11.2%}"
          f"{'10^' + str(d):>14}")

print("\n피적분함수의 변동계수 (오차 상수를 결정한다)")
for d in (1, 10, 100):
    v = np.exp(-(rng.random((200_000, d)) ** 2).sum(1))
    print(f"  d={d:>3}: sd/mean = {v.std() / v.mean():>8.4f}")

I_d = ∫_[0,1]^d exp(-|x|^2) dx = c^d,   c = 0.746824

   d          참값        몬테카를로 n=10^5       상대오차        격자 점 수
   1   7.468e-01           7.471e-01      0.04%          10^1
   3   4.165e-01           4.165e-01      0.00%          10^3
  10   5.397e-02           5.379e-02      0.34%         10^10
  30   1.572e-04           1.579e-04      0.41%         10^30
 100   2.098e-13           1.925e-13      8.23%        10^100

피적분함수의 변동계수 (오차 상수를 결정한다)
  d=  1: sd/mean =   0.2690
  d= 10: sd/mean =   1.0087
  d=100: sd/mean =  21.3544

import numpy as np

def ar1(n, phi, rng):
    e = rng.normal(0, 1, n)
    x = np.empty(n)
    x[0] = e[0] / np.sqrt(1 - phi ** 2)          # 정상분포에서 출발
    for i in range(1, n):
        x[i] = phi * x[i - 1] + e[i]
    return x

rng = np.random.default_rng(0)
n = 10_000
print(f"{'phi':>6}{'표본평균 sd':>15}{'독립 가정 sd':>15}{'실제 n_eff':>13}{'이론 n_eff':>13}")
for phi in (0.0, 0.5, 0.9, 0.99):
    m = np.array([ar1(n, phi, rng).mean() for _ in range(300)])
    sd_indep = np.sqrt(1 / (1 - phi ** 2) / n)
    print(f"{phi:>6.2f}{m.std():>15.5f}{sd_indep:>15.5f}"
          f"{n * (sd_indep / m.std()) ** 2:>13.1f}"
          f"{n * (1 - phi) / (1 + phi):>13.1f}")

   phi        표본평균 sd       독립 가정 sd     실제 n_eff     이론 n_eff
  0.00        0.01011        0.01000       9789.6      10000.0
  0.50        0.01995        0.01155       3350.6       3333.3
  0.90        0.09962        0.02294        530.4        526.3
  0.99        1.03067        0.07089         47.3         50.3

import numpy as np
from scipy.special import logsumexp

grid = np.linspace(-2, 60, 200_001)
Q = -grid ** 2 / (1 + grid ** 2)                       # 최대는 t=0
print("Q(t) = -t^2/(1+t^2),   Q_n(t) = Q(t) + 3 exp(-4 (t-n)^2)")
for n in (5, 10, 50):
    Qn = Q + 3 * np.exp(-4 * (grid - n) ** 2)
    print(f"  n={n:>3}:  argmax Q_n = {grid[Qn.argmax()]:>6.2f}")
print(f"  참 argmax Q = {grid[Q.argmax()]:.2f}"
      f"   (그런데 고정된 t 마다 Q_n(t) -> Q(t) 이다)")

rng = np.random.default_rng(0)
x = rng.normal(0, 1, 20)

def loglik(mu1, s1, mu2=0.0, s2=1.0, w=0.5):
    lp = np.log(2 * np.pi) / 2
    a = np.log(w) - np.log(s1) - lp - 0.5 * ((x - mu1) / s1) ** 2
    b = np.log(1 - w) - np.log(s2) - lp - 0.5 * ((x - mu2) / s2) ** 2
    return logsumexp(np.vstack([a, b]), axis=0).sum()

print("\n혼합 0.5 N(mu1, s1^2) + 0.5 N(0, 1),  자료 20개 ~ N(0,1)")
print(f"  정직한 해 (mu1=0, s1=1):  로그가능도 {loglik(0.0, 1.0):>9.2f}")
print(f"{'s1':>10}{'로그가능도':>16}   <- mu1 을 자료점 하나에 고정")
for s1 in (1e-1, 1e-3, 1e-6, 1e-10, 1e-20):
    print(f"{s1:>10.0e}{loglik(x[0], s1):>16.2f}")

Q(t) = -t^2/(1+t^2),   Q_n(t) = Q(t) + 3 exp(-4 (t-n)^2)
  n=  5:  argmax Q_n =   5.00
  n= 10:  argmax Q_n =  10.00
  n= 50:  argmax Q_n =  50.00
  참 argmax Q = 0.00   (그런데 고정된 t 마다 Q_n(t) -> Q(t) 이다)

혼합 0.5 N(mu1, s1^2) + 0.5 N(0, 1),  자료 20개 ~ N(0,1)
  정직한 해 (mu1=0, s1=1):  로그가능도    -25.95
        s1           로그가능도   <- mu1 을 자료점 하나에 고정
     1e-01          -31.94
     1e-03          -32.90
     1e-06          -25.99
     1e-10          -16.78
     1e-20            6.24