유도와 카이제곱 근사¶

import numpy as np
from scipy import stats

# 원자료 없이 집단 크기와 표본분산만으로 통계량을 만들 수 있다.
n = np.array([10, 12, 8])
s2 = np.array([5.2, 8.1, 4.7])
k = len(n)
nu = n - 1
N = n.sum()

# 귀무가설 아래의 공통 분산 추정값
s2_pooled = np.sum(nu * s2) / np.sum(nu)

# 분자는 합동분산의 로그와 각 분산 로그의 가중평균의 차이다.
# 분산들이 고를수록 0 에 가까워진다.
numerator = np.sum(nu) * np.log(s2_pooled) - np.sum(nu * np.log(s2))

# 보정인자. 표본이 작을수록 1 보다 눈에 띄게 커져 통계량을 낮춘다.
C = 1 + (1 / (3 * (k - 1))) * (np.sum(1 / nu) - 1 / np.sum(nu))

# 검정통계량
T = numerator / C

# p-값
p_value = stats.chi2.sf(T, k - 1)

print(f"Pooled variance: {s2_pooled:.4f}")
print(f"Numerator: {numerator:.4f}")
print(f"Correction factor C: {C:.4f}")
print(f"Test statistic T: {T:.4f}")
print(f"P-value: {p_value:.4f}")


def bart_from_summary(n, s2):
    """(n_i, S_i^2) 만으로 바틀렛 통계량과 p 값을 돌려준다."""
    n, s2 = np.asarray(n, float), np.asarray(s2, float)
    k, nu = len(n), np.asarray(n, float) - 1
    sp2 = np.sum(nu * s2) / np.sum(nu)
    num = np.sum(nu) * np.log(sp2) - np.sum(nu * np.log(s2))
    C = 1 + (1 / (3 * (k - 1))) * (np.sum(1 / nu) - 1 / np.sum(nu))
    return num / C, stats.chi2.sf(num / C, k - 1)


# 요약값에 맞는 원자료를 아무렇게나 만들어 scipy 에 넣어 본다.
# 평균을 어디에 두든, 어떤 분포에서 뽑든 분산만 맞추면 같은 값이 나와야 한다.
def fake(n, s2, rng, draw):
    out = []
    for ni, v in zip(n, s2):
        z = draw(ni)
        z = (z - z.mean()) / z.std(ddof=1) * np.sqrt(v)
        out.append(z + rng.normal() * 10)
    return out


rng = np.random.default_rng(99)
g_norm = fake(n, s2, rng, lambda m: rng.normal(size=m))
rng2 = np.random.default_rng(7)
g_t3 = [a - 123.0 for a in fake(n, s2, rng2, lambda m: rng2.standard_t(3, size=m))]

print(f"\n만든 자료의 표본분산: {[round(float(g.var(ddof=1)), 10) for g in g_norm]}")
print(f"{'요약값만으로':22}{bart_from_summary(n, s2)[0]:.13f}")
print(f"{'정규에서 만든 자료':22}{stats.bartlett(*g_norm).statistic:.13f}")
print(f"{'t_3 에서 만든 자료':22}{stats.bartlett(*g_t3).statistic:.13f}")

# 척도불변성: 분산에 공통 상수를 곱해도 T 가 바뀌지 않는다.
print(f"\n분산을 1000배 하면:  T = {bart_from_summary(n, s2 * 1000)[0]:.13f}")

# ddof 함정. n 으로 나눈 분산을 받아 그대로 넣으면?
s2_pop = s2 * nu / n
print(f"\nn 으로 나눈 분산 = {np.round(s2_pop, 4)}")
print(f"  불균형 n=(10,12,8):  ddof=1 -> T = {bart_from_summary(n, s2)[0]:.4f}"
      f"   ddof=0 -> T = {bart_from_summary(n, s2_pop)[0]:.4f}")
ne = np.array([10, 10, 10])
print(f"  균형   n=(10,10,10): ddof=1 -> T = {bart_from_summary(ne, s2)[0]:.13f}")
print(f"                       ddof=0 -> T = {bart_from_summary(ne, s2 * 9 / 10)[0]:.13f}")

Pooled variance: 6.2519
Numerator: 0.8063
Correction factor C: 1.0513
Test statistic T: 0.7670
P-value: 0.6815

만든 자료의 표본분산: [5.2, 8.1, 4.7]
요약값만으로                0.7669773746398
정규에서 만든 자료            0.7669773746398
t_3 에서 만든 자료          0.7669773746398

분산을 1000배 하면:  T = 0.7669773746398

n 으로 나눈 분산 = [4.68   7.425  4.1125]
  불균형 n=(10,12,8):  ddof=1 -> T = 0.7670   ddof=0 -> T = 0.8743
  균형   n=(10,10,10): ddof=1 -> T = 0.7478100156592
                       ddof=0 -> T = 0.7478100156592

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
k, n, R, alpha = 3, 5, 20000, 0.05
nu = n - 1
N = k * n
C = 1 + (1 / (3 * (k - 1))) * (k / nu - 1 / (N - k))
crit = stats.chi2.ppf(1 - alpha, k - 1)

rej_raw = rej_corr = 0
for _ in range(R):
    s2 = np.array([rng.normal(0, 1, n).var(ddof=1) for _ in range(k)])
    sp2 = s2.mean()                      # equal nu, so simple mean
    raw = (N - k) * np.log(sp2) - nu * np.sum(np.log(s2))
    rej_raw += (raw > crit)
    rej_corr += (raw / C > crit)

print(f"Correction factor C = {C:.4f}")
print(f"Uncorrected size: {rej_raw / R:.4f}")
print(f"Corrected size:   {rej_corr / R:.4f}")

Correction factor C = 1.1111
Uncorrected size: 0.0672
Corrected size:   0.0505

import numpy as np
from scipy import stats

n = np.array([10, 12, 8])
k, nu, N = 3, n - 1, n.sum()
C = 1 + (1 / (3 * (k - 1))) * (np.sum(1 / nu) - 1 / np.sum(nu))
crit = stats.chi2.ppf(0.95, k - 1)

print(f"{'s2_2':>8} {'T':>8} {'p':>8} {'ratio to others':>18}")
for s2_2 in [8.1, 15, 21, 25, 40, 60, 80]:
    s2 = np.array([5.2, s2_2, 4.7])
    sp2 = np.sum(nu * s2) / np.sum(nu)
    T = (np.sum(nu) * np.log(sp2) - np.sum(nu * np.log(s2))) / C
    p = stats.chi2.sf(T, k - 1)
    print(f"{s2_2:>8.1f} {T:>8.3f} {p:>8.4f} {s2_2 / 5.0:>18.1f}")

    s2_2        T        p    ratio to others
     8.1    0.767   0.6815                1.6
    15.0    3.856   0.1454                3.0
    21.0    6.468   0.0394                4.2
    25.0    8.045   0.0179                5.0
    40.0   12.938   0.0016                8.0
    60.0   17.761   0.0001               12.0
    80.0   21.438   0.0000               16.0