F 분포¶

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

x = np.linspace(0.01, 4, 400)

fig, ax = plt.subplots(figsize=(12, 3))
# 분자 자유도 d1이 0 근처의 모양을, 분모 자유도 d2가 꼬리의 무게를 정한다.
#   (1, 10)  : d1=1 이라 0에서 발산한다. 이것이 t(10)^2 의 분포다.
#   (5, 20)  : 봉우리가 0.545, 평균이 1.111. 오른쪽으로 치우쳐 있다.
#   (20, 20) : 두 자유도가 같아 1 근처에 모인다.
#   (50, 50) : 더 좁아진다. 자유도가 크면 비가 1에서 벗어나기 어렵다.
for d1, d2 in [(1, 10), (5, 20), (20, 20), (50, 50)]:
    ax.plot(x, stats.f(d1, d2).pdf(x), lw=2, label=f'F({d1}, {d2})')
ax.axvline(1, color='gray', ls=':', lw=1)   # 비가 1이면 두 분산이 같다는 뜻
ax.set_ylim(0, 2.2)
ax.set_xlabel('x')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

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


def log_const(d1, d2):
    """밀도 앞의 베타 상수의 로그. gammaln 을 쓰지 않으면 자유도가 조금만 커도 넘친다."""
    return (special.gammaln((d1 + d2) / 2)
            - special.gammaln(d1 / 2) - special.gammaln(d2 / 2))


def head_const(d1, d2):     # x -> 0+ 에서  f(x) ~ C x^(d1/2 - 1)
    return np.exp(log_const(d1, d2) + (d1 / 2) * np.log(d1 / d2))


def tail_const(d1, d2):     # x -> inf 에서 f(x) ~ A x^-(d2/2 + 1)
    return np.exp(log_const(d1, d2) + (d2 / 2) * np.log(d2 / d1))


# 그림에 그린 네 곡선에 대해 두 근사식을 양 끝에서 실제 밀도와 견준다.
print(f"{'자유도':>9}{'d1/2-1':>8}{'C':>12}{'비(x=1e-8)':>12}"
      f"{'-(d2/2+1)':>11}{'A':>12}{'비(x=1e6)':>11}")
for d1, d2 in [(1, 10), (2, 10), (5, 20), (20, 20), (50, 50)]:
    f = stats.f(d1, d2)
    C, A = head_const(d1, d2), tail_const(d1, d2)
    r0 = f.pdf(1e-8) / (C * (1e-8) ** (d1 / 2 - 1))
    r1 = f.pdf(1e6) / (A * (1e6) ** (-(d2 / 2 + 1)))
    print(f"({d1:2d},{d2:3d}){d1/2-1:>8.1f}{C:>12.4g}{r0:>12.6f}"
          f"{-(d2/2+1):>11.1f}{A:>12.4g}{r1:>11.6f}")

# d1 = 2 는 경계다. 지수가 0 이 되고, 상수가 d2 와 무관하게 정확히 1 이 된다.
print("\nd1 = 2 에서 f(0+) = 1 (d2 와 무관):")
for d2 in (2, 10, 50, 1000):
    print(f"  d2 = {d2:5d}:  f(1e-12) = {stats.f(2, d2).pdf(1e-12):.12f}")


def moment(k, d1, d2):
    """E[X^k] 의 닫힌 꼴. 2k >= d2 면 Gamma(d2/2 - k) 의 인수가 양수가 아니다."""
    if 2 * k >= d2:
        return np.inf
    return (d2 / d1) ** k * np.exp(
        special.gammaln(d1 / 2 + k) + special.gammaln(d2 / 2 - k)
        - special.gammaln(d1 / 2) - special.gammaln(d2 / 2))


# quad 는 [0, inf] 를 한 번에 주면 0 근처의 질량을 놓친다. 1 에서 끊어 두 조각으로 쟌다.
d1, d2 = 5, 20
print(f"\nF({d1},{d2}) 의 적률 — 닫힌 꼴과 수치적분:")
for k in (1, 2, 5, 9):
    q = (integrate.quad(lambda x: x ** k * stats.f(d1, d2).pdf(x), 0, 1)[0]
         + integrate.quad(lambda x: x ** k * stats.f(d1, d2).pdf(x), 1, np.inf)[0])
    print(f"  k = {k}:  닫힌꼴 {moment(k, d1, d2):13.6f}   quad {q:13.6f}")

# 2k >= d2 면 꼬리 적분이 멈추지 않는다. T 를 키우며 본다.
print(f"\n꼬리 적분 ∫_1^T x^k f(x) dx   (d2/2 = {d2 // 2} 이 경계):")
for k in (9, 10, 11):
    vals = [integrate.quad(lambda x: x ** k * stats.f(d1, d2).pdf(x), 1, T)[0]
            for T in (1e2, 1e4, 1e6)]
    tag = "수렴" if 2 * k < d2 else "발산"
    print(f"  k = {k:2d} (2k = {2 * k:2d}, {tag}): "
          + "  ".join(f"T=1e{e}  {v:10.4g}" for e, v in zip((2, 4, 6), vals)))

# 그림의 테두리가 무엇을 자르는지 센다. 그림은 x in [0.01, 4], y in [0, 2.2] 다.
print("\n그림의 테두리 밖에 있는 것:")
for d1, d2 in [(1, 10), (5, 20), (20, 20), (50, 50)]:
    f = stats.f(d1, d2)
    print(f"  F({d1:2d},{d2:3d}): f(0.01) = {f.pdf(0.01):9.4f}   f(4) = {f.pdf(4):.3e}"
          f"   P(X > 4) = {f.sf(4):.4f}")
grid = np.linspace(1e-6, 4, 4_000_001)
over = grid[stats.f(1, 10).pdf(grid) > 2.2]
print(f"  F(1,10) 의 밀도가 2.2(세로축 상한)를 넘는 구간은 x < {over.max():.4f} 뿐이다")

      자유도  d1/2-1           C   비(x=1e-8)  -(d2/2+1)           A   비(x=1e6)
( 1, 10)    -0.5      0.3891    1.000000       -6.0    1.23e+05   0.999945
( 2, 10)     0.0           1    1.000000       -6.0   1.562e+04   0.999970
( 5, 20)     1.5       8.865    1.000000      -11.0   2.975e+08   0.999950
(20, 20)     9.0   9.238e+05    1.000000      -11.0   9.238e+05   0.999980
(50, 50)    24.0    1.58e+15    1.000000      -26.0    1.58e+15   0.999950

d1 = 2 에서 f(0+) = 1 (d2 와 무관):
  d2 =     2:  f(1e-12) = 0.999999999998
  d2 =    10:  f(1e-12) = 0.999999999999
  d2 =    50:  f(1e-12) = 0.999999999999
  d2 =  1000:  f(1e-12) = 0.999999999999

F(5,20) 의 적률 — 닫힌 꼴과 수치적분:
  k = 1:  닫힌꼴      1.111111   quad      1.111111
  k = 2:  닫힌꼴      1.944444   quad      1.944444
  k = 5:  닫힌꼴     95.333333   quad     95.333333
  k = 9:  닫힌꼴 6466460.000000   quad 6466460.000000

꼬리 적분 ∫_1^T x^k f(x) dx   (d2/2 = 10 이 경계):
  k =  9 (2k = 18, 수렴): T=1e2   4.119e+06  T=1e4   6.437e+06  T=1e6   6.466e+06
  k = 10 (2k = 20, 발산): T=1e2   1.775e+08  T=1e4   1.418e+09  T=1e6   2.786e+09
  k = 11 (2k = 22, 발산): T=1e2   1.007e+10  T=1e4    2.89e+12  T=1e6   2.973e+14

그림의 테두리 밖에 있는 것:
  F( 1, 10): f(0.01) =    3.8698   f(4) = 3.057e-02   P(X > 4) = 0.0734
  F( 5, 20): f(0.01) =    0.0086   f(4) = 1.224e-02   P(X > 4) = 0.0112
  F(20, 20): f(0.01) =    0.0000   f(4) = 2.539e-03   P(X > 4) = 0.0016
  F(50, 50): f(0.01) =    0.0000   f(4) = 5.008e-06   P(X > 4) = 0.0000
  F(1,10) 의 밀도가 2.2(세로축 상한)를 넘는 구간은 x < 0.0303 뿐이다

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

# F 분포는 자유도가 둘이다. 두 카이제곱을 각자의 자유도로 나눈 뒤의 비율이다.
#   dfn = 분자 자유도, dfd = 분모 자유도
d1, d2 = 5, 12
f_dist = stats.f(dfn=d1, dfd=d2)

x = np.linspace(f_dist.ppf(1e-6), f_dist.ppf(1 - 1e-6), 600)
y = f_dist.pdf(x)

# 평균은 분모 자유도만으로 정해지며 d2 > 2 일 때만 존재한다.
# d2가 작으면 꼬리가 매우 무거워 평균이 아예 없다.
mean = d2 / (d2 - 2)
mode = ((d1 - 2) / d1) * (d2 / (d2 + 2))
# 최빈값 0.514 < 1 < 평균 1.2 다. 둘이 1을 사이에 두고 갈라져 있다.

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, y, lw=2, label=f"F PDF (d1={d1}, d2={d2})")
ax.axvline(mean, linestyle='--', alpha=0.85, label=f"mean = {mean:.3f}")
ax.axvline(mode, linestyle=':', alpha=0.85, label=f"mode = {mode:.3f}")
ax.set_title("F Distribution — PDF")
ax.set_xlabel("x")
ax.set_ylabel("density")
ax.legend()
ax.grid(True, linestyle=":")
plt.tight_layout()
plt.show()

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

d1, d2 = 5, 12
F = stats.f(d1, d2)

# (1) 에서 유도한 닫힌 꼴들
mode_cf = d2 * (d1 - 2) / (d1 * (d2 + 2))
mean_cf = d2 / (d2 - 2)
var_cf = 2 * d2 ** 2 * (d1 + d2 - 2) / (d1 * (d2 - 2) ** 2 * (d2 - 4))

# 최빈값: 격자와 수치최적화 둘로 확인한다.
grid = np.linspace(1e-9, 6, 6_000_001)
mode_grid = grid[F.pdf(grid).argmax()]
mode_opt = optimize.minimize_scalar(lambda x: -F.logpdf(x), bounds=(1e-9, 6),
                                    method="bounded", options={"xatol": 1e-12}).x

# 평균·분산: quad 로 확인한다. 0 근처에 질량이 몰려 있으므로 1 에서 끊는다.
m1 = (integrate.quad(lambda x: x * F.pdf(x), 0, 1)[0]
      + integrate.quad(lambda x: x * F.pdf(x), 1, np.inf)[0])
m2 = (integrate.quad(lambda x: x * x * F.pdf(x), 0, 1)[0]
      + integrate.quad(lambda x: x * x * F.pdf(x), 1, np.inf)[0])

print(f"F({d1},{d2}) — 닫힌 꼴과 수치값")
print(f"  최빈값  닫힌꼴 {mode_cf:.9f}   격자 {mode_grid:.9f}   최적화 {mode_opt:.9f}")
print(f"  평균    닫힌꼴 {mean_cf:.9f}   quad {m1:.9f}")
print(f"  분산    닫힌꼴 {var_cf:.9f}   quad {m2 - m1 ** 2:.9f}")
print(f"  E[1/V] (V~chi2_{d2}):  닫힌꼴 1/(d2-2) = {1 / (d2 - 2):.9f}   quad "
      f"{integrate.quad(lambda v: stats.chi2(d2).pdf(v) / v, 0, np.inf)[0]:.9f}")

print(f"\n네 위치의 순서  (표준편차 {np.sqrt(var_cf):.4f})")
print(f"  최빈값 {mode_cf:.4f}  <  중앙값 {F.median():.4f}  <  1  <  평균 {mean_cf:.4f}")
print(f"  순서가 맞는가: {mode_cf < F.median() < 1 < mean_cf}")

print("\nd1 <= 2 에서는 최빈값 공식이 내부의 봉우리를 가리키지 않는다")
for a in (1, 2, 3, 5):
    print(f"  d1 = {a}:  공식값 {d2 * (a - 2) / (a * (d2 + 2)):+.6f}"
          f"   (d1 > 2 인가: {a > 2})")

print("\n두 모분산이 정말 같을 때 관측된 비가 1 아래로 떨어질 확률")
# 베타 경유로도 같은 값이 나오는지 함께 확인한다.
for a, b in [(5, 12), (12, 5), (5, 20), (20, 5), (10, 10), (3, 3)]:
    p = stats.f(a, b).cdf(1)
    p_beta = stats.beta(a / 2, b / 2).cdf(a / (a + b))
    print(f"  F({a:2d},{b:2d}): P(X < 1) = {p:.6f}   베타 경유 {p_beta:.6f}"
          f"   중앙값 {stats.f(a, b).median():.6f}")

print("\n역수 관계가 주는 항등식  P(F_{d1,d2} < 1) = 1 - P(F_{d2,d1} < 1)")
for a, b in [(5, 12), (5, 20), (3, 34)]:
    lhs = stats.f(a, b).cdf(1)
    rhs = 1 - stats.f(b, a).cdf(1)
    print(f"  ({a},{b}): {lhs:.12f} vs {rhs:.12f}   차 {lhs - rhs:+.2e}")

print("\nd1 = d2 이면 중앙값이 정확히 1 이다")
for d in (3, 5, 12, 40):
    print(f"  F({d},{d}).median() = {stats.f(d, d).median():.12f}   "
          f"P(X < 1) = {stats.f(d, d).cdf(1):.12f}")

print(f"\n평균보다 작을 확률:  P(X < {mean_cf:.4f}) = {F.cdf(mean_cf):.6f}")
print("자유도를 함께 키우면 최빈값과 평균이 둘 다 1 로 모인다")
print(f"{'d1=d2':>8}{'최빈값':>12}{'평균':>10}{'평균-최빈값':>13}")
for d in (5, 10, 20, 50, 200, 1000):
    mo = d * (d - 2) / (d * (d + 2))
    me = d / (d - 2)
    print(f"{d:>8}{mo:>12.6f}{me:>10.6f}{me - mo:>13.6f}")

F(5,12) — 닫힌 꼴과 수치값
  최빈값  닫힌꼴 0.514285714   격자 0.514286001   최적화 0.514285678
  평균    닫힌꼴 1.200000000   quad 1.200000000
  분산    닫힌꼴 1.080000000   quad 1.080000000
  E[1/V] (V~chi2_12):  닫힌꼴 1/(d2-2) = 0.100000000   quad 0.100000000

네 위치의 순서  (표준편차 1.0392)
  최빈값 0.5143  <  중앙값 0.9212  <  1  <  평균 1.2000
  순서가 맞는가: True

d1 <= 2 에서는 최빈값 공식이 내부의 봉우리를 가리키지 않는다
  d1 = 1:  공식값 -0.857143   (d1 > 2 인가: False)
  d1 = 2:  공식값 +0.000000   (d1 > 2 인가: False)
  d1 = 3:  공식값 +0.285714   (d1 > 2 인가: True)
  d1 = 5:  공식값 +0.514286   (d1 > 2 인가: True)

두 모분산이 정말 같을 때 관측된 비가 1 아래로 떨어질 확률
  F( 5,12): P(X < 1) = 0.541803   베타 경유 0.541803   중앙값 0.921242
  F(12, 5): P(X < 1) = 0.458197   베타 경유 0.458197   중앙값 1.085492
  F( 5,20): P(X < 1) = 0.556975   베타 경유 0.556975   중앙값 0.900376
  F(20, 5): P(X < 1) = 0.443025   베타 경유 0.443025   중앙값 1.110647
  F(10,10): P(X < 1) = 0.500000   베타 경유 0.500000   중앙값 1.000000
  F( 3, 3): P(X < 1) = 0.500000   베타 경유 0.500000   중앙값 1.000000

역수 관계가 주는 항등식  P(F_{d1,d2} < 1) = 1 - P(F_{d2,d1} < 1)
  (5,12): 0.541803300486 vs 0.541803300486   차 +1.11e-16
  (5,20): 0.556974815315 vs 0.556974815315   차 +2.22e-16
  (3,34): 0.595308357974 vs 0.595308357974   차 -8.88e-16

d1 = d2 이면 중앙값이 정확히 1 이다
  F(3,3).median() = 1.000000000000   P(X < 1) = 0.500000000000
  F(5,5).median() = 1.000000000000   P(X < 1) = 0.500000000000
  F(12,12).median() = 1.000000000000   P(X < 1) = 0.500000000000
  F(40,40).median() = 1.000000000000   P(X < 1) = 0.500000000000

평균보다 작을 확률:  P(X < 1.2000) = 0.633977
자유도를 함께 키우면 최빈값과 평균이 둘 다 1 로 모인다
   d1=d2         최빈값        평균       평균-최빈값
       5    0.428571  1.666667     1.238095
      10    0.666667  1.250000     0.583333
      20    0.818182  1.111111     0.292929
      50    0.923077  1.041667     0.118590
     200    0.980198  1.010101     0.029903
    1000    0.996008  1.002004     0.005996

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

np.random.seed(42)
d1, d2 = 5, 20
# 정의를 그대로 실행한다. 두 카이제곱을 **따로** 뽑아야 독립이 보장된다.
# 각각을 자기 자유도로 나누는 것이 핵심이다. 그래야 비가 1 근처에 놓인다.
u = stats.chi2(d1).rvs(200_000)
v = stats.chi2(d2).rvs(200_000)
x = (u / d1) / (v / d2)

grid = np.linspace(0.01, 5, 400)
fig, ax = plt.subplots(figsize=(12, 3))
# 구간을 [0, 5]로 고정한다. 오른쪽 꼬리가 두꺼워 x > 20 인 값도 나온다.
ax.hist(x, bins=200, range=(0, 5), density=True, alpha=0.5, label='(U/d1) / (V/d2)')
ax.plot(grid, stats.f(d1, d2).pdf(grid), 'r-', lw=2, label='F(5, 20) pdf')
ax.set_xlabel('x')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

print(f"sample mean {x.mean():.4f} (theory {d2/(d2-2):.4f})")
print(f"P(X > 2.711) = {np.mean(x > 2.711):.4f} (theory 0.05)")

sample mean 1.1128 (theory 1.1111)
P(X > 2.711) = 0.0507 (theory 0.05)

import numpy as np
from scipy import special, stats

d1, d2 = 5, 20

# ── (1) 베타에서 변수변환으로 얻은 밀도가 정리 1 의 밀도와 같은가 ──
a, b, c = d1 / 2, d2 / 2, d1 / d2


def pdf_from_beta(x):
    """Beta(a,b) 를 X = (b/a)Y/(1-Y) 로 보낸 밀도. 유도 결과를 그대로 적었다."""
    return c ** a / special.beta(a, b) * x ** (a - 1) * (1 + c * x) ** (-(a + b))


print("(1) 베타 경유로 유도한 밀도 vs scipy 의 F 밀도")
for x in (0.05, 0.3, 1.0, 2.711, 7.0):
    mine, ref = pdf_from_beta(x), stats.f(d1, d2).pdf(x)
    print(f"  x = {x:6.3f}:  유도 {mine:.12f}   scipy {ref:.12f}   차 {mine - ref:+.2e}")

# ── 쪽의 코드와 똑같은 씨앗·순서로 표본을 되살린다 ──
np.random.seed(42)
u = stats.chi2(d1).rvs(200_000)
v = stats.chi2(d2).rvs(200_000)
x = (u / d1) / (v / d2)
n = x.size

# ── (2) 같은 표본을 베타 쪽으로 보낸다 ──
y = u / (u + v)
print(f"\n(2) 단조변환 Y = d1*X/(d1*X+d2) 가 되는가: 최대오차 "
      f"{np.abs(y - d1 * x / (d1 * x + d2)).max():.3e}")
ks_f = stats.kstest(x, stats.f(d1, d2).cdf)
ks_b = stats.kstest(y, stats.beta(a, b).cdf)
print(f"  KS(X,  F({d1},{d2}))         D = {ks_f.statistic:.17f}  p = {ks_f.pvalue:.6f}")
print(f"  KS(Y,  Beta({a}, {b}))    D = {ks_b.statistic:.17f}  p = {ks_b.pvalue:.6f}")
print(f"  두 D 의 차 = {ks_b.statistic - ks_f.statistic:+.3e}"
      f"   상대차 = {abs(ks_b.statistic / ks_f.statistic - 1):.3e}"
      f"   배정도 간격 몇 개인가 = {round(abs(ks_b.statistic - ks_f.statistic) / np.spacing(ks_f.statistic))}")

# ── (3) 몬테카를로 오차 ──
mean_th = d2 / (d2 - 2)
var_th = 2 * d2 ** 2 * (d1 + d2 - 2) / (d1 * (d2 - 2) ** 2 * (d2 - 4))
se_mean = np.sqrt(var_th / n)
p_th = stats.f(d1, d2).sf(2.711)
se_p = np.sqrt(p_th * (1 - p_th) / n)
print("\n(3) 어긋난 몫이 몬테카를로 오차로 설명되는가")
print(f"  표본평균      {x.mean():.6f}   이론 {mean_th:.6f}   "
      f"차 {x.mean() - mean_th:+.6f} = {(x.mean() - mean_th) / se_mean:+.2f} SE  (SE {se_mean:.6f})")
print(f"  P(X > 2.711)  {np.mean(x > 2.711):.6f}   이론 {p_th:.6f}   "
      f"차 {np.mean(x > 2.711) - p_th:+.6f} = {(np.mean(x > 2.711) - p_th) / se_p:+.2f} SE  (SE {se_p:.6f})")
print(f"  KS 임계값(5%) ≈ 1.36/sqrt(n) = {1.36 / np.sqrt(n):.5f},  관측 D = {ks_f.statistic:.5f}")
print(f"  표본분산      {x.var(ddof=1):.6f}   이론 {var_th:.6f}")

# 그림의 히스토그램은 range=(0,5) 안에서만 정규화된다. 그래서 막대가 조금 높다.
h, edges = np.histogram(x, bins=200, range=(0, 5), density=True)
mid = (edges[:-1] + edges[1:]) / 2
sel = (mid > 0.2) & (mid < 3)
print("\n그림의 히스토그램은 구간 안에서만 정규화된다")
print(f"  P(X > 5): 모의 {np.mean(x > 5):.5f}   이론 {stats.f(d1, d2).sf(5):.5f}")
print(f"  막대/참밀도 비의 중앙값 (0.2 < x < 3) = {np.median(h[sel] / stats.f(d1, d2).pdf(mid[sel])):.6f}")
print(f"  예상되는 들뜸 1/P(X <= 5) = {1 / stats.f(d1, d2).cdf(5):.6f}")

(1) 베타 경유로 유도한 밀도 vs scipy 의 F 밀도
  x =  0.050:  유도 0.084857775634   scipy 0.084857775634   차 +4.16e-16
  x =  0.300:  유도 0.589862265493   scipy 0.589862265493   차 +3.00e-15
  x =  1.000:  유도 0.544878125182   scipy 0.544878125182   차 +2.55e-15
  x =  2.711:  유도 0.061415532019   scipy 0.061415532019   차 +2.64e-16
  x =  7.000:  유도 0.000529252506   scipy 0.000529252506   차 +6.29e-18

(2) 단조변환 Y = d1*X/(d1*X+d2) 가 되는가: 최대오차 2.220e-16
  KS(X,  F(5,20))         D = 0.00154584734914034  p = 0.725011
  KS(Y,  Beta(2.5, 10.0))    D = 0.00154584734914043  p = 0.725011
  두 D 의 차 = +8.327e-17   상대차 = 5.396e-14   배정도 간격 몇 개인가 = 384

(3) 어긋난 몫이 몬테카를로 오차로 설명되는가
  표본평균      1.112822   이론 1.111111   차 +0.001711 = +0.91 SE  (SE 0.001884)
  P(X > 2.711)  0.050735   이론 0.049993   차 +0.000742 = +1.52 SE  (SE 0.000487)
  KS 임계값(5%) ≈ 1.36/sqrt(n) = 0.00304,  관측 D = 0.00155
  표본분산      0.708868   이론 0.709877

그림의 히스토그램은 구간 안에서만 정규화된다
  P(X > 5): 모의 0.00378   이론 0.00393
  막대/참밀도 비의 중앙값 (0.2 < x < 3) = 1.003307
  예상되는 들뜸 1/P(X <= 5) = 1.003946

from scipy import stats

# (1) 역수 관계: F_0.05(5, 20) = 1 / F_0.95(20, 5)
print(f"F_0.05(5,20)     = {stats.f(5, 20).ppf(0.05):.6f}")
print(f"1 / F_0.95(20,5) = {1 / stats.f(20, 5).ppf(0.95):.6f}")

# (2) t^2 = F(1, d): 양측 97.5백분위점의 제곱이 F의 95백분위점이다.
print(f"t(10).ppf(0.975)^2 = {stats.t(10).ppf(0.975)**2:.6f}")
print(f"F(1,10).ppf(0.95)  = {stats.f(1, 10).ppf(0.95):.6f}")

# (3) d2가 크면 d1*F 가 카이제곱에 가까워진다.
print(f"4 * F(4,1000).ppf(0.95) = {4 * stats.f(4, 1000).ppf(0.95):.4f}")
print(f"chi2(4).ppf(0.95)       = {stats.chi2(4).ppf(0.95):.4f}")

F_0.05(5,20)     = 0.219388
1 / F_0.95(20,5) = 0.219388
t(10).ppf(0.975)^2 = 4.964603
F(1,10).ppf(0.95)  = 4.964603
4 * F(4,1000).ppf(0.95) = 9.5233
chi2(4).ppf(0.95)       = 9.4877

import numpy as np
from scipy import stats

print("(i) 역수 관계는 정확한 등식이다")
for al, (a, b) in zip((0.05, 0.10, 0.25, 0.01), [(5, 20), (12, 15), (3, 34), (4, 8)]):
    lhs = stats.f(a, b).ppf(al)
    rhs = 1 / stats.f(b, a).ppf(1 - al)
    print(f"  F_{al:<4}({a:2d},{b:2d}) = {lhs:.14f}   1/F_{1-al:.2f}({b:2d},{a:2d}) = {rhs:.14f}"
          f"   차 {lhs - rhs:+.1e}")

print("\n(ii) t^2 = F(1,d) 도 정확한 등식이다")
for d in (1, 2, 5, 10, 34, 200):
    lhs = stats.t(d).ppf(0.975) ** 2
    rhs = stats.f(1, d).ppf(0.95)
    print(f"  d = {d:3d}:  t_{{d,0.975}}^2 = {lhs:.12f}   F_{{0.95}}(1,d) = {rhs:.12f}"
          f"   상대차 {abs(lhs / rhs - 1):.1e}")
xs = np.array([0.1, 0.5, 1.0, 2.5, 4.9646, 12.0])
gap = np.max(np.abs(2 * stats.t(10).cdf(np.sqrt(xs)) - 1 - stats.f(1, 10).cdf(xs)))
print(f"  분포함수 항등식 2*T_10(sqrt(x)) - 1 = F(1,10).cdf(x) 의 최대오차: {gap:.1e}")

print("\n(iii) 카이제곱 극한만 근사다. 1차 보정항의 예측과 견준다")
d1, alpha = 4, 0.05
q = stats.chi2(d1).ppf(1 - alpha)
# delta = -(q^2/d2) * dlog f/dq = (q/d2)*(q/2 - d1/2 + 1)
coef = q * (q / 2 - d1 / 2 + 1)
print(f"  d1 = {d1}, q = chi2_{{0.95}}({d1}) = {q:.6f},  보정계수 q(q/2 - d1/2 + 1) = {coef:.4f}")
print(f"{'d2':>8}{'d1*F_0.95':>12}{'예측 q+c/d2':>14}{'실제오차':>12}"
      f"{'예측오차':>12}{'상대오차%':>11}{'예측%':>9}")
for d2 in (50, 200, 1000, 10_000, 100_000):
    act = d1 * stats.f(d1, d2).ppf(1 - alpha)
    pred = q + coef / d2
    print(f"{d2:>8}{act:>12.6f}{pred:>14.6f}{act - q:>12.6f}"
          f"{coef / d2:>12.6f}{100 * (act - q) / q:>11.4f}{100 * coef / (q * d2):>9.4f}")

print("\n  오차가 1/d2 로 줄어드는가 (d2 를 열 배 할 때 오차의 비)")
errs = [(d2, d1 * stats.f(d1, d2).ppf(1 - alpha) - q) for d2 in (100, 1_000, 10_000, 100_000)]
for (a, e1), (b, e2) in zip(errs, errs[1:]):
    print(f"    d2 {a:6d} -> {b:6d}:  오차 {e1:.6f} -> {e2:.6f}   비 {e1 / e2:.3f}")

(i) 역수 관계는 정확한 등식이다
  F_0.05( 5,20) = 0.21938814195492   1/F_0.95(20, 5) = 0.21938814195492   차 +5.6e-17
  F_0.1 (12,15) = 0.47509302181881   1/F_0.90(15,12) = 0.47509302181881   차 +0.0e+00
  F_0.25( 3,34) = 0.40549210425717   1/F_0.75(34, 3) = 0.40549210425717   차 -5.6e-17
  F_0.01( 4, 8) = 0.06757264103728   1/F_0.99( 8, 4) = 0.06757264103728   차 +0.0e+00

(ii) t^2 = F(1,d) 도 정확한 등식이다
  d =   1:  t_{d,0.975}^2 = 161.447638804129   F_{0.95}(1,d) = 161.447638797588   상대차 4.1e-11
  d =   2:  t_{d,0.975}^2 = 18.512820512362   F_{0.95}(1,d) = 18.512820512820   상대차 2.5e-11
  d =   5:  t_{d,0.975}^2 = 6.607890973703   F_{0.95}(1,d) = 6.607890973703   상대차 6.7e-16
  d =  10:  t_{d,0.975}^2 = 4.964602743636   F_{0.95}(1,d) = 4.964602743731   상대차 1.9e-11
  d =  34:  t_{d,0.975}^2 = 4.130017745652   F_{0.95}(1,d) = 4.130017745652   상대차 8.9e-16
  d = 200:  t_{d,0.975}^2 = 3.888374716773   F_{0.95}(1,d) = 3.888374716782   상대차 2.3e-12
  분포함수 항등식 2*T_10(sqrt(x)) - 1 = F(1,10).cdf(x) 의 최대오차: 2.2e-16

(iii) 카이제곱 극한만 근사다. 1차 보정항의 예측과 견준다
  d1 = 4, q = chi2_{0.95}(4) = 9.487729,  보정계수 q(q/2 - d1/2 + 1) = 35.5208
      d2   d1*F_0.95     예측 q+c/d2        실제오차        예측오차      상대오차%      예측%
      50   10.228717     10.198144    0.740988    0.710415     7.8100   7.4877
     200    9.667199      9.665333    0.179470    0.177604     1.8916   1.8719
    1000    9.523324      9.523250    0.035595    0.035521     0.3752   0.3744
   10000    9.491282      9.491281    0.003553    0.003552     0.0374   0.0374
  100000    9.488084      9.488084    0.000355    0.000355     0.0037   0.0037

  오차가 1/d2 로 줄어드는가 (d2 를 열 배 할 때 오차의 비)
    d2    100 ->   1000:  오차 0.362731 -> 0.035595   비 10.191
    d2   1000 ->  10000:  오차 0.035595 -> 0.003553   비 10.019
    d2  10000 -> 100000:  오차 0.003553 -> 0.000355   비 10.002