동전 던지기 모의실험¶

import numpy as np

np.random.seed(42)

TOTAL_TOSSES = 30
OBSERVED_HEADS = 24
PROB_HEAD_FAIR = 0.5
NUM_SIMULATIONS = 100_000

def single_experiment(n_tosses=TOTAL_TOSSES, p=PROB_HEAD_FAIR):
    """공정한 동전을 n_tosses번 던진 한 판을 흉내 내고 앞면 횟수를 돌려준다."""
    # 0/1을 30개 뽑아 더할 필요가 없다. 앞면 횟수의 분포가 곧 Bin(30, 0.5)다.
    return np.random.binomial(n_tosses, p)

# H0가 참일 때 한 판을 돌리면 무엇이 나오는지 몇 번 본다.
print([single_experiment() for _ in range(10)])

[14, 20, 17, 16, 12, 12, 11, 18, 16, 17]

def simulate_coin_tosses(n_simulations=NUM_SIMULATIONS,
                         n_tosses=TOTAL_TOSSES,
                         p=PROB_HEAD_FAIR):
    """실험을 n_simulations번 반복하고 앞면 횟수 배열을 돌려준다."""
    return np.random.binomial(n_tosses, p, size=n_simulations)

# 여기서 세는 것은 "H0가 참일 때 관측값만큼 극단적인 일이 얼마나 자주 일어나는가"다.
# 그것이 p-값의 정의다. 이항분포 공식을 몰라도 이 논리는 그대로 성립한다.

head_counts = simulate_coin_tosses()
extreme = np.sum(head_counts >= OBSERVED_HEADS)
pct = extreme / NUM_SIMULATIONS * 100

print(f"Times with >= {OBSERVED_HEADS} heads: {extreme:,}")
print(f"Percentage: {pct:.4f}%")

# (1) 의 예측과 견준다.
p_star = 0.000715453
mean_y = NUM_SIMULATIONS * p_star
sd_y = np.sqrt(NUM_SIMULATIONS * p_star * (1 - p_star))
print(f"예측 개수 {mean_y:.3f} ± {sd_y:.4f}"
      f"   관측과의 거리 {(extreme - mean_y) / sd_y:+.3f} SD")

# 모의실험 p-값의 표준오차는 관측된 비율로도 잴 수 있다.
p_hat = extreme / NUM_SIMULATIONS
se_hat = np.sqrt(p_hat * (1 - p_hat) / NUM_SIMULATIONS)
print(f"p-hat = {p_hat:.6f},  SE = {se_hat:.3e},"
      f"  95% 구간 ({p_hat - 1.96 * se_hat:.6f}, {p_hat + 1.96 * se_hat:.6f})")

Times with >= 24 heads: 71
Percentage: 0.0710%
예측 개수 71.545 ± 8.4554   관측과의 거리 -0.064 SD
p-hat = 0.000710,  SE = 8.423e-05,  95% 구간 (0.000545, 0.000875)

from scipy.stats import binom

# P(X >= 24) = 1 - P(X <= 23) 이다. cdf에 24가 아니라 **23**을 넣어야 한다.
# 이산분포에서 부등호를 하나 어긋나게 쓰는 것이 가장 흔한 실수다.
p_exact = 1 - binom.cdf(OBSERVED_HEADS - 1, TOTAL_TOSSES, PROB_HEAD_FAIR)
print(f"Exact binomial P(X >= {OBSERVED_HEADS}): {p_exact:.6f}")

# 한 칸 어긋나게 쓰면 어떻게 되는가.
p_wrong = 1 - binom.cdf(OBSERVED_HEADS, TOTAL_TOSSES, PROB_HEAD_FAIR)
pmf24 = binom.pmf(OBSERVED_HEADS, TOTAL_TOSSES, PROB_HEAD_FAIR)
print(f"틀린 1 - F(24)          : {p_wrong:.6f}   (= P(X >= 25))")
print(f"빠뜨린 P(X = 24)        : {pmf24:.6f}"
      f"   꼬리에서 차지하는 몫 {pmf24 / p_exact:.1%}")
print(f"비 (옳은 값)/(틀린 값)  : {p_exact / p_wrong:.3f}")
print(f"모의실험 값             : {pct / 100:.6f}"
      f"   차이 {abs(pct / 100 - p_exact) / se_hat:.3f} SE")

Exact binomial P(X >= 24): 0.000715
틀린 1 - F(24)          : 0.000162   (= P(X >= 25))
빠뜨린 P(X = 24)        : 0.000553   꼬리에서 차지하는 몫 77.3%
비 (옳은 값)/(틀린 값)  : 4.404
모의실험 값             : 0.000710   차이 0.065 SE

import matplotlib.pyplot as plt

# 그림에서 읽을 수치를 미리 찍어 둔다.
counts = np.bincount(head_counts, minlength=TOTAL_TOSSES + 1)
print(f"모의 평균 {head_counts.mean():.4f}  표준편차 {head_counts.std(ddof=1):.4f}"
      f"   (이론 15, {np.sqrt(7.5):.4f})")
print(f"최솟값 {head_counts.min()}  최댓값 {head_counts.max()}")
print(f"가장 높은 막대: x = {counts.argmax()}, 높이 {counts.max():,}")
for x in (22, 23, 24, 25, 26):
    print(f"  x = {x}: 높이 {counts[x]:>5,}"
          f"   가장 높은 막대의 {counts[x] / counts.max():.3%}")

# 앞면 수가 정수이므로 계급 경계를 반 칸씩 밀어 막대 하나가 값 하나를 담게 한다.
fig, ax = plt.subplots(figsize=(8, 5))
bins = np.arange(0, TOTAL_TOSSES + 2) - 0.5
ax.hist(head_counts, bins=bins, edgecolor="white", alpha=0.7,
        label="Simulated head counts")
# 관측값 자리에 세로선을 긋는다. 그 오른쪽 막대들의 넓이 비율이 곧 p-값이다.
ax.axvline(OBSERVED_HEADS, color="red", linestyle="--", linewidth=2,
           label=f"Observed = {OBSERVED_HEADS}")
ax.set_xlabel("Number of heads")
ax.set_ylabel("Frequency")
ax.set_title(f"Coin Toss Simulation ({NUM_SIMULATIONS:,} runs)")
ax.legend()
plt.tight_layout()
plt.show()

모의 평균 14.9923  표준편차 2.7302   (이론 15, 2.7386)
최솟값 4  최댓값 26
가장 높은 막대: x = 15, 높이 14,442
  x = 22: 높이   535   가장 높은 막대의 3.704%
  x = 23: 높이   180   가장 높은 막대의 1.246%
  x = 24: 높이    56   가장 높은 막대의 0.388%
  x = 25: 높이    12   가장 높은 막대의 0.083%
  x = 26: 높이     3   가장 높은 막대의 0.021%

head_counts = np.random.binomial(30, 0.5, size=100_000)
extreme_two_sided = np.sum((head_counts >= 24) | (head_counts <= 6))
p_two_sided = extreme_two_sided / 100_000
print(f"Two-sided simulated p-value: {p_two_sided:.4f}")

Two-sided simulated p-value: 0.0014

import numpy as np

print(f"{'B':>7s} {'단순 b/B':>10s} {'보정 (b+1)/(B+1)':>18s}")
for B in [20, 100, 1000, 10000]:
    a = 0.05
    print(f"{B:7d} {(np.floor(a * B) + 1) / (B + 1):10.4f} "
          f"{np.floor(a * (B + 1)) / (B + 1):18.4f}")
print()
for B in [19, 99, 999, 9999]:
    a = 0.05
    print(f"{B:7d} {(np.floor(a * B) + 1) / (B + 1):10.4f} "
          f"{np.floor(a * (B + 1)) / (B + 1):18.4f}")

      B     단순 b/B     보정 (b+1)/(B+1)
     20     0.0952             0.0476
    100     0.0594             0.0495
   1000     0.0509             0.0500
  10000     0.0501             0.0500

     19     0.0500             0.0500
     99     0.0500             0.0500
    999     0.0500             0.0500
   9999     0.0500             0.0500

import numpy as np
from scipy import stats

print(f"{'n':>5s} {'정확검정':>10s} {'중간-p':>10s}")
for n in [10, 20, 30, 50]:
    k = np.arange(n + 1)
    pmf = stats.binom.pmf(k, n, 0.5)
    exact = np.array([stats.binomtest(int(j), n, 0.5).pvalue for j in k])
    midp = np.empty(n + 1)
    for j in k:
        one = (stats.binom.cdf(j - 1, n, 0.5) + 0.5 * pmf[j] if j <= n / 2
               else stats.binom.sf(j, n, 0.5) + 0.5 * pmf[j])
        midp[j] = min(2 * one, 1)
    print(f"{n:5d} {pmf[exact <= 0.05].sum():10.4f} "
          f"{pmf[midp <= 0.05].sum():10.4f}")

    n       정확검정       중간-p
   10     0.0215     0.0215
   20     0.0414     0.0414
   30     0.0428     0.0428
   50     0.0328     0.0649

import numpy as np
from scipy import stats

def exact_power(n, p1, p0=0.5, alpha=0.05):
    k = np.arange(n + 1)
    pv = np.array([stats.binomtest(int(j), n, p0).pvalue for j in k])
    reject = pv <= alpha
    return stats.binom.pmf(k, n, p1)[reject].sum(), \
           stats.binom.pmf(k, n, p0)[reject].sum()

print(f"{'n':>5s} {'실제 수준':>10s} {'p=0.6':>8s} {'p=0.7':>8s} {'p=0.8':>8s}")
for n in [20, 30, 50, 100, 200]:
    lvl = exact_power(n, 0.5)[1]
    row = [exact_power(n, p)[0] for p in [0.6, 0.7, 0.8]]
    print(f"{n:5d} {lvl:10.4f} " + " ".join(f"{v:8.4f}" for v in row))

print()
for p1 in [0.6, 0.7, 0.8]:
    n = next(m for m in range(5, 2000) if exact_power(m, p1)[0] >= 0.80)
    print(f"p1 = {p1}: 검정력 80%를 위해 n = {n}")

    n      실제 수준    p=0.6    p=0.7    p=0.8
   20     0.0414   0.1272   0.4164   0.8042
   30     0.0428   0.1771   0.5888   0.9389
   50     0.0328   0.2371   0.7822   0.9937
  100     0.0352   0.4621   0.9790   1.0000
  200     0.0400   0.7868   0.9999   1.0000

p1 = 0.6: 검정력 80%를 위해 n = 199
p1 = 0.7: 검정력 80%를 위해 n = 49
p1 = 0.8: 검정력 80%를 위해 n = 20

import numpy as np
from scipy import stats
from scipy.special import betaln

def bf10(n, k, a=1.0, b=1.0):
    return np.exp(betaln(k + a, n - k + b) - betaln(a, b) - n * np.log(0.5))

print(f"{'n':>6s} {'k':>6s} {'p-값':>10s} {'BF10(균등)':>12s} "
      f"{'BF10(제프리스)':>14s}")
for n, k in [(30, 24), (30, 21), (100, 61), (1000, 531)]:
    pv = stats.binomtest(k, n, 0.5).pvalue
    print(f"{n:6d} {k:6d} {pv:10.5f} {bf10(n, k):12.3f} "
          f"{bf10(n, k, 0.5, 0.5):14.3f}")

     n      k        p-값     BF10(균등)     BF10(제프리스)
    30     24    0.00143       58.333         46.736
    30     21    0.04277        2.421          1.704
   100     61    0.03520        1.392          0.913
  1000    531    0.05368        0.270          0.173

import numpy as np

rng = np.random.default_rng(9)
Bmax, alpha = 100_000, 0.05
M = 2_000

def sequential(u, h=20):
    """b가 h에 도달하면 조기 중단. 반환: (판정, 사용한 모의실험 수)"""
    b = 0
    for i in range(1, Bmax + 1):
        b += rng.random() < u
        if b >= h:                       # 충분히 많이 쌓임 → p 는 크다
            return "비기각", i
        if i >= Bmax:
            break
    return ("기각" if (b + 1) / (i + 1) <= alpha else "비기각"), i

used = []
for u in rng.random(M):                  # H0 아래의 참 p-값
    _, i = sequential(u)
    used.append(i)
used = np.array(used)
print(f"평균 사용 횟수 {used.mean():10.1f}  (고정 B = {Bmax:,d})")
print(f"중앙값 {np.median(used):10.1f}   90 백분위 {np.percentile(used, 90):10.1f}")
print(f"절약 비율 {1 - used.mean() / Bmax:.3%}")

평균 사용 횟수      130.3  (고정 B = 100,000)
중앙값       39.0   90 백분위      176.1
절약 비율 99.870%