카이제곱 분포¶

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

# 자유도 5 인 카이제곱의 밀도함수와 분포함수를 겹쳐 그린다.
# 값이 0 이상에서만 정의되고 오른쪽으로 길게 늘어진 모양임을 확인한다.
df = 5
x = np.linspace(0, 30, 100)
pdf = stats.chi2(df=df).pdf(x)
cdf = stats.chi2(df=df).cdf(x)

# (1) 에서 유도한 최빈값 d-2 를 격자가 고른 최대점과 견준다.
print(f"해석적 최빈값  d - 2      = {df - 2}")
print(f"격자 최대점              = {x[np.argmax(pdf)]:.4f}  (격자 간격 {x[1] - x[0]:.4f})")
print(f"F(15) = P(Q <= 15)       = {stats.chi2(df=df).cdf(15):.4f}")
print(f"상위 5% 임계값            = {stats.chi2(df=df).ppf(0.95):.4f}")

fig, ax = plt.subplots()
ax.plot(x, pdf, label="PDF")
ax.plot(x, cdf, label="CDF")
ax.legend()
ax.set_title(f"PDF and CDF of chi-squared({df})")
plt.show()

해석적 최빈값  d - 2      = 3
격자 최대점              = 3.0303  (격자 간격 0.3030)
F(15) = P(Q <= 15)       = 0.9896
상위 5% 임계값            = 11.0705

import numpy as np
import scipy.stats as stats

df, seed = 5, 1
data_direct = stats.chi2(df=df).rvs(10_000, random_state=seed)
# 표준정규를 df개 제곱해서 더한다. 이것이 카이제곱분포의 정의다.
# axis=0으로 합해야 열마다(표본마다) 제곱합이 하나씩 나온다.
data_from_norm = np.sum(
    stats.norm().rvs(size=(df, 10_000), random_state=seed) ** 2,
    axis=0
)

for name, d in [("직접 표집", data_direct), ("제곱합 구성", data_from_norm)]:
    print(f"{name:<12} mean={d.mean():.4f}  var={d.var(ddof=1):.4f}")
print(f"{'이론값':<12} mean={df:.4f}  var={2*df:.4f}")

# 몬테카를로 오차. 평균은 SD/sqrt(N), 분산은 Var * sqrt((beta2-1)/N) 이며
# 카이제곱(d) 의 첨도는 beta2 = 3 + 12/d 다.
N, beta2 = 10_000, 3 + 12 / df
se_mean = np.sqrt(2 * df) / np.sqrt(N)
se_var = 2 * df * np.sqrt((beta2 - 1) / N)
print(f"{'MC 오차':<12} mean={se_mean:.4f}  var={se_var:.4f}")

직접 표집        mean=4.9953  var=10.0317
제곱합 구성       mean=5.0068  var=9.9999
이론값          mean=5.0000  var=10.0000
MC 오차        mean=0.0316  var=0.2098

import numpy as np
from scipy import stats

rng = np.random.default_rng(303)
print("정의 검증: ΣZ² 와 χ²(d) 가 같은 분포인가")
print(f"{'d':>4s} {'평균(이론 d)':>12s} {'분산(이론 2d)':>13s} "
      f"{'왜도(이론 √(8/d))':>16s} {'KS p':>8s}")
for d in [1, 3, 10, 30]:
    Q = (rng.normal(size=(200_000, d))**2).sum(1)
    ks = stats.kstest(Q[:5_000], "chi2", args=(d,))
    print(f"{d:4d} {Q.mean():12.4f} {Q.var(ddof=1):13.4f} "
          f"{stats.skew(Q):16.4f} {ks.pvalue:8.4f}"
          f"   (이론 {d}, {2 * d}, {np.sqrt(8 / d):.4f})")

print("\n가법성: Q1~χ²(3), Q2~χ²(7) 독립 → Q1+Q2 ~ χ²(10)?")
Q1 = rng.chisquare(3, 200_000)
Q2 = rng.chisquare(7, 200_000)
print(f"  합의 평균 {(Q1 + Q2).mean():.4f} (이론 10),  "
      f"분산 {(Q1 + Q2).var(ddof=1):.4f} (이론 20)")
print(f"  KS 검정 p = "
      f"{stats.kstest((Q1 + Q2)[:5_000], 'chi2', args=(10,)).pvalue:.4f}")
print(f"  P(Q1+Q2 > 18.31) 경험 {np.mean(Q1 + Q2 > 18.31):.4f}, "
      f"이론 {stats.chi2.sf(18.31, 10):.4f}")

정의 검증: ΣZ² 와 χ²(d) 가 같은 분포인가
   d     평균(이론 d)     분산(이론 2d)    왜도(이론 √(8/d))     KS p
   1       0.9956        1.9932           2.8438   0.0451   (이론 1, 2, 2.8284)
   3       2.9999        6.0025           1.6445   0.9745   (이론 3, 6, 1.6330)
  10       9.9882       19.9640           0.9013   0.7920   (이론 10, 20, 0.8944)
  30      30.0207       59.7564           0.5061   0.6276   (이론 30, 60, 0.5164)

가법성: Q1~χ²(3), Q2~χ²(7) 독립 → Q1+Q2 ~ χ²(10)?
  합의 평균 10.0069 (이론 10),  분산 20.0144 (이론 20)
  KS 검정 p = 0.4704
  P(Q1+Q2 > 18.31) 경험 0.0502, 이론 0.0500

import numpy as np
from scipy import stats

z = stats.norm.ppf(0.95)
print("상위 5% 임계값")
print(f"{'d':>5s} {'정확':>9s} {'정규근사':>9s} {'피셔 √2Q':>10s} {'윌슨-힐퍼티':>12s}")
for d in [5, 10, 30, 50, 100, 300]:
    ex = stats.chi2.ppf(0.95, d)
    norm = d + z * np.sqrt(2 * d)
    fisher = 0.5 * (z + np.sqrt(2 * d - 1))**2
    wh = d * (1 - 2 / (9 * d) + z * np.sqrt(2 / (9 * d)))**3
    print(f"{d:5d} {ex:9.4f} {norm:9.4f} {fisher:10.4f} {wh:12.4f}")

print("\n상대오차 (%)")
print(f"{'d':>5s} {'정규근사':>9s} {'피셔':>9s} {'윌슨-힐퍼티':>12s}")
for d in [5, 10, 30, 50, 100, 300]:
    ex = stats.chi2.ppf(0.95, d)
    norm = d + z * np.sqrt(2 * d)
    fisher = 0.5 * (z + np.sqrt(2 * d - 1))**2
    wh = d * (1 - 2 / (9 * d) + z * np.sqrt(2 / (9 * d)))**3
    print(f"{d:5d} {100 * (norm - ex) / ex:9.3f} {100 * (fisher - ex) / ex:9.3f} "
          f"{100 * (wh - ex) / ex:12.3f}")

상위 5% 임계값
    d        정확      정규근사     피셔 √2Q       윌슨-힐퍼티
    5   11.0705   10.2015    10.7873      11.0439
   10   18.3070   17.3560    18.0225      18.2918
   30   43.7730   42.7410    43.4871      43.7666
   50   67.5048   66.4485    67.2189      67.5006
  100  124.3421  123.2617   124.0563     124.3396
  300  341.3951  340.2905   341.1097     341.3940

상대오차 (%)
    d      정규근사        피셔       윌슨-힐퍼티
    5    -7.850    -2.558       -0.240
   10    -5.195    -1.554       -0.083
   30    -2.358    -0.653       -0.014
   50    -1.565    -0.424       -0.006
  100    -0.869    -0.230       -0.002
  300    -0.324    -0.084       -0.000

import numpy as np
from scipy import stats

x = np.linspace(0.5, 15, 5)

print("(가) χ²(d) = Gamma(shape=d/2, scale=2)")
print(f"{'x':>6s} {'chi2.pdf':>12s} {'gamma.pdf':>12s}")
for xi in x:
    print(f"{xi:6.2f} {stats.chi2.pdf(xi, 6):12.8f} "
          f"{stats.gamma.pdf(xi, 3, scale=2):12.8f}")

print("\n(나) χ²(2) = Exp(평균 2)")
for xi in x:
    print(f"{xi:6.2f} {stats.chi2.pdf(xi, 2):12.8f} "
          f"{stats.expon.pdf(xi, scale=2):12.8f}")

print("\n(다) Q~χ²(1) ⟺ Q = Z²")
print(f"  P(χ²(1) ≤ 3.8415) = {stats.chi2.cdf(3.8415, 1):.6f}")
print(f"  P(|Z| ≤ √3.8415)  = {2 * stats.norm.cdf(np.sqrt(3.8415)) - 1:.6f}")

print("\n(라) d2 → ∞ 이면 d1·F(d1, d2) → χ²(d1)")
print(f"  χ²(3) 의 0.95 분위수 = {stats.chi2.ppf(0.95, 3):.4f}")
print(f"{'d2':>8s} {'3·F(3,d2) 의 0.95 분위수':>24s}")
for d2 in [10, 50, 200, 10_000]:
    print(f"{d2:8d} {3 * stats.f.ppf(0.95, 3, d2):24.4f}")

print("\n(마) t(ν)² = F(1, ν),  ν → ∞ 이면 χ²(1)")
print(f"  t(30)² 의 0.95 분위수 = {stats.t.ppf(0.975, 30)**2:.4f},  "
      f"χ²(1) 의 0.95 = {stats.chi2.ppf(0.95, 1):.4f}")

(가) χ²(d) = Gamma(shape=d/2, scale=2)
     x     chi2.pdf    gamma.pdf
  0.50   0.01216876   0.01216876
  4.12   0.13520587   0.13520587
  7.75   0.07790984   0.07790984
 11.38   0.02739889   0.02739889
 15.00   0.00777775   0.00777775

(나) χ²(2) = Exp(평균 2)
  0.50   0.38940039   0.38940039
  4.12   0.06356787   0.06356787
  7.75   0.01037717   0.01037717
 11.38   0.00169403   0.00169403
 15.00   0.00027654   0.00027654

(다) Q~χ²(1) ⟺ Q = Z²
  P(χ²(1) ≤ 3.8415) = 0.950001
  P(|Z| ≤ √3.8415)  = 0.950001

(라) d2 → ∞ 이면 d1·F(d1, d2) → χ²(d1)
  χ²(3) 의 0.95 분위수 = 7.8147
      d2     3·F(3,d2) 의 0.95 분위수
      10                  11.1248
      50                   8.3700
     200                   7.9493
   10000                   7.8174

(마) t(ν)² = F(1, ν),  ν → ∞ 이면 χ²(1)
  t(30)² 의 0.95 분위수 = 4.1709,  χ²(1) 의 0.95 = 3.8415

import numpy as np
from scipy import stats

rng = np.random.default_rng(404)
B = 40_000

# (1) 표본분산: σ=2, n=8
n = 8
T1 = np.array([(n - 1) * rng.normal(0, 2, n).var(ddof=1) / 4
               for _ in range(B)])

# (2) 피어슨 적합도: 균등한 4칸, N=200
k, N = 4, 200
T2 = []
for _ in range(B // 4):
    obs = rng.multinomial(N, [1 / k] * k)
    E = N / k
    T2.append(((obs - E)**2 / E).sum())
T2 = np.array(T2)

# (3) 우도비: 알려진 분산의 정규 평균 검정 (모수 1개)
T3 = np.array([(rng.normal(0, 1, 20).mean() * np.sqrt(20))**2
               for _ in range(B)])

for lab, T, df in [("(n-1)S²/σ² (n=8)", T1, 7),
                   ("피어슨 X² (k=4)", T2, 3),
                   ("우도비 2ln Λ (모수 1개)", T3, 1)]:
    print(f"  {lab:26s} 평균 {T.mean():7.4f} (이론 {df})  "
          f"분산 {T.var(ddof=1):8.4f} (이론 {2 * df})  "
          f"KS p = {stats.kstest(T[:4_000], 'chi2', args=(df,)).pvalue:.4f}")

  (n-1)S²/σ² (n=8)           평균  7.0215 (이론 7)  분산  14.1117 (이론 14)  KS p = 0.8943
  피어슨 X² (k=4)               평균  3.0222 (이론 3)  분산   5.9384 (이론 6)  KS p = 0.0558
  우도비 2ln Λ (모수 1개)          평균  0.9801 (이론 1)  분산   1.9085 (이론 2)  KS p = 0.1084

                Gamma(d/2, 척도 2)
                       │
                    χ²(d)
                 ╱     │     ╲
       d=1 ↙       d=2 ↓       ↘ 비 χ²/χ²
      Z²        Exp(평균 2)        F(d1, d2)
                                      │ d1=1
                                      ↓
                                   t(ν)²