도박사의 역설: 큰수의 법칙이 실패할 때¶

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def st_petersburg_sample_means(n_max=10_000, tries=100, n_grid=200):
    """상트페테르부르크 게임: 처음 앞면이 k번째에 나오면 2^k 을 받는다.

    상금이 2^k 이고 그 확률이 (1/2)^k 이므로 각 항의 기여가 1이고,
    항이 무한히 많아 E[X] = 1 + 1 + 1 + ... = 무한대다.
    기댓값이 없으므로 대수의 법칙이 성립하지 않는다.
    """
    # 표본 크기를 로그 눈금으로 200개 잡는다(10부터 n_max까지).
    # 로그 눈금이라야 "n이 10배가 될 때마다" 무슨 일이 생기는지 보인다.
    n_vals = np.unique(np.logspace(1, np.log10(n_max), n_grid).astype(int))
    results = []
    for n in n_vals:
        # 기하분포: 처음 성공(앞면)까지의 시행 횟수 k
        flips = np.random.geometric(0.5, size=(n, tries))
        winnings = 2.0 ** flips              # 상금
        # 같은 n에 대해 tries(=100)번 독립적으로 되풀이해 표본평균을 낸다.
        # 100개의 평균이 얼마나 흩어지는지가 이 그림의 핵심이다.
        means = winnings.mean(axis=0)
        results.append((n, means))
    return results

infinite_results = st_petersburg_sample_means()

# n이 커져도 표본평균이 한 값으로 모이지 않는다는 것을 숫자로 확인한다.
for n, means in [infinite_results[0], infinite_results[len(infinite_results)//2],
                 infinite_results[-1]]:
    print(f"n = {n:>6,}: 표본평균 100개의 중앙값 {np.median(means):8.1f}, "
          f"최댓값 {means.max():10.1f}")

n =     10: 표본평균 100개의 중앙값      5.8, 최댓값      415.4
n =    396: 표본평균 100개의 중앙값     10.2, 최댓값     1340.6
n = 10,000: 표본평균 100개의 중앙값     16.9, 최댓값      223.0

# 위에서 만든 infinite_results 를 그대로 이어받는다.
for lo, hi in [(10, 20), (95, 105), (950, 1050), (9000, 10_000)]:
    sel = [(n, m) for n, m in infinite_results if lo <= n <= hi]
    pooled = np.concatenate([m for _, m in sel])
    nbar = np.mean([n for n, _ in sel])
    q10, q50, q90 = np.percentile(pooled, [10, 50, 90])
    print(f"n ~ {nbar:>7.0f} (격자 {len(sel):>2}개):  log2 n = {np.log2(nbar):5.2f}"
          f"   중앙값 {q50:6.2f}   비 {q50 / np.log2(nbar):.2f}"
          f"   10~90% {q10:6.2f}–{q90:6.2f} ({q90/q10:.2f}배)")

n ~      15 (격자 11개):  log2 n =  3.91   중앙값   6.42   비 1.64   10~90%   3.50– 19.18 (5.48배)
n ~     100 (격자  4개):  log2 n =  6.64   중앙값   9.75   비 1.47   10~90%   5.96– 22.72 (3.81배)
n ~    1012 (격자  3개):  log2 n =  9.98   중앙값  12.48   비 1.25   10~90%   9.19– 24.56 (2.67배)
n ~    9500 (격자  4개):  log2 n = 13.21   중앙값  15.66   비 1.18   10~90%  12.41– 32.62 (2.63배)

def bounded_game_sample_means(n_max=10_000, tries=100, n_grid=200):
    """같은 게임이되 상금을 2^10 = 1024 로 잘라 낸다.

    이 한 줄의 차이로 E[X]가 유한해지고, 대수의 법칙이 되살아난다.
    상금이 유계이면 기댓값이 반드시 존재하기 때문이다.
    """
    n_vals = np.unique(np.logspace(1, np.log10(n_max), n_grid).astype(int))
    results = []
    for n in n_vals:
        flips = np.random.geometric(0.5, size=(n, tries))
        # 위 함수와 다른 곳은 이 minimum 한 군데뿐이다
        winnings = np.minimum(2.0 ** flips, 1024.0)
        means = winnings.mean(axis=0)
        results.append((n, means))
    return results

bounded_results = bounded_game_sample_means()

# 이번에는 n이 커질수록 100개의 표본평균이 한 값으로 모여든다.
# 흩어짐(표준편차)이 줄어드는지 보라.
for n, means in [bounded_results[0], bounded_results[len(bounded_results)//2],
                 bounded_results[-1]]:
    print(f"n = {n:>6,}: 표본평균 100개의 중앙값 {np.median(means):7.2f}, "
          f"표준편차 {means.std():7.2f}")

n =     10: 표본평균 100개의 중앙값    5.50, 표준편차   17.34
n =    396: 표본평균 100개의 중앙값   10.50, 표준편차    2.61
n = 10,000: 표본평균 100개의 중앙값   11.01, 표준편차    0.59

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 왼쪽 패널: 기댓값이 무한한 게임 — 발산
ax = axes[0]
for n, means in infinite_results:
    # 각 n마다 100개의 표본평균을 세로로 흩뿌린다.
    # n이 커져도 이 흩어짐이 좁아지지 않는 것이 핵심이다.
    # 세로축도 로그(loglog)로 두어야 큰 값들이 화면에 들어온다.
    ax.loglog(n * np.ones(len(means)), means, ".", color="black", ms=2, alpha=0.5)
ax.set_xlabel("n (number of rounds)")
ax.set_ylabel("Sample mean of winnings")
ax.set_title("St. Petersburg Game (E[X] = ∞)\nSample mean does NOT converge")

# 오른쪽 패널: 기댓값이 유한한 게임 — 수렴
ax = axes[1]
for n, means in bounded_results:
    # 이쪽은 세로축을 로그로 하지 않는다(semilogx). 값이 좁은 범위에 모이기 때문이다.
    # n이 커질수록 점들이 깔때기처럼 한 선으로 좁아진다. 이것이 큰수의 법칙이다.
    ax.semilogx(n * np.ones(len(means)), means, ".", color="steelblue", ms=2, alpha=0.5)

# 이론적 기댓값. 상금이 2^10 에서 잘리므로
#   E[X] = (k=1..10 각각 1씩) + 1024 * P(k > 10) = 10 + 1 = 11
true_mean = 11.0
ax.axhline(true_mean, color="red", linestyle="--", lw=2,
           label=f"E[X] = {true_mean:.1f}")
ax.set_xlabel("n (number of rounds)")
ax.set_ylabel("Sample mean of winnings")
ax.set_title("Bounded Game (E[X] < ∞)\nSample mean converges (LLN)")
ax.legend()

plt.tight_layout()
plt.show()

for tag, res in [("E = 무한", infinite_results), ("E = 11  ", bounded_results)]:
    row = []
    for lo, hi in [(10, 20), (95, 105), (950, 1050), (9000, 10_000)]:
        pooled = np.concatenate([m for n, m in res if lo <= n <= hi])
        q10, q50, q90 = np.percentile(pooled, [10, 50, 90])
        row.append(f"중앙값 {q50:6.2f} 폭 {q90/q10:4.2f}배")
    print(f"{tag}  " + " | ".join(row))

E = 무한  중앙값   6.42 폭 5.48배 | 중앙값   9.75 폭 3.81배 | 중앙값  12.48 폭 2.67배 | 중앙값  15.66 폭 2.63배
E = 11    중앙값   6.73 폭 5.55배 | 중앙값   8.80 폭 3.31배 | 중앙값  10.67 폭 1.46배 | 중앙값  11.02 폭 1.14배

import numpy as np

rng = np.random.default_rng(0)
print(f"{'n':>9}{'S_n/n 중앙값':>16}{'log2 n':>10}"
      f"{'S_n/(n log2 n) 중앙값':>24}{'사분위범위':>18}")
for n in (10 ** 2, 10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
    G = rng.geometric(0.5, (400, n))          # 첫 앞면까지 던진 횟수
    S = (2.0 ** G).sum(1)
    r = S / (n * np.log2(n))
    print(f"{n:>9}{np.median(S / n):>16.2f}{np.log2(n):>10.2f}"
          f"{np.median(r):>24.4f}"
          f"{str(np.round(np.percentile(r, [25, 75]), 3)):>18}")

        n       S_n/n 중앙값    log2 n      S_n/(n log2 n) 중앙값             사분위범위
      100            9.46      6.64                  1.4239     [1.068 2.058]
     1000           12.90      9.97                  1.2946     [1.022 1.786]
    10000           16.41     13.29                  1.2348     [1.016 1.63 ]
   100000           19.13     16.61                  1.1519     [1.034 1.363]
  1000000           22.79     19.93                  1.1436     [1.028 1.375]

import numpy as np

rng = np.random.default_rng(0)
up, down = 1.5, 0.6
print(f"산술평균 배수 {(up + down) / 2:.4f}   기하평균 배수 {np.sqrt(up * down):.4f}")
for T in (100, 1000):
    W = np.prod(np.where(rng.random((20_000, T)) < 0.5, up, down), axis=1)
    print(f"\nT = {T}")
    print(f"  평균 자산   모의 {W.mean():>11.3e}   이론 {((up + down) / 2) ** T:>11.3e}")
    print(f"  중앙값 자산 모의 {np.median(W):>11.3e}   이론 {(up * down) ** (T / 2):>11.3e}")
    print(f"  원금 미만으로 끝난 비율 {np.mean(W < 1):.4f}")
    print(f"  상위 0.1%가 전체 자산에서 차지하는 비중 "
          f"{np.sort(W)[-20:].sum() / W.sum():.4f}")

print("\n자산의 f 비율만 거는 경우의 장기 성장률")
fs = np.linspace(0.01, 1.0, 1000)
g = 0.5 * np.log(1 + 0.5 * fs) + 0.5 * np.log(1 - 0.4 * fs)
print(f"  최적 f* = {fs[g.argmax()]:.3f}   성장률 {g.max():.5f}"
      f"   (회당 배수 {np.exp(g.max()):.5f})")
for f in (0.25, 0.5, 1.0):
    gg = 0.5 * np.log(1 + 0.5 * f) + 0.5 * np.log(1 - 0.4 * f)
    print(f"  f = {f:.2f}: 성장률 {gg:>9.5f}   (회당 배수 {np.exp(gg):.5f})")

산술평균 배수 1.0500   기하평균 배수 0.9487

T = 100
  평균 자산   모의   6.622e+01   이론   1.315e+02
  중앙값 자산 모의   5.154e-03   이론   5.154e-03
  원금 미만으로 끝난 비율 0.8658
  상위 0.1%가 전체 자산에서 차지하는 비중 0.8353

T = 1000
  평균 자산   모의   3.656e-03   이론   1.546e+21
  중앙값 자산 모의   1.322e-23   이론   1.322e-23
  원금 미만으로 끝난 비율 0.9999
  상위 0.1%가 전체 자산에서 차지하는 비중 1.0000

자산의 f 비율만 거는 경우의 장기 성장률
  최적 f* = 0.250   성장률 0.00621   (회당 배수 1.00623)
  f = 0.25: 성장률   0.00621   (회당 배수 1.00623)
  f = 0.50: 성장률   0.00000   (회당 배수 1.00000)
  f = 1.00: 성장률  -0.05268   (회당 배수 0.94868)

import numpy as np

rng = np.random.default_rng(0)
n = 1_000_000
print(f"파레토 표본 n = {n:,}: 총합에서 상위 몇 건이 차지하는 비중")
print(f"{'alpha':>7}{'평균':>8}{'분산':>8}{'상위 1건':>11}{'상위 10건':>12}"
      f"{'상위 0.1%':>12}{'상위 1%':>11}")
for a in (1.05, 1.5, 2.5, 5.0):
    x = np.sort(rng.pareto(a, n) + 1)[::-1]
    tot = x.sum()
    print(f"{a:>7.2f}{'유한':>8}{('유한' if a > 2 else '무한'):>8}"
          f"{x[0] / tot:>11.4f}{x[:10].sum() / tot:>12.4f}"
          f"{x[:n // 1000].sum() / tot:>12.4f}{x[:n // 100].sum() / tot:>11.4f}")

rng = np.random.default_rng(1)
print("\nalpha = 1.5 (참 평균 3): 표본평균은 어느 쪽으로 틀리는가")
print(f"{'n':>9}{'P(표본평균 < 3)':>18}{'표본평균의 중앙값':>20}"
      f"{'표본평균의 평균':>18}")
for n in (10 ** 3, 10 ** 4, 10 ** 5, 10 ** 6):
    m = np.array([rng.pareto(1.5, n).mean() + 1 for _ in range(2000)])
    print(f"{n:>9}{np.mean(m < 3):>18.4f}{np.median(m):>20.4f}{m.mean():>18.4f}")

파레토 표본 n = 1,000,000: 총합에서 상위 몇 건이 차지하는 비중
  alpha      평균      분산      상위 1건      상위 10건     상위 0.1%      상위 1%
   1.05      유한      무한     0.0167      0.1044      0.3957     0.5777
   1.50      유한      무한     0.0045      0.0168      0.0958     0.2113
   2.50      유한      유한     0.0004      0.0012      0.0158     0.0630
   5.00      유한      유한     0.0000      0.0001      0.0040     0.0251

alpha = 1.5 (참 평균 3): 표본평균은 어느 쪽으로 틀리는가
        n       P(표본평균 < 3)           표본평균의 중앙값          표본평균의 평균
     1000            0.6750              2.8655            3.1038
    10000            0.6705              2.9402            3.0097
   100000            0.6580              2.9723            2.9994
  1000000            0.6605              2.9868            2.9978