카이제곱분포¶

import numpy as np
from scipy import optimize, special, stats

# (1) q -> 0+ 에서의 거동. 지수 d/2 - 1 의 부호가 세 갈래를 만든다.
print(f"{'d':>4}{'지수 d/2-1':>12}{'f(1e-6)':>14}{'f(0.01)':>12}")
for d in (1, 2, 3, 4, 8):
    print(f"{d:>4}{d / 2 - 1:>12.1f}{stats.chi2(d).pdf(1e-6):>14.6g}"
          f"{stats.chi2(d).pdf(0.01):>12.6g}")

# d=1 은 밀도가 발산하지만 질량은 0 으로 간다.
print(f"\nd=1:  f(0.01) = {stats.chi2(1).pdf(0.01):.4f}  <- 그림의 ylim 0.6 을 훌쩍 넘는다")
print(f"      P(Q <= 0.01) = {stats.chi2(1).cdf(0.01):.6f}"
      f"  = 2*Phi(0.1) - 1 = {2 * stats.norm.cdf(0.1) - 1:.6f}")

# (2) 교차점의 닫힌 꼴. 비가 q 의 양의 거듭제곱이므로 영점은 하나뿐이다.
def q_star(d1, d2):
    return 2 * (special.gamma(d2 / 2) / special.gamma(d1 / 2)) ** (2 / (d2 - d1))

print(f"\n{'(d1,d2)':>10}{'닫힌 꼴':>14}{'수치해':>14}{'두 밀도 값':>26}")
for d1, d2 in [(1, 2), (2, 4), (4, 8), (1, 8)]:
    q = q_star(d1, d2)
    root = optimize.brentq(
        lambda x: stats.chi2(d2).pdf(x) - stats.chi2(d1).pdf(x), 1e-9, 200)
    print(f"{f'({d1},{d2})':>10}{q:>14.9f}{root:>14.9f}"
          f"{stats.chi2(d1).pdf(q):>13.9f}{stats.chi2(d2).pdf(q):>13.9f}")

print(f"\n손으로 얻은 값:  2/pi = {2 / np.pi:.9f},   2 = 2,   2*sqrt(6) = {2 * np.sqrt(6):.9f}")

   d    지수 d/2-1       f(1e-6)     f(0.01)
   1        -0.5       398.942     3.96953
   2         0.0           0.5    0.497506
   3         0.5   0.000398942   0.0396953
   4         1.0       2.5e-07  0.00248753
   8         3.0   1.04167e-20 1.03647e-08

d=1:  f(0.01) = 3.9695  <- 그림의 ylim 0.6 을 훌쩍 넘는다
      P(Q <= 0.01) = 0.079656  = 2*Phi(0.1) - 1 = 0.079656

   (d1,d2)          닫힌 꼴           수치해                    두 밀도 값
     (1,2)   0.636619772   0.636619772  0.363688675  0.363688675
     (2,4)   2.000000000   2.000000000  0.183939721  0.183939721
     (4,8)   4.898979486   4.898979486  0.105741569  0.105741569
     (1,8)   2.833593280   2.833593280  0.057469093  0.057469093

손으로 얻은 값:  2/pi = 0.636619772,   2 = 2,   2*sqrt(6) = 4.898979486

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

x = np.linspace(0.01, 20, 400)   # 0에서 시작하면 d=1에서 발산해 그림이 깨진다

fig, ax = plt.subplots(figsize=(12, 3))
# 자유도 = 더한 제곱의 개수.
#   d=1 : 0 근처에서 치솟는다(정규분포가 0 근처에 몰려 있으므로)
#   d=2 : 지수분포 Exp(1/2)와 정확히 같다. 0에서 높이 0.5로 시작한다.
#   d>=3: 최빈값 d-2 에 봉우리가 생기고, 평균 d 는 그보다 오른쪽에 있다.
for d in [1, 2, 4, 8]:
    ax.plot(x, stats.chi2(d).pdf(x), lw=2, label=f'd={d}')
ax.axhline(0.5, color='gray', ls=':', lw=1)   # d=2가 0에서 닿는 높이
ax.set_ylim(0, 0.6)
ax.set_xlabel('q')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

import numpy as np
from scipy import integrate, stats

d = 5
chi2 = stats.chi2(d)

# (1) 최빈값을 격자로 확인한다. 유도한 값은 d - 2 = 3 이다.
grid = np.linspace(1e-9, 20, 2_000_001)
print(f"d = {d}:  유도한 최빈값 d-2 = {d - 2},"
      f"  격자 최대점 = {grid[chi2.pdf(grid).argmax()]:.6f}")

# 세 중심의 순서와 간격. 평균 - 최빈값은 자유도와 무관하게 2 다.
print(f"\n{'d':>5}{'최빈값 d-2':>12}{'중앙값':>12}{'WH 근사':>12}{'평균 d':>9}"
      f"{'평균-중앙':>11}{'중앙-최빈':>11}{'비':>7}")
for k in (3, 5, 10, 20, 50, 100, 500):
    med = stats.chi2(k).median()
    wh = k * (1 - 2 / (9 * k)) ** 3
    print(f"{k:>5}{k - 2:>12}{med:>12.6f}{wh:>12.6f}{k:>9}"
          f"{k - med:>11.6f}{med - (k - 2):>11.6f}{(med - (k - 2)) / (k - med):>7.3f}")

# 치우침을 표준편차 단위로 재면 평균-최빈값 간격이 왜도의 절반이다.
print(f"\n{'d':>5}{'(평균-최빈)/sd':>16}{'sqrt(2/d)':>12}{'왜도/2':>10}{'P(Q<=d)':>11}")
for k in (3, 5, 10, 20, 50, 100, 500):
    print(f"{k:>5}{2 / np.sqrt(2 * k):>16.6f}{np.sqrt(2 / k):>12.6f}"
          f"{np.sqrt(8 / k) / 2:>10.6f}{stats.chi2(k).cdf(k):>11.6f}")

# d=5 에서 세 자리의 밀도 높이와, 평균 왼쪽에 놓인 확률
print(f"\nd = 5:  f(최빈 3) = {chi2.pdf(3):.6f},  f(중앙 {chi2.median():.4f})"
      f" = {chi2.pdf(chi2.median()):.6f},  f(평균 5) = {chi2.pdf(5):.6f}")
print(f"        P(Q <= 5) = {chi2.cdf(5):.6f}"
      f"   (quad 로 재확인 {integrate.quad(chi2.pdf, 0, 5)[0]:.6f})")

d = 5:  유도한 최빈값 d-2 = 3,  격자 최대점 = 3.000000

    d     최빈값 d-2         중앙값       WH 근사     평균 d      평균-중앙      중앙-최빈      비
    3           1    2.365974    2.381497        3   0.634026   1.365974  2.154
    5           3    4.351460    4.362524        5   0.648540   1.351460  2.084
   10           8    9.341818    9.348038       10   0.658182   1.341818  2.039
   20          18   19.337429   19.340713       20   0.662571   1.337429  2.019
   50          48   49.334937   49.336292       50   0.665063   1.334937  2.007
  100          98   99.334129   99.334814      100   0.665871   1.334129  2.004
  500         498  499.333492  499.333630      500   0.666508   1.333492  2.001

    d      (평균-최빈)/sd   sqrt(2/d)      왜도/2    P(Q<=d)
    3        0.816497    0.816497  0.816497   0.608375
    5        0.632456    0.632456  0.632456   0.584120
   10        0.447214    0.447214  0.447214   0.559507
   20        0.316228    0.316228  0.316228   0.542070
   50        0.200000    0.200000  0.200000   0.526602
  100        0.141421    0.141421  0.141421   0.518808
  500        0.063246    0.063246  0.063246   0.508411

d = 5:  f(최빈 3) = 0.154180,  f(중앙 4.3515) = 0.137036,  f(평균 5) = 0.122042
        P(Q <= 5) = 0.584120   (quad 로 재확인 0.584120)

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

k = 5                      # 자유도. 표준정규 k개를 제곱해 더한 것의 분포다.
chi2 = stats.chi2(df=k)

# x 범위를 분위수로 정한다. 눈대중으로 (0, 20) 같은 범위를 쓰면
# 자유도가 바뀔 때마다 그림이 잘리거나 남는다.
# ppf(1e-6)부터 ppf(1-1e-6)까지 잡으면 어떤 k에서도 꼬리까지 알맞게 담긴다.
x = np.linspace(chi2.ppf(1e-6), chi2.ppf(1 - 1e-6), 600)
y = chi2.pdf(x)

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, y, lw=2, label=f"χ² PDF (k={k})")
# 평균은 정확히 k, 최빈값은 k-2 (k >= 2일 때).
# 둘이 다르다는 것이 곧 이 분포가 오른쪽으로 치우쳐 있다는 뜻이다.
ax.axvline(k, linestyle='--', alpha=0.8, label=f"mean = {k}")
ax.axvline(max(k - 2, 0), linestyle=':', alpha=0.8, label=f"mode = {max(k-2, 0)}")
ax.set_title("Chi-square Distribution — PDF")
ax.set_xlabel("x")
ax.set_ylabel("density")
ax.legend()
ax.grid(True, linestyle=":")
plt.tight_layout()
plt.show()

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

np.random.seed(42)
d = 5
# 정의를 그대로 실행한다. 표준정규 5개를 뽑아 제곱해 더하기를 5만 번.
#   rvs((d, 50000)) 이 (5, 50000) 배열을 주고
#   axis=0 으로 더하면 열마다(= 시행마다) 5개의 제곱합이 나온다.
z = stats.norm.rvs(size=(d, 50_000))
q = (z ** 2).sum(axis=0)

x = np.linspace(0.01, 25, 400)
fig, ax = plt.subplots(figsize=(12, 3))
ax.hist(q, bins=80, density=True, alpha=0.5, label='sum of 5 squared normals')
ax.plot(x, stats.chi2(d).pdf(x), 'r-', lw=2, label='chi2(5) pdf')
ax.set_xlabel('q')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

print(f"sample mean {q.mean():.3f} (theory {d})")
print(f"sample var  {q.var():.3f} (theory {2*d})")

sample mean 5.000 (theory 5)
sample var  10.022 (theory 10)

import numpy as np
from scipy import stats

np.random.seed(42)
d, N = 5, 50_000
q = (stats.norm.rvs(size=(d, N)) ** 2).sum(axis=0)

# 누율은 kappa_n = 2^(n-1) (n-1)! d 다(연습문제 9). 쓸 것은 kappa_2, kappa_4.
k2, k4 = 2 * d, 48 * d

# 표본평균과 표본분산의 이론 표준오차. 이것이 "얼마나 맞아야 맞는 것인가"의 기준이다.
se_mean = np.sqrt(k2 / N)
se_var = np.sqrt(k4 / N + 2 * k2 ** 2 / (N - 1))

print(f"표본평균 {q.mean():.4f}   이론 {d}   SE = sqrt(2d/N) = {se_mean:.4f}"
      f"   편차/SE = {abs(q.mean() - d) / se_mean:.2f}")
print(f"표본분산 {q.var():.4f}  이론 {k2}  SE = {se_var:.4f}"
      f"   편차/SE = {abs(q.var() - k2) / se_var:.2f}")

ks = stats.kstest(q, "chi2", args=(d,))
print(f"\nKS 거리 = {ks.statistic:.5f}   (5% 임계 ~ 1.36/sqrt(N) = {1.36 / np.sqrt(N):.5f}),"
      f"  p = {ks.pvalue:.4f}")

# MGF 를 직접 재 본다. 이론은 (1-2t)^(-d/2) 이고 t < 1/2 에서만 정의된다.
# 그런데 추정량 mean(exp(tQ)) 의 **분산**은 E[exp(2tQ)] 를 요구하므로
# t < 1/4 에서만 유한하다. 그 경계를 넘으면 수가 맞지 않기 시작한다.
print(f"\n{'t':>7}{'표본 MGF':>13}{'이론 MGF':>13}{'SE':>12}{'편차/SE':>10}")
for t in (-0.50, -0.10, 0.10, 0.20, 0.24, 0.25, 0.30, 0.40):
    theory = (1 - 2 * t) ** (-d / 2)
    sample = np.mean(np.exp(t * q))
    if 2 * t < 0.5:
        se = np.sqrt(((1 - 4 * t) ** (-d / 2) - theory ** 2) / N)
        print(f"{t:>7.2f}{sample:>13.5f}{theory:>13.5f}{se:>12.5f}"
              f"{abs(sample - theory) / se:>10.2f}")
    else:
        print(f"{t:>7.2f}{sample:>13.5f}{theory:>13.5f}{'무한':>11}{'--':>10}")

표본평균 4.9999   이론 5   SE = sqrt(2d/N) = 0.0141   편차/SE = 0.01
표본분산 10.0223  이론 10  SE = 0.0938   편차/SE = 0.24

KS 거리 = 0.00234   (5% 임계 ~ 1.36/sqrt(N) = 0.00608),  p = 0.9473

      t       표본 MGF       이론 MGF          SE     편차/SE
  -0.50      0.17678      0.17678     0.00081      0.00
  -0.10      0.63398      0.63394     0.00077      0.05
   0.10      1.74732      1.74693     0.00327      0.12
   0.20      3.59425      3.58610     0.02934      0.28
   0.24      5.15315      5.12852     0.24895      0.10
   0.25      5.68914      5.65685         무한        --
   0.30      9.99557      9.88212         무한        --
   0.40     51.29793     55.90170         무한        --

import numpy as np
from scipy import stats

x = np.array([0.5, 1.0, 2.0, 4.0])

# chi2(d) = Gamma(shape=d/2, scale=2). scipy의 gamma는 a=shape, scale=scale 이다.
print("chi2(5) :", np.round(stats.chi2(5).pdf(x), 6))
print("gamma   :", np.round(stats.gamma(a=2.5, scale=2).pdf(x), 6))

# d=2 이면 비율 1/2 인 지수분포다. scipy의 expon은 scale=1/rate 를 받는다.
print("chi2(2) :", np.round(stats.chi2(2).pdf(x), 6))
print("expon   :", np.round(stats.expon(scale=2).pdf(x), 6))

chi2(5) : [0.036616 0.080657 0.138369 0.143976]
gamma   : [0.036616 0.080657 0.138369 0.143976]
chi2(2) : [0.3894   0.303265 0.18394  0.067668]
expon   : [0.3894   0.303265 0.18394  0.067668]

import numpy as np
from scipy import stats

x = np.array([0.5, 1.0, 2.0, 4.0])

# 밀도만이 아니라 CDF 와 분위수까지 같아야 "같은 분포"다.
for d in (1, 2, 3, 5, 8):
    g = stats.gamma(a=d / 2, scale=2)
    c = stats.chi2(d)
    print(f"d={d}:  gamma(a={d/2}, scale=2) 와 pdf {np.allclose(c.pdf(x), g.pdf(x))},"
          f"  cdf {np.allclose(c.cdf(x), g.cdf(x))},"
          f"  ppf {np.allclose(c.ppf([.05, .5, .95]), g.ppf([.05, .5, .95]))}")

# 함정: scipy 의 gamma·expon 은 **척도**를 받는다. 비율을 넣으면 조용히 틀린다.
print(f"\n척도 2 를 바르게 준 경우   : 평균 {stats.gamma(a=2.5, scale=2).mean():.4f}"
      f"  (chi2(5) 의 평균 {stats.chi2(5).mean():.4f})")
print(f"비율 1/2 를 잘못 준 경우   : 평균 {stats.gamma(a=2.5, scale=0.5).mean():.4f}"
      f"  <- 4배 작다. 예외도 경고도 없다")
print(f"expon(scale=2)  평균 {stats.expon(scale=2).mean():.4f}"
      f"   expon(scale=0.5) 평균 {stats.expon(scale=0.5).mean():.4f}")

# 짝수 자유도는 얼랑분포, 곧 독립 Exp(1/2) 를 d/2 개 더한 것이다.
rng = np.random.default_rng(7)
print(f"\n{'k':>3}{'chi2_(2k)':>12}{'KS 거리':>12}{'p값':>9}{'표본평균':>12}{'이론':>7}")
for k in (1, 2, 3, 5):
    s = rng.exponential(scale=2, size=(200_000, k)).sum(axis=1)
    r = stats.kstest(s, "chi2", args=(2 * k,))
    print(f"{k:>3}{f'chi2_{2*k}':>12}{r.statistic:>12.5f}{r.pvalue:>9.4f}"
          f"{s.mean():>12.4f}{2*k:>7}")
print(f"KS 의 5% 임계 ~ 1.36/sqrt(200000) = {1.36 / np.sqrt(200_000):.5f}")

d=1:  gamma(a=0.5, scale=2) 와 pdf True,  cdf True,  ppf True
d=2:  gamma(a=1.0, scale=2) 와 pdf True,  cdf True,  ppf True
d=3:  gamma(a=1.5, scale=2) 와 pdf True,  cdf True,  ppf True
d=5:  gamma(a=2.5, scale=2) 와 pdf True,  cdf True,  ppf True
d=8:  gamma(a=4.0, scale=2) 와 pdf True,  cdf True,  ppf True

척도 2 를 바르게 준 경우   : 평균 5.0000  (chi2(5) 의 평균 5.0000)
비율 1/2 를 잘못 준 경우   : 평균 1.2500  <- 4배 작다. 예외도 경고도 없다
expon(scale=2)  평균 2.0000   expon(scale=0.5) 평균 0.5000

  k   chi2_(2k)       KS 거리       p값        표본평균     이론
  1      chi2_2     0.00230   0.2421      1.9993      2
  2      chi2_4     0.00148   0.7761      3.9948      4
  3      chi2_6     0.00157   0.7083      6.0015      6
  5     chi2_10     0.00324   0.0296      9.9865     10
KS 의 5% 임계 ~ 1.36/sqrt(200000) = 0.00304

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

z = np.linspace(-4, 4, 400)

fig, ax = plt.subplots(figsize=(12, 3))
# (Q - d)/sqrt(2d) 를 그린다. 중심극한정리에 따라 N(0,1)로 가야 한다.
# 왜도가 sqrt(8/d) 이므로 d=5에서 1.26, d=20에서 0.63, d=100에서 0.28 이다.
# 두 가지를 보라.
#   Q >= 0 이므로 왼쪽은 -sqrt(d/2) 에서 **뚝 잘린다**(d=5면 -1.58).
#   대신 오른쪽 꼬리가 정규보다 두껍다. 치우침이 사라지는 속도가 느리다.
for d in [5, 20, 100]:
    ax.plot(z, stats.chi2(d).pdf(d + z * np.sqrt(2 * d)) * np.sqrt(2 * d),
            lw=2, label=f'standardized chi2({d})')
ax.plot(z, stats.norm.pdf(z), 'k--', lw=2, label='N(0, 1)')
ax.set_xlabel('(q - d) / sqrt(2d)')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

import numpy as np
from scipy import stats

# 표준화한 밀도. Y = (Q - d)/sqrt(2d) 이므로 f_Y(z) = sqrt(2d) * f_Q(d + z*sqrt(2d)).
def std_pdf(z, d):
    return stats.chi2(d).pdf(d + z * np.sqrt(2 * d)) * np.sqrt(2 * d)

print(f"{'d':>7}{'왼쪽끝':>10}{'최빈(유도)':>12}{'최빈(격자)':>12}"
      f"{'봉우리':>9}{'왜도':>9}{'초과첨도':>11}")
for d in (5, 20, 100, 400, 1600):
    left, mode = -np.sqrt(d / 2), -np.sqrt(2 / d)
    zs = np.linspace(left + 1e-9, 8, 400_001)
    pk = std_pdf(zs, d)
    g1, g2 = stats.chi2.stats(d, moments="sk")
    print(f"{d:>7}{left:>10.4f}{mode:>12.4f}{zs[pk.argmax()]:>12.4f}"
          f"{pk.max():>9.4f}{float(g1):>9.4f}{float(g2):>11.4f}")
print(f"{'N(0,1)':>7}{'-inf':>10}{0.0:>12.4f}{0.0:>12.4f}"
      f"{1 / np.sqrt(2 * np.pi):>9.4f}{0.0:>9.4f}{0.0:>11.4f}")

# 정규분포까지의 KS 거리. 유도에 따르면 1/(3 sqrt(pi d)) 이고
# 가장 많이 벌어지는 자리는 꼬리가 아니라 z = 0, 곧 평균이다.
print(f"\n{'d':>7}{'KS 거리':>11}{'1/(3 sqrt(pi d))':>19}{'비':>8}{'최대 벌어짐 z':>15}")
for d in (5, 20, 100, 400, 1600, 6400):
    zz = np.linspace(-np.sqrt(d / 2), 10, 400_001)
    gap = np.abs(stats.chi2(d).cdf(d + zz * np.sqrt(2 * d)) - stats.norm.cdf(zz))
    ks, pred = gap.max(), 1 / (3 * np.sqrt(np.pi * d))
    print(f"{d:>7}{ks:>11.6f}{pred:>19.6f}{ks / pred:>8.4f}{zz[gap.argmax()]:>15.4f}")

# 그림 창은 [-4, 4] 다. 왼쪽 절단이 그 안에 들어오는 자유도만 잘린 모습이 보인다.
print("\n그림 창 [-4, 4] 안에 왼쪽 끝 -sqrt(d/2) 가 들어오는가:")
for d in (5, 20, 100):
    left = -np.sqrt(d / 2)
    print(f"  d={d:>4}: {left:>8.4f}  ->  "
          f"{'창 안이라 절단이 보인다' if left > -4 else '창 밖이라 절단이 안 보인다'}")

      d       왼쪽끝      최빈(유도)      최빈(격자)      봉우리       왜도       초과첨도
      5   -1.5811     -0.6325     -0.6325   0.4876   1.2649     2.4000
     20   -3.1623     -0.3162     -0.3162   0.4166   0.6325     0.6000
    100   -7.0711     -0.1414     -0.1414   0.4023   0.2828     0.1200
    400  -14.1421     -0.0707     -0.0707   0.3998   0.1414     0.0300
   1600  -28.2843     -0.0354     -0.0353   0.3992   0.0707     0.0075
 N(0,1)      -inf      0.0000      0.0000   0.3989   0.0000     0.0000

      d      KS 거리   1/(3 sqrt(pi d))       비       최대 벌어짐 z
      5   0.084459           0.084104  1.0042        -0.0516
     20   0.042114           0.042052  1.0015        -0.0262
    100   0.018812           0.018806  1.0003        -0.0118
    400   0.009404           0.009403  1.0001        -0.0059
   1600   0.004702           0.004702  1.0000        -0.0029
   6400   0.002351           0.002351  1.0000        -0.0014

그림 창 [-4, 4] 안에 왼쪽 끝 -sqrt(d/2) 가 들어오는가:
  d=   5:  -1.5811  ->  창 안이라 절단이 보인다
  d=  20:  -3.1623  ->  창 안이라 절단이 보인다
  d= 100:  -7.0711  ->  창 밖이라 절단이 안 보인다

import numpy as np
from scipy import stats

print(f"{'d':>5}{'sqrt(8/d)':>13}{'scipy 왜도':>13}{'12/d':>10}{'scipy 초과첨도':>16}")
for d in (1, 2, 5, 10, 50, 100):
    s, kk = stats.chi2.stats(d, moments="sk")
    print(f"{d:>5}{np.sqrt(8 / d):>13.4f}{float(s):>13.4f}"
          f"{12 / d:>10.4f}{float(kk):>16.4f}")

    d    sqrt(8/d)     scipy 왜도      12/d      scipy 초과첨도
    1       2.8284       2.8284   12.0000         12.0000
    2       2.0000       2.0000    6.0000          6.0000
    5       1.2649       1.2649    2.4000          2.4000
   10       0.8944       0.8944    1.2000          1.2000
   50       0.4000       0.4000    0.2400          0.2400
  100       0.2828       0.2828    0.1200          0.1200

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = np.linspace(0.1, 20, 7)
print("chi2(d=5) 와 Gamma(a=2.5, scale=2) 의 밀도가 같은가:",
      np.allclose(stats.chi2.pdf(x, 5), stats.gamma.pdf(x, 2.5, scale=2)))

n, mu, sig = 8, 3.0, 2.0
X = rng.normal(mu, sig, (400_000, n))
m, s2 = X.mean(1), X.var(1, ddof=1)
Q = (n - 1) * s2 / sig ** 2

print(f"\n(n-1)s^2/sigma^2:  평균 {Q.mean():.4f} (df={n-1}),"
      f"  분산 {Q.var():.4f} (2df={2*(n-1)})")
print(f"  chi2_{n-1} 까지의 KS 거리 = "
      f"{stats.kstest(Q, 'chi2', args=(n - 1,)).statistic:.4f}")
print(f"  corr(Xbar, s^2) = {np.corrcoef(m, s2)[0, 1]:+.5f}   <- 독립")

left = ((X - mu) ** 2).sum(1) / sig ** 2
right1 = n * (m - mu) ** 2 / sig ** 2
print(f"\n분해:  좌변 평균 {left.mean():.4f} (df={n})"
      f"  =  {right1.mean():.4f} (df=1)  +  {Q.mean():.4f} (df={n-1})")

chi2(d=5) 와 Gamma(a=2.5, scale=2) 의 밀도가 같은가: True

(n-1)s^2/sigma^2:  평균 6.9973 (df=7),  분산 14.0410 (2df=14)
  chi2_7 까지의 KS 거리 = 0.0017
  corr(Xbar, s^2) = +0.00095   <- 독립

분해:  좌변 평균 7.9963 (df=8)  =  0.9990 (df=1)  +  6.9973 (df=7)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
for m, n in [(4, 50), (4, 500), (10, 500)]:
    p = np.full(m, 1 / m)
    O = rng.multinomial(n, p, size=200_000)
    X2 = ((O - n * p) ** 2 / (n * p)).sum(1)
    print(f"범주 {m:>2}, n={n:>4}:  평균 {X2.mean():>7.4f} (df={m-1}),"
          f"  분산 {X2.var():>7.4f} (2df={2*(m-1)}),"
          f"  KS = {stats.kstest(X2, 'chi2', args=(m - 1,)).statistic:.4f}")

범주  4, n=  50:  평균  3.0022 (df=3),  분산  5.9008 (2df=6),  KS = 0.0507
범주  4, n= 500:  평균  3.0051 (df=3),  분산  5.9967 (2df=6),  KS = 0.0084
범주 10, n= 500:  평균  8.9948 (df=9),  분산 17.9558 (2df=18),  KS = 0.0031

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
p1 = np.array([0.30, 0.25, 0.25, 0.20])
p0 = np.full(4, 0.25)
crit = stats.chi2.ppf(0.95, 3)
print(f"기각역: X^2 > {crit:.4f}")
print(f"{'n':>7}{'lambda':>12}{'이론 검정력':>15}{'모의 검정력':>15}")
for n in (100, 300, 1000, 3000):
    lam = n * ((p1 - p0) ** 2 / p0).sum()
    O = rng.multinomial(n, p1, size=100_000)
    X2 = ((O - n * p0) ** 2 / (n * p0)).sum(1)
    print(f"{n:>7}{lam:>12.4f}{stats.ncx2.sf(crit, 3, lam):>15.4f}"
          f"{np.mean(X2 > crit):>15.4f}")

기각역: X^2 > 7.8147
      n      lambda         이론 검정력         모의 검정력
    100      2.0000         0.1922         0.1903
    300      6.0000         0.5181         0.5168
   1000     20.0000         0.9751         0.9765
   3000     60.0000         1.0000         1.0000