Shapiro-Wilk 검정력 모의실험¶

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

def power_for_n(n, sims=500, sigma_ln=0.6, alpha=0.05, seed=0):
    """주어진 표본크기에서 Shapiro-Wilk 검정의 검정력을 모의실험으로 잰다.

    대립가설 쪽 모집단은 로그정규로 고정한다. 정규에서 벗어난 정도는
    내내 같으므로, 검정력의 변화는 오직 표본크기 때문이다.
    """
    rng = np.random.default_rng(seed)
    rejections = 0
    for _ in range(sims):
        x = rng.lognormal(mean=0.0, sigma=sigma_ln, size=n)
        W, p = stats.shapiro(x)
        if p < alpha:
            rejections += 1
    return rejections / sims

# n 을 키워 가며 검정력 곡선을 그린다. 같은 이탈인데도 n=20 에서는
# 절반도 잡아내지 못하고 n=300 에서는 거의 놓치지 않는다.
ns = [20, 30, 50, 80, 120, 200, 300]
powers = [power_for_n(n, sims=400, sigma_ln=0.6, alpha=0.05,
                      seed=42 + n) for n in ns]

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(ns, powers, marker="o")
ax.set_ylim(0, 1)
ax.set_xlabel("Sample size (n)")
ax.set_ylabel("Empirical power (alpha = 0.05)")
ax.set_title("Shapiro-Wilk power vs n (lognormal alt, sigma = 0.6)")
plt.tight_layout()
plt.show()

for n, pw in zip(ns, powers):
    print(f"n = {n:>3}: power = {pw:.3f}")

n =  20: power = 0.647
n =  30: power = 0.833
n =  50: power = 0.978
n =  80: power = 1.000
n = 120: power = 1.000
n = 200: power = 1.000
n = 300: power = 1.000

import numpy as np
from scipy import stats

sigma_ln, alpha = 0.6, 0.05

# 1) 대립가설의 이탈 크기는 닫힌 꼴로 정해져 있다.
m = np.exp(sigma_ln ** 2)
g1 = (m + 2) * np.sqrt(m - 1)
g2 = m ** 4 + 2 * m ** 3 + 3 * m ** 2 - 6
print(f"Lognormal(0, {sigma_ln}):  m = exp(sigma^2) = {m:.5f}")
print(f"  왜도      (m+2) sqrt(m-1)        = {g1:.4f}")
print(f"  초과첨도  m^4 + 2m^3 + 3m^2 - 6  = {g2:.4f}")
lg = stats.lognorm(s=sigma_ln)
print(f"  scipy 확인:  왜도 {lg.stats('s'):.4f}   초과첨도 {lg.stats('k'):.4f}")

# 2) 반복 400 회가 주는 정밀도. 기각 횟수까지 함께 본다.
print(f"\n보기의 표를 기각 횟수와 몬테카를로 오차까지 적으면 (sims = 400)")
sims = 400
for n in (20, 30, 50, 80, 120, 200, 300):
    rng = np.random.default_rng(42 + n)
    rej = sum(1 for _ in range(sims)
              if stats.shapiro(rng.lognormal(0.0, sigma_ln, n)).pvalue < alpha)
    ph = rej / sims
    se = np.sqrt(ph * (1 - ph) / sims)
    if rej == sims:
        lower = alpha ** (1 / sims)       # 클로퍼-피어슨 95% 하한
        note = f"-> 95% 하한 {lower:.4f}"
    else:
        note = f"+- {1.96 * se:.4f}  (MC SE {se:.4f})"
    print(f"  n = {n:4d}  기각 {rej:3d}/{sims}  검정력 {ph:.4f}  {note}")

# 3) 0.8 선을 더 촘촘히, 반복을 늘려 좁힌다.
print("\n0.8 선 좁히기 (sims = 20000)")
for n in (22, 24, 26, 28, 30):
    rng = np.random.default_rng(7000 + n)
    rej = sum(1 for _ in range(20000)
              if stats.shapiro(rng.lognormal(0.0, sigma_ln, n)).pvalue < alpha)
    print(f"  n = {n:4d}  검정력 {rej / 20000:.4f}  (MC SE "
          f"{np.sqrt(rej / 20000 * (1 - rej / 20000) / 20000):.4f})")

Lognormal(0, 0.6):  m = exp(sigma^2) = 1.43333
  왜도      (m+2) sqrt(m-1)        = 2.2601
  초과첨도  m^4 + 2m^3 + 3m^2 - 6  = 10.2734
  scipy 확인:  왜도 2.2601   초과첨도 10.2734

보기의 표를 기각 횟수와 몬테카를로 오차까지 적으면 (sims = 400)
  n =   20  기각 259/400  검정력 0.6475  +- 0.0468  (MC SE 0.0239)
  n =   30  기각 333/400  검정력 0.8325  +- 0.0366  (MC SE 0.0187)
  n =   50  기각 391/400  검정력 0.9775  +- 0.0145  (MC SE 0.0074)
  n =   80  기각 400/400  검정력 1.0000  -> 95% 하한 0.9925
  n =  120  기각 400/400  검정력 1.0000  -> 95% 하한 0.9925
  n =  200  기각 400/400  검정력 1.0000  -> 95% 하한 0.9925
  n =  300  기각 400/400  검정력 1.0000  -> 95% 하한 0.9925

0.8 선 좁히기 (sims = 20000)
  n =   22  검정력 0.6968  (MC SE 0.0033)
  n =   24  검정력 0.7439  (MC SE 0.0031)
  n =   26  검정력 0.7819  (MC SE 0.0029)
  n =   28  검정력 0.8154  (MC SE 0.0027)
  n =   30  검정력 0.8436  (MC SE 0.0026)

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

def power_curve(sigma, ns, sims=500, alpha=0.05):
    powers = []
    for n in ns:
        rng = np.random.default_rng(42 + n)
        rej = sum(1 for _ in range(sims)
                  if stats.shapiro(rng.lognormal(0, sigma, n))[1] < alpha)
        powers.append(rej / sims)
    return powers

ns = [20, 30, 50, 80, 120, 200, 300, 500]
p03 = power_curve(0.3, ns)
p06 = power_curve(0.6, ns)

for n, a, b in zip(ns, p03, p06):
    print(f"n = {n:>3}: sigma=0.3 -> {a:.3f}, sigma=0.6 -> {b:.3f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(ns, p03, marker="o", label="sigma = 0.3")
ax.plot(ns, p06, marker="s", label="sigma = 0.6")
ax.axhline(0.8, color="gray", linestyle="--", alpha=0.5)
ax.set_xlabel("Sample size")
ax.set_ylabel("Power")
ax.set_title("Power Curves: Lognormal Alternatives")
ax.legend()
plt.tight_layout()
plt.show()

n =  20: sigma=0.3 -> 0.246, sigma=0.6 -> 0.658
n =  30: sigma=0.3 -> 0.364, sigma=0.6 -> 0.830
n =  50: sigma=0.3 -> 0.540, sigma=0.6 -> 0.980
n =  80: sigma=0.3 -> 0.786, sigma=0.6 -> 1.000
n = 120: sigma=0.3 -> 0.936, sigma=0.6 -> 1.000
n = 200: sigma=0.3 -> 0.998, sigma=0.6 -> 1.000
n = 300: sigma=0.3 -> 0.998, sigma=0.6 -> 1.000
n = 500: sigma=0.3 -> 1.000, sigma=0.6 -> 1.000

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

ns = [20, 30, 50, 80, 120, 200, 300]
sims, alpha, sigma = 500, 0.05, 0.6
pw_sw, pw_k2 = [], []

for n in ns:
    rng = np.random.default_rng(42 + n)
    r_sw = r_k2 = 0
    for _ in range(sims):
        x = rng.lognormal(0, sigma, n)
        if stats.shapiro(x)[1] < alpha:
            r_sw += 1
        if stats.normaltest(x)[1] < alpha:
            r_k2 += 1
    pw_sw.append(r_sw / sims)
    pw_k2.append(r_k2 / sims)
    print(f"n = {n:>3}: SW = {r_sw/sims:.3f}, K2 = {r_k2/sims:.3f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(ns, pw_sw, marker="o", label="Shapiro-Wilk")
ax.plot(ns, pw_k2, marker="s", label="D'Agostino K^2")
ax.set_xlabel("Sample size")
ax.set_ylabel("Power")
ax.legend()
plt.tight_layout()
plt.show()

n =  20: SW = 0.658, K2 = 0.552
n =  30: SW = 0.830, K2 = 0.702
n =  50: SW = 0.980, K2 = 0.912
n =  80: SW = 1.000, K2 = 0.990
n = 120: SW = 1.000, K2 = 1.000
n = 200: SW = 1.000, K2 = 1.000
n = 300: SW = 1.000, K2 = 1.000

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

rng = np.random.default_rng(0)
n, sims, alpha = 100, 1000, 0.05
dfs = [3, 5, 10, 30, 100]
powers = []

for df in dfs:
    rej = sum(1 for _ in range(sims)
              if stats.shapiro(rng.standard_t(df, n))[1] < alpha)
    powers.append(rej / sims)
    print(f"df = {df:>3}: power = {rej/sims:.3f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(dfs, powers, marker="o")
ax.set_xscale("log")
ax.axhline(0.05, color="red", linestyle="--", alpha=0.5,
           label="Nominal alpha")
ax.set_xlabel("Degrees of freedom (nu, log scale)")
ax.set_ylabel("Power")
ax.set_title("SW Power vs t(nu) at n = 100")
ax.legend()
plt.tight_layout()
plt.show()

df =   3: power = 0.878
df =   5: power = 0.567
df =  10: power = 0.263
df =  30: power = 0.071
df = 100: power = 0.048