모분산에 대한 카이제곱 검정¶

import numpy as np
from scipy import stats

LSL, USL = 9.5, 10.5
n, xbar, s = 30, 10.05, 0.12

Cp = (USL - LSL) / (6 * s)
Cpk = min((USL - xbar) / (3 * s), (xbar - LSL) / (3 * s))
print(f"Ĉp = {Cp:.4f},  Ĉpk = {Cpk:.4f}")

# Cp 의 신뢰구간: σ 의 카이제곱 구간을 뒤집는다
nu = n - 1
lo = Cp * np.sqrt(stats.chi2.ppf(0.025, nu) / nu)
hi = Cp * np.sqrt(stats.chi2.ppf(0.975, nu) / nu)
print(f"Cp 의 95% 구간 ({lo:.4f}, {hi:.4f})")

# 규격 이탈 비율
below = stats.norm.cdf(LSL, xbar, s)
above = stats.norm.sf(USL, xbar, s)
print(f"추정 불량률 {1e6 * (below + above):.1f} ppm")

Ĉp = 1.3889,  Ĉpk = 1.2500
Cp 의 95% 구간 (1.0332, 1.7439)
추정 불량률 90.7 ppm

for n in [30, 50, 100, 200, 500]:
    nu = n - 1
    # Ĉp 가 얼마여야 구간 하한이 1.33 을 넘는가
    need = 1.33 / np.sqrt(stats.chi2.ppf(0.025, nu) / nu)
    print(f"n={n:4d}: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ {need:.3f}")

n=  30: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ 1.788
n=  50: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ 1.657
n= 100: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ 1.545
n= 200: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ 1.475
n= 500: 구간 하한이 1.33 을 넘으려면 Ĉp ≥ 1.418

import numpy as np
from scipy import stats

rng = np.random.default_rng(72)
M, k, n = 3_000, 4, 20          # 집단 4개, 각 20개

print(f"{'분포':>7s} {'바틀렛':>9s} {'레빈(평균)':>11s} "
      f"{'BF(중앙값)':>11s}")
for name, gen in [("정규", lambda s: rng.normal(0, 1, s)),
                  ("t(5)", lambda s: rng.standard_t(5, s)),
                  ("지수", lambda s: rng.exponential(1, s)),
                  ("균등", lambda s: rng.uniform(-1.7, 1.7, s))]:
    c = np.zeros(3)
    for _ in range(M):
        g = [gen(n) for _ in range(k)]
        c[0] += stats.bartlett(*g).pvalue < 0.05
        c[1] += stats.levene(*g, center="mean").pvalue < 0.05
        c[2] += stats.levene(*g, center="median").pvalue < 0.05
    print(f"{name:>7s} {c[0] / M:9.4f} {c[1] / M:11.4f} "
          f"{c[2] / M:11.4f}")

     분포       바틀렛      레빈(평균)     BF(중앙값)
     정규    0.0460      0.0553      0.0327
   t(5)    0.2440      0.0523      0.0340
     지수    0.4837      0.2233      0.0513
     균등    0.0000      0.0557      0.0227

import numpy as np
from scipy import stats

def k2(n, P=0.95, conf=0.95):
    """양측 허용구간의 계수 (하워 근사)"""
    nu = n - 1
    z = stats.norm.ppf(0.5 + P / 2)
    chi = stats.chi2.ppf(1 - conf, nu)
    return z * np.sqrt(nu * (1 + 1 / n) / chi)

n, xbar, s = 30, 10.05, 0.12
t = stats.t.ppf(0.975, n - 1)
print(f"신뢰구간(μ)     ±{t * s / np.sqrt(n):.4f}  → "
      f"({xbar - t * s / np.sqrt(n):.4f}, {xbar + t * s / np.sqrt(n):.4f})")
hp = t * s * np.sqrt(1 + 1 / n)
print(f"예측구간(다음 1개) ±{hp:.4f}  → "
      f"({xbar - hp:.4f}, {xbar + hp:.4f})")
for P in [0.90, 0.95, 0.99]:
    k = k2(n, P)
    print(f"허용구간(P={P:.2f})  ±{k * s:.4f}  → "
          f"({xbar - k * s:.4f}, {xbar + k * s:.4f})")

신뢰구간(μ)     ±0.0448  → (10.0052, 10.0948)
예측구간(다음 1개) ±0.2495  → (9.8005, 10.2995)
허용구간(P=0.90)  ±0.2568  → (9.7932, 10.3068)
허용구간(P=0.95)  ±0.3060  → (9.7440, 10.3560)
허용구간(P=0.99)  ±0.4021  → (9.6479, 10.4521)

for n in [10, 30, 100, 1000]:
    print(f"n={n:5d}: k2(95%, 95%) = {k2(n):.3f}  "
          f"(n→∞ 이면 z = {stats.norm.ppf(0.975):.3f})")

n=   10: k2(95%, 95%) = 3.382  (n→∞ 이면 z = 1.960)
n=   30: k2(95%, 95%) = 2.550  (n→∞ 이면 z = 1.960)
n=  100: k2(95%, 95%) = 2.233  (n→∞ 이면 z = 1.960)
n= 1000: k2(95%, 95%) = 2.036  (n→∞ 이면 z = 1.960)

import numpy as np

rng = np.random.default_rng(11)
n_part, n_op, n_rep = 10, 3, 3
sd_part, sd_op, sd_rep = 1.0, 0.25, 0.30

part = rng.normal(0, sd_part, n_part)
op = rng.normal(0, sd_op, n_op)
y = np.empty((n_part, n_op, n_rep))
for i in range(n_part):
    for j in range(n_op):
        y[i, j] = part[i] + op[j] + rng.normal(0, sd_rep, n_rep)

# 분산분석식 적률법
grand = y.mean()
ms_part = n_op * n_rep * ((y.mean((1, 2)) - grand)**2).sum() / (n_part - 1)
ms_op = n_part * n_rep * ((y.mean((0, 2)) - grand)**2).sum() / (n_op - 1)
cell = y.mean(2)
ms_int = (n_rep * ((cell - y.mean((1, 2))[:, None]
                    - y.mean((0, 2))[None, :] + grand)**2).sum()
          / ((n_part - 1) * (n_op - 1)))
ms_err = ((y - cell[:, :, None])**2).sum() / (n_part * n_op * (n_rep - 1))

v_rep = ms_err
v_op = max((ms_op - ms_int) / (n_part * n_rep), 0)
v_int = max((ms_int - ms_err) / n_rep, 0)
v_part = max((ms_part - ms_int) / (n_op * n_rep), 0)
v_grr = v_rep + v_op + v_int
v_tot = v_grr + v_part

print(f"{'성분':>10s} {'분산':>9s} {'표준편차':>9s} {'총 변동 비중':>12s}")
for name, v in [("반복성", v_rep), ("재현성(검사자)", v_op),
                ("상호작용", v_int), ("GRR 합계", v_grr),
                ("부품", v_part), ("총계", v_tot)]:
    print(f"{name:>10s} {v:9.4f} {np.sqrt(v):9.4f} "
          f"{100 * v / v_tot:11.1f}%")
print(f"\n%GRR (표준편차 기준) = "
      f"{100 * np.sqrt(v_grr / v_tot):.1f}%")

      성분        분산      표준편차       총 변동 비중
     반복성    0.0696    0.2638         7.2%
재현성(검사자)    0.0259    0.1610         2.7%
    상호작용    0.0066    0.0810         0.7%
  GRR 합계    0.1021    0.3195        10.5%
      부품    0.8676    0.9315        89.5%
      총계    0.9697    0.9847       100.0%

import numpy as np
from scipy import stats

def power_upper(n, r, alpha=0.05):
    """상단 단측: H0: σ² ≤ σ0², 참 비율이 r 일 때의 검정력"""
    nu = n - 1
    return stats.chi2.sf(stats.chi2.ppf(1 - alpha, nu) / r, nu)

print(f"{'r = σ²/σ0²':>11s} {'검정력 80% 에 필요한 n':>22s} "
      f"{'검정력 90%':>12s}")
for r in [1.25, 1.5, 2.0, 2.5, 3.0]:
    n80 = next(m for m in range(5, 5000) if power_upper(m, r) >= 0.80)
    n90 = next(m for m in range(5, 5000) if power_upper(m, r) >= 0.90)
    print(f"{r:11.2f} {n80:22d} {n90:12d}")

 r = σ²/σ0²        검정력 80% 에 필요한 n      검정력 90%
       1.25                    245          343
       1.50                     74          105
       2.00                     26           37
       2.50                     16           22
       3.00                     11           16

for n in [30, 50, 100, 200, 500, 1000]:
    nu = n - 1
    need = 1.33 / np.sqrt(stats.chi2.ppf(0.025, nu) / nu)
    print(f"n={n:5d}: Ĉp ≥ {need:.3f} 이어야 'Cp ≥ 1.33' 을 주장 가능")

n=   30: Ĉp ≥ 1.788 이어야 'Cp ≥ 1.33' 을 주장 가능
n=   50: Ĉp ≥ 1.657 이어야 'Cp ≥ 1.33' 을 주장 가능
n=  100: Ĉp ≥ 1.545 이어야 'Cp ≥ 1.33' 을 주장 가능
n=  200: Ĉp ≥ 1.475 이어야 'Cp ≥ 1.33' 을 주장 가능
n=  500: Ĉp ≥ 1.418 이어야 'Cp ≥ 1.33' 을 주장 가능
n= 1000: Ĉp ≥ 1.391 이어야 'Cp ≥ 1.33' 을 주장 가능