베리–에센 정리¶

import numpy as np
from scipy import stats

def berry_esseen_bound(sigma, rho, n, C=0.4748):
    """베리-에센 상계를 계산한다.

        |F_n(x) - Phi(x)|  <=  C * rho / (sigma^3 * sqrt(n))

    이 값은 **모든 x에 대해 한꺼번에** 성립하는 오차 상한이다.
    중심극한정리가 "n이 크면 정규분포에 가까워진다"고만 말하는 데 비해,
    베리-에센은 "얼마나 가까운가"를 수치로 준다.
    """
    return C * rho / (sigma**3 * np.sqrt(n))

# 베르누이(0.5): 대칭이라 정규근사가 빠르다
sigma_b = np.sqrt(0.25)     # 표준편차 0.5
rho_b = 0.125               # E|X - mu|^3 = 0.5^3 = 0.125
print("=== Bernoulli(0.5) ===")
for n in [10, 30, 100, 1000]:
    bound = berry_esseen_bound(sigma_b, rho_b, n)
    print(f"n = {n:5d}: Berry–Esseen bound = {bound:.4f}")

print()

# 지수분포(1): 오른쪽으로 크게 치우쳐 정규근사가 훨씬 느리다.
# rho/sigma^3 이 2.415 대 1 로 2.4배 크므로 상계도 그만큼 커진다.
# rho = E|X - 1|^3 = 12/e - 2 (지수분포의 3차 절대적률)
sigma_e = 1.0
rho_e = 12 / np.e - 2
print("=== Exponential(1) ===")
for n in [10, 30, 100, 1000]:
    bound = berry_esseen_bound(sigma_e, rho_e, n)
    print(f"n = {n:5d}: Berry–Esseen bound = {bound:.4f}")

=== Bernoulli(0.5) ===
n =    10: Berry–Esseen bound = 0.1501
n =    30: Berry–Esseen bound = 0.0867
n =   100: Berry–Esseen bound = 0.0475
n =  1000: Berry–Esseen bound = 0.0150

=== Exponential(1) ===
n =    10: Berry–Esseen bound = 0.3625
n =    30: Berry–Esseen bound = 0.2093
n =   100: Berry–Esseen bound = 0.1146
n =  1000: Berry–Esseen bound = 0.0363

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

plt.rcParams['font.family'] = 'Apple SD Gothic Neo'
plt.rcParams['axes.unicode_minus'] = False

def berry_esseen_visualization(dist_name, rvs_fn, mu, sigma, rho, sample_sizes):
    """실제 오차와 베리-에센 상계를 나란히 그린다.

    상계는 "이보다 나쁘지는 않다"는 보장일 뿐이므로,
    실제 오차가 상계보다 훨씬 작게 나오는 것이 정상이다.
    """
    C = 0.4748
    x_grid = np.linspace(-4, 4, 1000)     # 오차를 재는 x 격자
    n_sim = 50_000                        # 표본평균을 5만 번 만든다

    fig, axes = plt.subplots(1, len(sample_sizes), figsize=(12, 3),
                             sharey=True)
    fig.suptitle(f'베리–에센: {dist_name}', fontsize=14)

    for ax, n in zip(axes, sample_sizes):
        np.random.seed(42)
        # 크기 n인 표본을 5만 개 만들어 각각의 평균을 낸다
        samples = rvs_fn(size=(n_sim, n))
        x_bar = samples.mean(axis=1)

        # 표준화: (X-bar - mu) / (sigma/sqrt(n)).
        # 중심극한정리는 이 z가 N(0,1)로 간다고 말한다.
        z = (x_bar - mu) / (sigma / np.sqrt(n))

        # 경험적 CDF와 표준정규 CDF의 차이를 x마다 잰다
        ecdf = np.array([np.mean(z <= x) for x in x_grid])
        ncdf = stats.norm.cdf(x_grid)
        actual_error = np.abs(ecdf - ncdf)
        max_error = actual_error.max()      # 베리-에센이 상한을 주는 대상

        bound = C * rho / (sigma**3 * np.sqrt(n))

        ax.plot(x_grid, actual_error, lw=1.5, label=f'실제 최대오차: {max_error:.4f}')
        ax.axhline(bound, color='r', linestyle='--', lw=1.5,
                   label=f'베리–에센 상한: {bound:.4f}')
        ax.set_title(f'n = {n}')
        ax.set_xlabel('x')
        ax.legend(fontsize=8)
        ax.spines[['top', 'right']].set_visible(False)

    axes[0].set_ylabel(r'$|F_n(x) - \Phi(x)|$')
    plt.tight_layout()
    plt.show()

# Exp(1): 오른쪽으로 치우친 분포
berry_esseen_visualization(
    'Exponential(1)',
    lambda size: np.random.exponential(1, size),
    mu=1.0, sigma=1.0, rho=12 / np.e - 2,      # E|X-1|^3 = 12/e - 2
    sample_sizes=[5, 30, 100]
)

print(" n   상한     정확(감마)  에지워스  모의(5만)  상한/정확")
for n in (5, 30, 100):
    x = np.linspace(-np.sqrt(n), 12, 200_001)
    exact = np.abs(stats.gamma.cdf(n + x * np.sqrt(n), n)
                   - stats.norm.cdf(x)).max()
    edge = (2 / 6) * stats.norm.pdf(0) / np.sqrt(n)     # 에지워스 첫 항
    bound = 0.4748 * (12 / np.e - 2) / np.sqrt(n)

    np.random.seed(42)
    z = (np.random.exponential(1, (50_000, n)).mean(axis=1) - 1) * np.sqrt(n)
    grid = np.linspace(-4, 4, 1000)
    mc = np.abs(np.array([np.mean(z <= xx) for xx in grid])
                - stats.norm.cdf(grid)).max()

    print(f"{n:4d}  {bound:.4f}    {exact:.5f}    {edge:.5f}   {mc:.4f}"
          f"      {bound/exact:.2f}")

 n   상한     정확(감마)  에지워스  모의(5만)  상한/정확
   5  0.5127    0.05963    0.05947   0.0570      8.60
  30  0.2093    0.02429    0.02428   0.0271      8.62
 100  0.1146    0.01330    0.01330   0.0141      8.62

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.family'] = 'Apple SD Gothic Neo'
plt.rcParams['axes.unicode_minus'] = False

def convergence_rate_comparison():
    """분포마다 정규근사가 얼마나 빨리 좋아지는지 비교한다."""
    C = 0.4748           # 이 쪽의 계산에 쓰는 값. 가장 날카로운 것은 0.4690 이다(정리 1)
    ns = np.arange(5, 501)

    # 각 분포에 필요한 것은 두 값뿐이다.
    #   sigma: 표준편차
    #   rho  : 세 번째 절대적률 E|X - mu|^3
    # rho/sigma^3 이 클수록(= 비대칭하거나 꼬리가 무거울수록) 수렴이 느리다.
    distributions = {
        'Bernoulli(0.5)': {'sigma': np.sqrt(0.25), 'rho': 0.125},
        'Uniform(0,1)':   {'sigma': 1/np.sqrt(12), 'rho': 1/32},
        'Exponential(1)': {'sigma': 1.0,            'rho': 12/np.e - 2},
    }

    fig, ax = plt.subplots(figsize=(12, 4))
    for name, params in distributions.items():
        # 베리-에센 상계: |F_n(x) - Phi(x)| <= C * rho / (sigma^3 * sqrt(n))
        # n 의존성이 1/sqrt(n) 뿐이라는 점이 핵심이다.
        # 오차를 절반으로 줄이려면 표본을 네 배로 늘려야 한다.
        bounds = C * params['rho'] / (params['sigma']**3 * np.sqrt(ns))
        ax.plot(ns, bounds, label=name, lw=2)

        # 곡선의 높이를 정하는 것은 rho/sigma^3 하나뿐이다.
        # 상한을 0.1 이하로 만들려면 n >= (C rho / sigma^3 / 0.1)^2 이어야 한다.
        g = params['rho'] / params['sigma']**3
        print(f"{name:<16} rho/sigma^3 = {g:.4f}   상한 = {C*g:.4f}/sqrt(n)"
              f"   상한 <= 0.1 이 되는 n = {int(np.ceil((C * g / 0.1)**2))}")

    ax.set_xlabel('표본크기 $n$')
    ax.set_ylabel('베리–에센 상한')
    ax.set_title('정규분포로의 수렴 속도 — 베리–에센 상한')
    ax.legend()
    ax.spines[['top', 'right']].set_visible(False)
    plt.tight_layout()
    plt.show()

convergence_rate_comparison()

Bernoulli(0.5)   rho/sigma^3 = 1.0000   상한 = 0.4748/sqrt(n)   상한 <= 0.1 이 되는 n = 23
Uniform(0,1)     rho/sigma^3 = 1.2990   상한 = 0.6168/sqrt(n)   상한 <= 0.1 이 되는 n = 39
Exponential(1)   rho/sigma^3 = 2.4146   상한 = 1.1464/sqrt(n)   상한 <= 0.1 이 되는 n = 132

import numpy as np
from scipy import stats

def sup_err(n, p):                      # sup_x |F_n(x) - Phi(x)| 를 정확히
    k = np.arange(0, n + 1)
    sd = np.sqrt(p * (1 - p))
    F = stats.binom.cdf(k, n, p)
    Phi = stats.norm.cdf((k - n * p) / (sd * np.sqrt(n)))
    F_left = np.concatenate([[0.0], F[:-1]])          # 계단 직전 값
    return max(np.abs(F - Phi).max(), np.abs(F_left - Phi).max())

print("Bernoulli(0.5):  sigma^2 = 0.25,  rho = 0.125,  rho/sigma^3 = 1")
print(f"{'n':>7}{'실제 sup 오차':>17}{'베리-에센 한계':>17}{'한계/실제':>12}")
for n in (10, 30, 100, 1000, 10_000):
    actual, bound = sup_err(n, 0.5), 0.4748 / np.sqrt(n)
    print(f"{n:>7}{actual:>17.6f}{bound:>17.6f}{bound / actual:>12.2f}")

Bernoulli(0.5):  sigma^2 = 0.25,  rho = 0.125,  rho/sigma^3 = 1
      n        실제 sup 오차         베리-에센 한계       한계/실제
     10         0.123047         0.150145        1.22
     30         0.072232         0.086686        1.20
    100         0.039795         0.047480        1.19
   1000         0.012613         0.015014        1.19
  10000         0.003989         0.004748        1.19

import numpy as np
from scipy import stats

def ratio(n, p):                # 실제오차 * sqrt(n) * sigma^3 / rho
    k = np.arange(0, n + 1)
    sd = np.sqrt(p * (1 - p))
    F = stats.binom.cdf(k, n, p)
    Phi = stats.norm.cdf((k - n * p) / (sd * np.sqrt(n)))
    F_left = np.concatenate([[0.0], F[:-1]])
    err = max(np.abs(F - Phi).max(), np.abs(F_left - Phi).max())
    rho = p * (1 - p) ** 3 + (1 - p) * p ** 3
    return err * np.sqrt(n) / (rho / sd ** 3)

print(f"{'p':>7}{'rho/sigma^3':>14}"
      + "".join(f"{'n=' + str(n):>12}" for n in (100, 1000, 10_000, 100_000)))
for p in (0.5, 0.2, 0.05, 0.01):
    sd = np.sqrt(p * (1 - p))
    r = (p * (1 - p) ** 3 + (1 - p) * p ** 3) / sd ** 3
    vals = [ratio(n, p) for n in (100, 1000, 10_000, 100_000)]
    print(f"{p:>7}{r:>14.4f}" + "".join(f"{v:>12.4f}" for v in vals))

print(f"\n에세인의 하한 = (sqrt(10)+3)/(6 sqrt(2 pi)) = "
      f"{(np.sqrt(10) + 3) / (6 * np.sqrt(2 * np.pi)):.6f}")

      p   rho/sigma^3       n=100      n=1000     n=10000    n=100000
    0.5        1.0000      0.3979      0.3988      0.3989      0.3989
    0.2        1.7000      0.3498      0.3518      0.3520      0.3520
   0.05        4.1524      0.2794      0.2858      0.2865      0.2865
   0.01        9.8514      0.2393      0.2666      0.2696      0.2699

에세인의 하한 = (sqrt(10)+3)/(6 sqrt(2 pi)) = 0.409732

import numpy as np
from scipy import stats

print(f"{'n':>7}{'균등 오차':>14}{'|z|<1 구간':>15}{'|z|>2 구간':>16}"
      f"{'(1+|z|)^3 가중 최대':>22}")
for n in (100, 1000, 10_000):
    k = np.arange(0, n + 1)
    z = (k - n * 0.5) / (0.5 * np.sqrt(n))
    d = np.abs(stats.binom.cdf(k, n, 0.5) - stats.norm.cdf(z))
    print(f"{n:>7}{d.max():>14.6f}{d[np.abs(z) < 1].max():>15.6f}"
          f"{d[np.abs(z) > 2].max():>16.8f}{(d * (1 + np.abs(z)) ** 3).max():>22.6f}")

      n         균등 오차       |z|<1 구간        |z|>2 구간       (1+|z|)^3 가중 최대
    100      0.039795       0.039795      0.00369665              0.220041
   1000      0.012613       0.012613      0.00165320              0.067135
  10000      0.003989       0.003989      0.00052121              0.020971

import numpy as np
from scipy import stats

def compare(n, p):
    k = np.arange(0, n + 1)
    sd = np.sqrt(n * p * (1 - p))
    F = stats.binom.cdf(k, n, p)
    raw = np.abs(F - stats.norm.cdf((k - n * p) / sd)).max()
    cc = np.abs(F - stats.norm.cdf((k + 0.5 - n * p) / sd)).max()
    return raw, cc

print("대칭:  Binomial(n, 0.5)")
print(f"{'n':>7}{'보정 없음':>15}{'보정 있음':>15}{'개선 배수':>12}")
for n in (10, 30, 100, 1000, 10_000):
    raw, cc = compare(n, 0.5)
    print(f"{n:>7}{raw:>15.6f}{cc:>15.6f}{raw / cc:>12.1f}")

print("\n치우침:  Binomial(n, 0.1)")
print(f"{'n':>7}{'보정 없음':>15}{'보정 있음':>15}{'개선 배수':>12}")
for n in (30, 100, 1000, 10_000):
    raw, cc = compare(n, 0.1)
    print(f"{n:>7}{raw:>15.6f}{cc:>15.6f}{raw / cc:>12.1f}")

대칭:  Binomial(n, 0.5)
      n          보정 없음          보정 있음       개선 배수
     10       0.123047       0.002686        45.8
     30       0.072232       0.000918        78.7
    100       0.039795       0.000273       146.0
   1000       0.012613       0.000027       464.3
  10000       0.003989       0.000003      1469.2

치우침:  Binomial(n, 0.1)
      n          보정 없음          보정 있음       개선 배수
     30       0.147439       0.030898         4.8
    100       0.083156       0.017474         4.8
   1000       0.026599       0.005599         4.8
  10000       0.008421       0.001773         4.8