중심극한정리¶

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

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

# 이름 -> (표본추출기, 평균 μ, 표준편차 σ, 라벨, 이산인가)
DISTRIBUTIONS = {
    "uniform": (lambda rng, s: rng.uniform(0, 1, s),
                0.5, np.sqrt(1 / 12), "Uniform(0,1)", False),
    "exponential": (lambda rng, s: rng.exponential(1.0, s),
                    1.0, 1.0, "Exponential(1)", False),
    "lognormal": (lambda rng, s: rng.lognormal(0, 0.75, s),
                  np.exp(0.75**2 / 2),
                  np.sqrt((np.exp(0.75**2) - 1) * np.exp(0.75**2)),
                  "LogNormal(0, 0.75)", False),
    "bernoulli": (lambda rng, s: rng.binomial(1, 0.3, s).astype(float),
                  0.3, np.sqrt(0.3 * 0.7), "Bernoulli(0.3)", True),
}


def lln_and_clt(dist):
    sampler, mu, sigma, label, discrete = DISTRIBUTIONS[dist]
    rng = np.random.default_rng(2026)
    eps = EPS * sigma
    X = sampler(rng, (M, max(N_LIST)))    # 가장 큰 n 으로 한 번만 뽑는다

    fig, axes = plt.subplots(2, 10, figsize=(22, 7.5), constrained_layout=True)
    lo, hi = np.percentile(X[:, :min(N_LIST)].mean(axis=1), [0.5, 99.5])
    pad = 0.1 * (hi - lo)
    lln_bins = np.linspace(lo - pad, hi + pad, 60)
    z_grid, z_bins = np.linspace(-4, 4, 400), np.linspace(-4, 4, 50)

    for j, n in enumerate(N_LIST):
        xbar = X[:, :n].mean(axis=1)

        # 정수값 분포에서는 X̄_n 이 격자 k/n 위에만 있다.
        # 격자점마다 막대 하나를 두어야 밀도로 읽을 수 있다.
        if discrete:
            s = np.rint(X[:, :n].sum(axis=1))
            lln_bins = (np.arange(s.min(), s.max() + 2) - 0.5) / n
            z_bins = np.sqrt(n) * (lln_bins - mu) / sigma

        # 윗줄 — 배율을 주지 않으면 μ 로 오그라든다 (큰수의 법칙)
        ax = axes[0, j]
        ax.hist(xbar, bins=lln_bins, density=True, color="tab:blue",
                alpha=0.6, edgecolor="white", linewidth=0.3)
        ax.axvline(mu, color="red", lw=2)
        ax.axvspan(mu - eps, mu + eps, color="orange", alpha=0.18)
        p_out = np.mean(np.abs(xbar - mu) > eps)
        ax.set_title(f"n = {n}\n"
                     rf"$\hat P(|\bar X_n-\mu|>\varepsilon)$ = {p_out:.3f}")
        ax.set_xlim(lo - pad, hi + pad)

        # 아랫줄 — √n/σ 를 곱해 확대하면 모양이 남는다 (중심극한정리)
        ax = axes[1, j]
        z = np.sqrt(n) * (xbar - mu) / sigma
        ax.hist(z, bins=z_bins, density=True, color="tab:green",
                alpha=0.6, edgecolor="white", linewidth=0.3)
        ax.plot(z_grid, stats.norm.pdf(z_grid), "k-", lw=2)
        ax.set_title(f"n = {n}\nKS = {stats.kstest(z, 'norm').statistic:.3f},"
                     f"  skew = {stats.skew(z):.2f}")
        ax.set_xlim(-4, 4)

    fig.suptitle(f"Weak LLN (top) and CLT (bottom) for {label}")
    plt.show()


for name in DISTRIBUTIONS:
    lln_and_clt(name)

for name, (sampler, mu, sigma, label, _) in DISTRIBUTIONS.items():
    rng = np.random.default_rng(2026)
    X = sampler(rng, (M, max(N_LIST)))
    for n in (5, 50):
        xbar = X[:, :n].mean(axis=1)
        z = np.sqrt(n) * (xbar - mu) / sigma
        print(f"{label:<18} n = {n:>2}   "
              f"P_out = {np.mean(np.abs(xbar - mu) > EPS * sigma):.3f}   "
              f"sd(Z) = {z.std(ddof=1):.3f}   "
              f"skew = {stats.skew(z):+.2f}   "
              f"KS = {stats.kstest(z, 'norm').statistic:.3f}")

Uniform(0,1)       n =  5   P_out = 0.673   sd(Z) = 1.008   skew = +0.03   KS = 0.019
Uniform(0,1)       n = 50   P_out = 0.164   sd(Z) = 1.005   skew = +0.04   KS = 0.018
Exponential(1)     n =  5   P_out = 0.661   sd(Z) = 1.007   skew = +0.95   KS = 0.057
Exponential(1)     n = 50   P_out = 0.147   sd(Z) = 0.987   skew = +0.32   KS = 0.029
LogNormal(0, 0.75) n =  5   P_out = 0.632   sd(Z) = 1.004   skew = +1.62   KS = 0.083
LogNormal(0, 0.75) n = 50   P_out = 0.149   sd(Z) = 0.986   skew = +0.41   KS = 0.029
Bernoulli(0.3)     n =  5   P_out = 1.000   sd(Z) = 1.012   skew = +0.46   KS = 0.228
Bernoulli(0.3)     n = 50   P_out = 0.170   sd(Z) = 1.009   skew = +0.15   KS = 0.085

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

np.random.seed(42)
N_REPS = 2000

def sample_means(dist_rvs, sample_sizes, n_reps=N_REPS):
    """표본 크기마다 n_reps개의 표본을 뽑아 각 표본평균을 모아 돌려준다.

    dist_rvs 는 "크기 n을 받아 표본을 돌려주는 함수"다.
    이렇게 함수를 인자로 받아 두면 균등·베르누이·감마 등
    어떤 모집단에도 같은 코드를 쓸 수 있다.

    돌려주는 results[n] 이 X-bar_n 의 표본분포 근사(2000개)다.
    """
    results = {}
    for n in sample_sizes:
        means = np.array([dist_rvs(n).mean() for _ in range(n_reps)])
        results[n] = means
    return results


# 균등분포 U(0,1)로 시험해 본다. 모평균 0.5, 모분산 1/12.
# 중심극한정리는 X-bar_n 의 표준편차가 sigma/sqrt(n) 이라고 예측한다.
demo = sample_means(lambda n: np.random.uniform(0, 1, n), [2, 10, 100])
sigma = (1 / 12) ** 0.5
for n, means in demo.items():
    theory = sigma / n**0.5
    print(f"n = {n:>3}: 평균 {means.mean():.4f}  "
          f"표준편차 {means.std():.4f}  (이론 {theory:.4f},"
          f"  몬테카를로 오차 {theory / (2 * N_REPS)**0.5:.4f})")

n =   2: 평균 0.4975  표준편차 0.2050  (이론 0.2041,  몬테카를로 오차 0.0032)
n =  10: 평균 0.5022  표준편차 0.0911  (이론 0.0913,  몬테카를로 오차 0.0014)
n = 100: 평균 0.5002  표준편차 0.0288  (이론 0.0289,  몬테카를로 오차 0.0005)

sample_sizes = [2, 10, 100]

# 모양이 서로 전혀 다른 모집단 셋을 준비한다.
# 각 항목은 표본을 뽑는 함수(rvs)와 모집단 밀도를 그릴 정보를 담는다.
#   Uniform: 평평하다        (봉우리가 없음)
#   Beta   : 왼쪽으로 치우침  (유계)
#   Gamma  : 오른쪽으로 치우침 (유계가 아님)
# 세 모집단이 이렇게 다른데도 표본평균은 모두 종 모양으로 간다는 것이 요점이다.
distributions = {
    "Uniform(2, 8)": {
        "rvs": lambda n: np.random.uniform(2, 8, n),
        "color": "tomato",
        "pop_x": np.linspace(2, 8, 200),
        "pop_pdf": lambda x: np.ones_like(x) / 6,
        "sd": 3 ** 0.5, "skew": 0.0,
    },
    "Beta(6, 2)": {
        "rvs": lambda n: stats.beta.rvs(6, 2, size=n),
        "color": "seagreen",
        "pop_x": np.linspace(0, 1, 200),
        "pop_pdf": lambda x: stats.beta.pdf(x, 6, 2),
        "sd": stats.beta.std(6, 2), "skew": float(stats.beta.stats(6, 2, moments="s")),
    },
    "Gamma(6, 1)": {
        "rvs": lambda n: stats.gamma.rvs(6, size=n),
        "color": "steelblue",
        "pop_x": np.linspace(0, 25, 200),
        "pop_pdf": lambda x: stats.gamma.pdf(x, 6),
        "sd": 6 ** 0.5, "skew": 2 / 6 ** 0.5,
    },
}

n_dists = len(distributions)
# 격자 구성: 열 = 모집단, 행 = (모집단 자체, n=2, n=10, n=100)
# 세로로 내려가며 읽으면 "n이 커질수록 어떻게 변하는가"가 보이고,
# 가로로 읽으면 "모집단이 달라도 결과가 같은가"가 보인다.
n_rows = 1 + len(sample_sizes)
fig, axes = plt.subplots(n_rows, n_dists, figsize=(6 * n_dists, 4 * n_rows))

for col, (name, d) in enumerate(distributions.items()):
    c = d["color"]

    # 0행: 모집단의 밀도함수. 셋이 얼마나 다른지 먼저 확인한다.
    ax = axes[0, col]
    ax.plot(d["pop_x"], d["pop_pdf"](d["pop_x"]), lw=3, color=c)
    ax.fill_between(d["pop_x"], d["pop_pdf"](d["pop_x"]), alpha=0.3, color=c)
    ax.set_title(name, fontsize=14, fontweight="bold")
    if col == 0:
        ax.set_ylabel("Population PDF", fontsize=11)

    # 1~3행: 표본평균의 표집분포.
    # 앞서 정의한 sample_means 로 각 n마다 2000개의 표본평균을 얻는다.
    means_dict = sample_means(d["rvs"], sample_sizes)
    for row, n in enumerate(sample_sizes, start=1):
        ax = axes[row, col]
        ax.hist(means_dict[n], bins=30, color=c, alpha=0.5,
                edgecolor="white", density=True)
        ax.set_title(f"n = {n}", fontsize=11)
        if col == 0:
            ax.set_ylabel(f"Sampling Dist (n={n})", fontsize=10)

        # 폭과 모양을 이론값과 나란히 적어 둔다.
        # 폭은 sigma/sqrt(n) 으로 등식이고, 왜도는 gamma/sqrt(n) 으로 줄어든다.
        m = means_dict[n]
        print(f"{name:<14} n = {n:>3}:  표준편차 {m.std(ddof=1):.4f}"
              f" (이론 {d['sd'] / n**0.5:.4f})"
              f"   왜도 {stats.skew(m):+.3f} (이론 {d['skew'] / n**0.5:+.3f})")

plt.suptitle("Central Limit Theorem: Sampling Distribution of x̄",
             fontsize=15, y=1.01)
plt.tight_layout()
plt.show()

Uniform(2, 8)  n =   2:  표준편차 1.2230 (이론 1.2247)   왜도 -0.074 (이론 +0.000)
Uniform(2, 8)  n =  10:  표준편차 0.5580 (이론 0.5477)   왜도 +0.012 (이론 +0.000)
Uniform(2, 8)  n = 100:  표준편차 0.1748 (이론 0.1732)   왜도 +0.055 (이론 +0.000)
Beta(6, 2)     n =   2:  표준편차 0.1058 (이론 0.1021)   왜도 -0.578 (이론 -0.490)
Beta(6, 2)     n =  10:  표준편차 0.0463 (이론 0.0456)   왜도 -0.189 (이론 -0.219)
Beta(6, 2)     n = 100:  표준편차 0.0146 (이론 0.0144)   왜도 -0.197 (이론 -0.069)
Gamma(6, 1)    n =   2:  표준편차 1.7247 (이론 1.7321)   왜도 +0.511 (이론 +0.577)
Gamma(6, 1)    n =  10:  표준편차 0.7902 (이론 0.7746)   왜도 +0.181 (이론 +0.258)
Gamma(6, 1)    n = 100:  표준편차 0.2404 (이론 0.2449)   왜도 +0.032 (이론 +0.082)

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

rng = np.random.default_rng(0)
dists = [("Uniform(0,1)", stats.uniform), ("Exponential(1)", stats.expon),
         ("Lognormal(0,1)", stats.lognorm(1.0)),
         ("Bernoulli(0.05)", stats.bernoulli(0.05))]
ns = (5, 30, 100, 1000)

print("표준화 표본평균과 N(0,1) 사이의 콜모고로프-스미르노프 거리")
print(f"{'분포':>18}{'왜도':>8}" + "".join(f"{'n=' + str(n):>11}" for n in ns))
ks = {}
for name, d in dists:
    mu, sd = d.mean(), d.std()
    ks[name] = [stats.kstest(
        np.sqrt(n) * (d.rvs(size=(100_000, n), random_state=rng).mean(1) - mu) / sd,
        "norm").statistic for n in ns]
    print(f"{name:>18}{float(d.stats(moments='s')):>8.2f}"
          + "".join(f"{v:>11.4f}" for v in ks[name]))

print(f"\n참고: 100000회 모의에서 KS 통계량 자체의 크기 ~ 0.83/sqrt(B) = "
      f"{0.83 / np.sqrt(100_000):.4f}")
print("\nKS ~ c/sqrt(n) 으로 보고 (n=100 기준) KS < 0.01 에 필요한 n")
print(f"{'분포':>18}{'c':>10}{'필요한 n':>12}")
for name, row in ks.items():
    c = row[2] * 10
    print(f"{name:>18}{c:>10.3f}{(c / 0.01) ** 2:>12.0f}")

plt.figure(figsize=(8, 6))
for name, row in ks.items():
    plt.loglog(ns, row, "o-", label=name)
plt.loglog(ns, [0.3 / np.sqrt(n) for n in ns], "k--", label="reference $n^{-1/2}$")
plt.axhline(0.01, color="gray", lw=1)
plt.xlabel("sample size n")
plt.ylabel("KS distance to N(0,1)")
plt.title("How fast does the CLT actually converge?")
plt.legend()
plt.grid(alpha=0.3, which="both")
plt.show()

표준화 표본평균과 N(0,1) 사이의 콜모고로프-스미르노프 거리
                분포      왜도        n=5       n=30      n=100     n=1000
      Uniform(0,1)    0.00     0.0079     0.0024     0.0020     0.0017
    Exponential(1)    2.00     0.0603     0.0243     0.0138     0.0036
    Lognormal(0,1)    6.18     0.1156     0.0590     0.0365     0.0110
   Bernoulli(0.05)    4.13     0.4716     0.2134     0.1165     0.0381

참고: 100000회 모의에서 KS 통계량 자체의 크기 ~ 0.83/sqrt(B) = 0.0026

KS ~ c/sqrt(n) 으로 보고 (n=100 기준) KS < 0.01 에 필요한 n
                분포         c       필요한 n
      Uniform(0,1)     0.020           4
    Exponential(1)     0.138         189
    Lognormal(0,1)     0.365        1333
   Bernoulli(0.05)     1.165       13563

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 30
Z = np.sqrt(n) * (rng.exponential(1.0, (4_000_000, n)).mean(1) - 1)

print(f"Exp(1), n = {n}.  Z = sqrt(n)(Xbar - 1) 를 N(0,1) 로 근사")
print(f"{'z':>4}{'참 P(Z>z)':>14}{'정규 근사':>14}{'절대오차':>12}{'상대오차':>11}")
for z in (0, 1, 2, 3, 4):
    emp, nor = np.mean(Z > z), stats.norm.sf(z)
    print(f"{z:>4}{emp:>14.6f}{nor:>14.6f}{abs(emp - nor):>12.6f}"
          f"{abs(emp - nor) / nor:>10.1%}")

print("\n왼쪽 꼬리")
print(f"{'z':>4}{'참 P(Z<-z)':>14}{'정규 근사':>14}{'절대오차':>12}{'상대오차':>11}")
for z in (1, 2, 3, 4):
    emp, nor = np.mean(Z < -z), stats.norm.sf(z)
    print(f"{z:>4}{emp:>14.6f}{nor:>14.6f}{abs(emp - nor):>12.6f}"
          f"{abs(emp - nor) / nor:>10.1%}")

Exp(1), n = 30.  Z = sqrt(n)(Xbar - 1) 를 N(0,1) 로 근사
   z      참 P(Z>z)         정규 근사        절대오차       상대오차
   0      0.475818      0.500000    0.024182      4.8%
   1      0.157599      0.158655    0.001056      0.7%
   2      0.031648      0.022750    0.008897     39.1%
   3      0.004181      0.001350    0.002831    209.7%
   4      0.000365      0.000032    0.000334   1053.3%

왼쪽 꼬리
   z     참 P(Z<-z)         정규 근사        절대오차       상대오차
   1      0.156857      0.158655    0.001799      1.1%
   2      0.012149      0.022750    0.010601     46.6%
   3      0.000079      0.001350    0.001270     94.1%
   4      0.000000      0.000032    0.000032    100.0%

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

rng = np.random.default_rng(0)
dists = [("Uniform(0,1)", stats.uniform), ("Exponential(1)", stats.expon),
         ("Lognormal(0,1)", stats.lognorm(1.0))]
ns = [5, 30, 100]

fig, axes = plt.subplots(len(dists), len(ns), figsize=(12, 11))
for i, (name, d) in enumerate(dists):
    mu, sd = d.mean(), d.std()
    for j, n in enumerate(ns):
        Z = np.sqrt(n) * (d.rvs(size=(4000, n), random_state=rng).mean(1) - mu) / sd
        ax = axes[i, j]
        stats.probplot(Z, dist="norm", plot=ax)
        ax.get_lines()[0].set_markersize(2)
        ax.get_lines()[1].set_color("red")
        ax.set_title(f"{name},  n = {n}", fontsize=10)
        ax.set_xlabel("Theoretical quantiles" if i == len(dists) - 1 else "")
        ax.set_ylabel("Sample quantiles" if j == 0 else "")
        ax.set_xlim(-4, 4)
        ax.set_ylim(-4, 6)
plt.suptitle("Q-Q plots of standardized sample means", fontsize=14)
plt.tight_layout()
plt.show()

print(f"{'분포':>16}{'n':>6}{'표준화 평균의 왜도':>20}{'P(Z>3)':>11}"
      f"{'정규값':>10}")
for name, d in dists:
    mu, sd = d.mean(), d.std()
    for n in ns:
        Z = np.sqrt(n) * (d.rvs(size=(200_000, n), random_state=rng).mean(1) - mu) / sd
        print(f"{name:>16}{n:>6}{stats.skew(Z):>20.3f}"
              f"{np.mean(Z > 3):>11.5f}{0.00135:>10.5f}")

              분포     n          표준화 평균의 왜도     P(Z>3)       정규값
    Uniform(0,1)     5               0.001    0.00047   0.00135
    Uniform(0,1)    30               0.010    0.00128   0.00135
    Uniform(0,1)   100              -0.008    0.00131   0.00135
  Exponential(1)     5               0.886    0.00935   0.00135
  Exponential(1)    30               0.359    0.00411   0.00135
  Exponential(1)   100               0.212    0.00288   0.00135
  Lognormal(0,1)     5               2.806    0.01650   0.00135
  Lognormal(0,1)    30               1.106    0.01041   0.00135
  Lognormal(0,1)   100               0.626    0.00684   0.00135

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
print("모집단이 N(0,1) 일 때 정규분포까지의 KS 거리")
print(f"{'n':>7}{'표본평균':>15}{'표본분산':>15}")
for n in (5, 10, 30, 100, 1000):
    X = rng.normal(0, 1, (200_000, n))
    Zm = np.sqrt(n) * X.mean(1)
    Zv = np.sqrt(n) * (X.var(1, ddof=1) - 1) / np.sqrt(2)
    print(f"{n:>7}{stats.kstest(Zm, 'norm').statistic:>15.4f}"
          f"{stats.kstest(Zv, 'norm').statistic:>15.4f}")

print("\n표본분산은 (X - mu)^2 들의 평균이다. 그 분포의 왜도를 보라.")
print(f"  X ~ N(0,1) 의 왜도          {float(stats.norm.stats(moments='s')):>8.4f}")
print(f"  (X-mu)^2 ~ chi^2_1 의 왜도  {float(stats.chi2.stats(1, moments='s')):>8.4f}"
      f"   (= sqrt(8))")
print(f"  참고: Exp(1) 의 왜도        {float(stats.expon.stats(moments='s')):>8.4f}")

모집단이 N(0,1) 일 때 정규분포까지의 KS 거리
      n           표본평균           표본분산
      5         0.0021         0.1002
     10         0.0012         0.0675
     30         0.0017         0.0372
    100         0.0024         0.0183
   1000         0.0015         0.0052

표본분산은 (X - mu)^2 들의 평균이다. 그 분포의 왜도를 보라.
  X ~ N(0,1) 의 왜도            0.0000
  (X-mu)^2 ~ chi^2_1 의 왜도    2.8284   (= sqrt(8))
  참고: Exp(1) 의 왜도          2.0000

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
print("자료 ~ Exp(1)  (mu = 1, sigma = 1)")
print(f"{'n':>6}{'sigma 사용 sd':>16}{'s 사용 sd':>14}"
      f"{'|T|>1.96 비율':>16}{'|T|>t 분위수':>15}")
for n in (5, 10, 30, 100, 1000):
    X = rng.exponential(1.0, (200_000, n))
    m, s = X.mean(1), X.std(1, ddof=1)
    Z = np.sqrt(n) * (m - 1) / 1.0
    T = np.sqrt(n) * (m - 1) / s
    tq = stats.t.ppf(0.975, n - 1)
    print(f"{n:>6}{Z.std():>16.4f}{T.std():>14.4f}"
          f"{np.mean(np.abs(T) > 1.96):>16.4f}{np.mean(np.abs(T) > tq):>15.4f}")
print("  (마지막 두 열의 목표는 0.05)")

자료 ~ Exp(1)  (mu = 1, sigma = 1)
     n     sigma 사용 sd       s 사용 sd     |T|>1.96 비율      |T|>t 분위수
     5          0.9989        2.3325          0.1882         0.1170
    10          0.9974        1.4935          0.1305         0.0990
    30          1.0015        1.1474          0.0824         0.0729
   100          1.0021        1.0459          0.0608         0.0582
  1000          1.0005        1.0046          0.0512         0.0510
  (마지막 두 열의 목표는 0.05)

import numpy as np

rng = np.random.default_rng(0)
print("반례:  P(X_n = 0) = 1 - 1/n,  P(X_n = n) = 1/n")
print(f"{'n':>7}{'P(X_n=0)':>12}{'E[X_n]':>10}{'E[X_n^2]':>12}")
for n in (10, 100, 1000, 10_000):
    print(f"{n:>7}{1 - 1 / n:>12.4f}{n * (1 / n):>10.2f}{n * n * (1 / n):>12.1f}")
print("  X_n -> 0 (분포수렴)  그런데 E[X_n] = 1 이고 E[X_n^2] -> 무한\n")

print("통계적 사례: Exp 비율의 최대가능도추정량 1/Xbar  (참값 1)")
print(f"{'n':>6}{'점근 sd = 1/sqrt(n)':>22}{'실제 sd':>12}{'실제 평균':>12}")
for n in (3, 5, 10, 30, 100):
    r = 1 / rng.exponential(1.0, (400_000, n)).mean(1)
    print(f"{n:>6}{1 / np.sqrt(n):>22.4f}{r.std():>12.4f}{r.mean():>12.4f}")

반례:  P(X_n = 0) = 1 - 1/n,  P(X_n = n) = 1/n
      n    P(X_n=0)    E[X_n]    E[X_n^2]
     10      0.9000      1.00        10.0
    100      0.9900      1.00       100.0
   1000      0.9990      1.00      1000.0
  10000      0.9999      1.00     10000.0
  X_n -> 0 (분포수렴)  그런데 E[X_n] = 1 이고 E[X_n^2] -> 무한

통계적 사례: Exp 비율의 최대가능도추정량 1/Xbar  (참값 1)
     n     점근 sd = 1/sqrt(n)       실제 sd       실제 평균
     3                0.5774      1.4360      1.5000
     5                0.4472      0.7238      1.2505
    10                0.3162      0.3937      1.1114
    30                0.1826      0.1955      1.0346
   100                0.1000      0.1019      1.0100

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n, B = 2000, 40_000
cases = [("정규", stats.norm, 0.50), ("정규", stats.norm, 0.90),
         ("라플라스", stats.laplace, 0.50),
         ("지수", stats.expon, 0.50), ("지수", stats.expon, 0.95)]

print(f"{'분포':>10}{'p':>7}{'q_p':>9}{'f(q_p)':>9}"
      f"{'이론 점근분산':>16}{'모의 n*Var':>14}")
for name, dist, p in cases:
    q, f = dist.ppf(p), dist.pdf(dist.ppf(p))
    qh = np.quantile(dist.rvs(size=(B, n), random_state=rng), p, axis=1)
    print(f"{name:>10}{p:>7.2f}{q:>9.4f}{f:>9.4f}"
          f"{p * (1 - p) / f ** 2:>16.4f}{n * qh.var():>14.4f}")

print("\n중앙값과 평균의 점근분산 비교")
for name, dist in [("정규", stats.norm), ("라플라스", stats.laplace)]:
    f = dist.pdf(dist.ppf(0.5))
    print(f"  {name:>8}: 중앙값 {0.25 / f ** 2:>7.4f}   평균 {dist.var():>7.4f}"
          f"   비 {0.25 / f ** 2 / dist.var():>6.4f}")

        분포      p      q_p   f(q_p)         이론 점근분산      모의 n*Var
        정규   0.50   0.0000   0.3989          1.5708        1.5712
        정규   0.90   1.2816   0.1755          2.9221        2.9198
      라플라스   0.50   0.0000   0.5000          1.0000        1.0507
        지수   0.50   0.6931   0.5000          1.0000        0.9987
        지수   0.95   2.9957   0.0500         19.0000       18.7912

중앙값과 평균의 점근분산 비교
        정규: 중앙값  1.5708   평균  1.0000   비 1.5708
      라플라스: 중앙값  1.0000   평균  2.0000   비 0.5000