이표본 t-검정 (대마초 가격)¶

import numpy as np
from scipy import stats

np.random.seed(42)

CA_JAN2014 = np.array([
    248.75, 248.59, 248.63, 248.37, 248.02, 247.68, 247.36,
    246.85, 246.44, 246.06, 245.81, 245.48, 245.18, 244.87,
    244.55, 244.23, 243.89, 243.60, 243.34, 243.08, 242.85,
    242.64, 242.36, 242.15, 241.88, 241.64, 241.40, 241.14,
    240.91, 240.65, 240.42,
])

CA_JAN2015 = np.array([
    245.02, 244.88, 244.76, 244.65, 244.53, 244.42, 244.30,
    244.18, 244.08, 243.97, 243.85, 243.74, 243.63, 243.52,
    243.40, 243.28, 243.17, 243.06, 242.95, 242.83, 242.72,
    242.61, 242.49, 242.38, 242.27, 242.15, 242.04, 241.93,
    241.81, 241.70, 241.59,
])

for label, data in [("Jan 2014", CA_JAN2014), ("Jan 2015", CA_JAN2015)]:
    # 여기서는 H0가 "정규분포를 따른다"이다. 보통과 방향이 반대다.
    # 그래서 p가 크면 "정규성에 반하는 증거가 없다"는 뜻이고,
    # 이것을 "정규성이 확인되었다"로 읽으면 안 된다.
    stat, p = stats.shapiro(data)
    print(f"Shapiro-Wilk ({label}): W={stat:.4f}, p={p:.4f}")

Shapiro-Wilk (Jan 2014): W=0.9402, p=0.0835
Shapiro-Wilk (Jan 2015): W=0.9587, p=0.2694

n = len(CA_JAN2014)
print(f"n = {n}")
print(f"완전한 등차수열이라면 DW = 12/(n(n+1)) = {12 / (n * (n + 1)):.6f}, "
      f"지연 1 상관 = 1")
print(f"독립 표본이라면       DW ~ 2,  지연 1 상관 ~ 0 "
      f"(표준오차 약 {1 / np.sqrt(n):.3f})")

print(f"\n{'자료':>9} {'차분<0':>7} {'지연1 r':>9} {'DW':>9} "
      f"{'추세 r^2':>9} {'기울기':>9} {'잔차 sd':>9}")
for label, data in [("Jan 2014", CA_JAN2014), ("Jan 2015", CA_JAN2015)]:
    d = np.diff(data)
    r_lag = np.corrcoef(data[:-1], data[1:])[0, 1]
    resid = data - data.mean()
    dw = np.sum(np.diff(resid) ** 2) / np.sum(resid ** 2)
    t = np.arange(n)
    fit = stats.linregress(t, data)
    resid_tr = data - (fit.intercept + fit.slope * t)
    print(f"{label:>9} {int((d < 0).sum()):4d}/{len(d):<3d} {r_lag:9.6f} "
          f"{dw:9.6f} {fit.rvalue**2:9.6f} {fit.slope:9.4f} "
          f"{resid_tr.std(ddof=1):9.4f}")
print(f"\n전체 표준편차: 2014 {CA_JAN2014.std(ddof=1):.4f}, "
      f"2015 {CA_JAN2015.std(ddof=1):.4f}")

n = 31
완전한 등차수열이라면 DW = 12/(n(n+1)) = 0.012097, 지연 1 상관 = 1
독립 표본이라면       DW ~ 2,  지연 1 상관 ~ 0 (표준오차 약 0.180)

       자료    차분<0     지연1 r        DW    추세 r^2       기울기     잔차 sd
 Jan 2014   29/30   0.999396  0.011950  0.993389   -0.2926    0.2170
 Jan 2015   30/30   0.999976  0.012322  0.999952   -0.1135    0.0072

전체 표준편차: 2014 2.6689, 2015 1.0320

def confidence_interval(data, confidence=0.95):
    n = len(data)
    m = data.mean()
    # stats.sem은 ddof=1로 계산한 s를 sqrt(n)으로 나눈 표준오차다.
    se = stats.sem(data)
    h = se * stats.t.ppf((1 + confidence) / 2, n - 1)
    return m - h, m + h

ci14 = confidence_interval(CA_JAN2014)
ci15 = confidence_interval(CA_JAN2015)
print(f"Jan 2014 95% CI: [{ci14[0]:.2f}, {ci14[1]:.2f}]")
print(f"Jan 2015 95% CI: [{ci15[0]:.2f}, {ci15[1]:.2f}]")

Jan 2014 95% CI: [243.50, 245.46]
Jan 2015 95% CI: [242.91, 243.67]

def f_inflation(rho, n):
    """AR(1) 에서 Var(X-bar) 가 독립일 때의 몇 배인가."""
    k = np.arange(1, n)
    return 1 + 2 * np.sum((1 - k / n) * rho ** k)


print("  rho      f_n     n_eff = n/f_n    점근식 (1+rho)/(1-rho)")
for rho in (0.5, 0.9, 0.99, 0.999396, 0.999976):
    f = f_inflation(rho, n)
    print(f"{rho:9.6f} {f:8.4f} {n / f:12.4f} {(1 + rho) / (1 - rho):18.1f}")

print("\n2015 자료의 산포는 어디서 오는가")
print(f"  2015 전체 표준편차 {CA_JAN2015.std(ddof=1):.4f} 가운데")
t = np.arange(n)
fit15 = stats.linregress(t, CA_JAN2015)
resid15 = CA_JAN2015 - (fit15.intercept + fit15.slope * t)
print(f"    추세가 설명하는 몫 {fit15.rvalue**2 * 100:.4f}% (r^2),")
print(f"    잔차에 남는 표준편차 {resid15.std(ddof=1):.4f}")

# 포함률 모의실험: AR(1) 자료에 소박한 t 구간을 씌운다
tc = stats.t.ppf(0.975, n - 1)
rng = np.random.default_rng(5)
B = 40_000
print(f"\n정상 AR(1) 자료에서 소박한 95% t 구간의 포함률 (n = {n}, 반복 {B:,})")
print("  rho    이론 예측   모의 포함률")
for rho in (0.0, 0.5, 0.9, 0.99):
    f = f_inflation(rho, n)
    arg = tc * np.sqrt((n - f) / ((n - 1) * f))
    pred = 2 * stats.t.cdf(arg, n - 1) - 1
    e = rng.standard_normal((B, n))
    x = np.empty((B, n))
    x[:, 0] = e[:, 0]
    for j in range(1, n):
        x[:, j] = rho * x[:, j - 1] + np.sqrt(1 - rho ** 2) * e[:, j]
    m, s = x.mean(1), x.std(1, ddof=1)
    cov = np.mean(np.abs(m) <= tc * s / np.sqrt(n))
    print(f"{rho:6.2f}   {pred:9.4f}   {cov:11.4f}")

  rho      f_n     n_eff = n/f_n    점근식 (1+rho)/(1-rho)
 0.500000   2.8710      10.7978                3.0
 0.900000  13.4151       2.3108               19.0
 0.990000  28.0196       1.1064              199.0
 0.999396  30.8076       1.0062             3310.3
 0.999976  30.9923       1.0002            83332.3

2015 자료의 산포는 어디서 오는가
  2015 전체 표준편차 1.0320 가운데
    추세가 설명하는 몫 99.9952% (r^2),
    잔차에 남는 표준편차 0.0072

정상 AR(1) 자료에서 소박한 95% t 구간의 포함률 (n = 31, 반복 40,000)
  rho    이론 예측   모의 포함률
  0.00      0.9500        0.9512
  0.50      0.7476        0.7414
  0.90      0.3275        0.3104
  0.99      0.0960        0.0917

# equal_var=True는 합동 t-검정이다. scipy의 기본값이기도 하다.
t_stat, p_val = stats.ttest_ind(CA_JAN2014, CA_JAN2015, equal_var=True)
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value:     {p_val:.6f}")
print("Reject H0" if p_val < 0.05 else "Fail to reject H0")

t-statistic: 2.3167
p-value:     0.023951
Reject H0

s1 = CA_JAN2014.std(ddof=1)
s2 = CA_JAN2015.std(ddof=1)
d1 = d2 = n - 1
se_w = np.sqrt(s1**2 / n + s2**2 / n)
sp2 = (d1 * s1**2 + d2 * s2**2) / (d1 + d2)
se_p = np.sqrt(sp2 * (2 / n))
diff = CA_JAN2014.mean() - CA_JAN2015.mean()
print(f"분산비 s1^2/s2^2 = {s1**2 / s2**2:.6f}  (표준편차로는 {s1 / s2:.6f} 배)")
print(f"SE_welch = {se_w:.12f}")
print(f"SE_pool  = {se_p:.12f}   차 = {se_p - se_w:.3e}")
print(f"t_welch = {diff / se_w:.12f}")
print(f"t_pool  = {diff / se_p:.12f}")

nu = (s1**2 / n + s2**2 / n) ** 2 / ((s1**2 / n) ** 2 / d1 + (s2**2 / n) ** 2 / d2)
print(f"\n자유도: 합동 {2 * n - 2},  Welch {nu:.6f}")
print(f"p-값:   합동 {2 * stats.t.sf(diff / se_p, 2 * n - 2):.6f},  "
      f"Welch {2 * stats.t.sf(diff / se_w, nu):.6f}")

# 균형 설계에서 합동 t 의 실제 수준 (5.3 절의 m)
rho = s2**2 / s1**2
m = d1 * (1 + rho) ** 2 / (1 + rho**2)
crit = stats.t.ppf(0.975, 2 * n - 2)
print(f"\n관측된 s 를 sigma 로 삼으면 m = {m:.6f}  (R = 1, 균형이므로)")
print(f"다듬은 어림 P(|t_m| > {crit:.6f}) = {2 * stats.t.sf(crit, m):.6f}")

rng = np.random.default_rng(23)
B = 200_000
x = rng.normal(0, s1, (B, n))
y = rng.normal(0, s2, (B, n))
dd = x.mean(1) - y.mean(1)
v1, v2 = x.var(1, ddof=1), y.var(1, ddof=1)
t_pool_sim = dd / np.sqrt(((d1 * v1 + d2 * v2) / (2 * n - 2)) * (2 / n))
a, b = v1 / n, v2 / n
t_welch_sim = dd / np.sqrt(a + b)
nu_sim = (a + b) ** 2 / (a**2 / d1 + b**2 / d2)
print(f"\n반복 {B:,} 회 (오차 약 ±0.0005)")
print(f"  합동 t 의 실제 수준  = {np.mean(np.abs(t_pool_sim) > crit):.5f}")
print(f"  Welch 의 실제 수준   = "
      f"{np.mean(np.abs(t_welch_sim) > stats.t(nu_sim).ppf(0.975)):.5f}")
print(f"  실현된 nu 의 평균 {nu_sim.mean():.4f}, "
      f"범위 [{nu_sim.min():.4f}, {nu_sim.max():.4f}]  (한계 {n - 1} 과 {2 * n - 2})")

분산비 s1^2/s2^2 = 6.687940  (표준편차로는 2.586105 배)
SE_welch = 0.513938489917
SE_pool  = 0.513938489917   차 = 0.000e+00
t_welch = 2.316707513931
t_pool  = 2.316707513931

자유도: 합동 60,  Welch 38.775185
p-값:   합동 0.023951,  Welch 0.025894

관측된 s 를 sigma 로 삼으면 m = 38.775185  (R = 1, 균형이므로)
다듬은 어림 P(|t_m| > 2.000298) = 0.052507

반복 200,000 회 (오차 약 ±0.0005)
  합동 t 의 실제 수준  = 0.05265
  Welch 의 실제 수준   = 0.05015
  실현된 nu 의 평균 39.2630, 범위 [31.6778, 59.8803]  (한계 30 과 60)

counts_2014 = np.array([453020, 688699, 271937])
observed_2015 = np.array([461900, 695432, 267120])

# 2014년 도수를 2015년 총합에 맞춰 축척한다.
# chisquare는 sum(E) == sum(O)를 요구하며, 이 줄을 빠뜨리면
# "총합이 다르다"는 사실이 통계량에 섞여 들어가 무의미한 값이 나온다.
expected_2015 = counts_2014 / counts_2014.sum() * observed_2015.sum()

chi2, p_chi = stats.chisquare(observed_2015, f_exp=expected_2015)
print(f"Chi-square stat: {chi2:.2f}, p-value: {p_chi:.6f}")

Chi-square stat: 240.93, p-value: 0.000000

p14 = counts_2014 / counts_2014.sum()
p15 = observed_2015 / observed_2015.sum()
N = observed_2015.sum()
labels = ("고품질", "중품질", "저품질")

print(f"N(2014) = {counts_2014.sum():,},  N(2015) = {N:,}")
print(f"{'등급':>7} {'2014 비율':>10} {'2015 비율':>10} {'차(%p)':>9} "
      f"{'기대도수':>13} {'기여':>10} {'몫':>7}")
contrib = (observed_2015 - expected_2015) ** 2 / expected_2015
for lab, a, b, e, c in zip(labels, p14, p15, expected_2015, contrib):
    print(f"{lab:>7} {a:10.6f} {b:10.6f} {(b - a) * 100:9.4f} "
          f"{e:13.2f} {c:10.4f} {c / contrib.sum():7.4f}")
print(f"{'합':>7} {p14.sum():10.6f} {p15.sum():10.6f} {'':9} "
      f"{expected_2015.sum():13.2f} {contrib.sum():10.4f} {1.0:7.4f}")

w2 = np.sum((p15 - p14) ** 2 / p14)
print(f"\nw^2 = sum (p15-p14)^2/p14 = {w2:.9f},  w = {np.sqrt(w2):.6f}")
print(f"N * w^2 = {N * w2:.6f}   (= chi2 = {chi2:.6f})")
print(f"크라메르 V = w/sqrt(k-1) = {np.sqrt(w2 / 2):.6f}")
print(f"자유도 = k - 1 = {len(p14) - 1}  "
      f"(기대비율을 2014 에서 가져왔으므로 추정한 모수가 없다)")
print(f"기각에 필요한 N = chi2_(0.95,2)/w^2 = "
      f"{stats.chi2.ppf(0.95, 2) / w2:.1f}")

N(2014) = 1,413,656,  N(2015) = 1,424,452
     등급    2014 비율    2015 비율     차(%p)          기대도수         기여       몫
    고품질   0.320460   0.324265    0.3805     456479.68    64.3617  0.2671
    중품질   0.487176   0.488210    0.1034     693958.55     3.1285  0.0130
    저품질   0.192364   0.187525   -0.4840     274013.77   173.4365  0.7199
      합   1.000000   1.000000              1424452.00   240.9268  1.0000

w^2 = sum (p15-p14)^2/p14 = 0.000169136,  w = 0.013005
N * w^2 = 240.926777   (= chi2 = 240.926777)
크라메르 V = w/sqrt(k-1) = 0.009196
자유도 = k - 1 = 2  (기대비율을 2014 에서 가져왔으므로 추정한 모수가 없다)
기각에 필요한 N = chi2_(0.95,2)/w^2 = 35423.8

ci99 = confidence_interval(CA_JAN2014, confidence=0.99)
print(f"99% CI: [{ci99[0]:.2f}, {ci99[1]:.2f}]")

99% CI: [243.16, 245.80]

t_w, p_w = stats.ttest_ind(CA_JAN2014, CA_JAN2015, equal_var=False)
print(f"Welch t = {t_w:.4f}, p = {p_w:.6f}")

Welch t = 2.3167, p = 0.025894

import numpy as np
from scipy import stats

rng = np.random.default_rng(1717)
M, n = 20_000, 31

def ci(d, conf):
    se = stats.sem(d)
    h = se * stats.t.ppf((1 + conf) / 2, len(d) - 1)
    return d.mean() - h, d.mean() + h

def disjoint(x, y, conf):
    l1, u1 = ci(x, conf)
    l2, u2 = ci(y, conf)
    return not (l1 <= u2 and l2 <= u1)

for shift, label in [(0.0, "수준 (참 차이 0)"), (0.8, "검정력 (참 차이 0.8)")]:
    a = b = c = 0
    for _ in range(M):
        x, y = rng.standard_normal(n) + shift, rng.standard_normal(n)
        a += disjoint(x, y, 0.95)
        b += disjoint(x, y, 0.84)
        c += stats.ttest_ind(x, y, equal_var=False).pvalue < 0.05
    print(f"{label}")
    print(f"  95% 구간 비겹침 {a / M:.4f}")
    print(f"  84% 구간 비겹침 {b / M:.4f}")
    print(f"  Welch t 기각률  {c / M:.4f}")

수준 (참 차이 0)
  95% 구간 비겹침 0.0053
  84% 구간 비겹침 0.0461
  Welch t 기각률  0.0490
검정력 (참 차이 0.8)
  95% 구간 비겹침 0.6068
  84% 구간 비겹침 0.8649
  Welch t 기각률  0.8706

import numpy as np
from scipy import stats

counts_2014 = np.array([453020, 688699, 271937])
observed_2015 = np.array([461900, 695432, 267120])
expected = counts_2014 / counts_2014.sum() * observed_2015.sum()
chi2, p = stats.chisquare(observed_2015, f_exp=expected)
N = observed_2015.sum()
k = len(observed_2015)

print(f"χ² = {chi2:.2f},  p = {p:.3g},  N = {N:,}")
print(f"코헨 w      = {np.sqrt(chi2 / N):.5f}")
print(f"크라메르 V  = {np.sqrt(chi2 / (N * (k - 1))):.5f}")
print(f"2014 비율 {np.round(counts_2014 / counts_2014.sum(), 5)}")
print(f"2015 비율 {np.round(observed_2015 / N, 5)}")
print(f"차이(%p)  {np.round((observed_2015 / N - counts_2014 / counts_2014.sum()) * 100, 3)}")

print("\n같은 비율을 유지하고 표본만 줄이면")
for scale in [1, 0.1, 0.01, 0.005, 0.001]:
    o = np.round(observed_2015 * scale)
    e = counts_2014 / counts_2014.sum() * o.sum()
    cs, ps = stats.chisquare(o, f_exp=e)
    print(f"  표본 {o.sum():>9,.0f}: χ² = {cs:8.3f},  p = {ps:.4f}")

for N0 in range(100, 200_000, 100):
    o = np.round(observed_2015 / N * N0)
    e = counts_2014 / counts_2014.sum() * o.sum()
    cs, ps = stats.chisquare(o, f_exp=e)
    if ps < 0.05:
        print(f"\n유의해지는 최소 표본 ≈ {o.sum():,.0f} "
              f"(χ² = {cs:.3f}, p = {ps:.4f})")
        break

χ² = 240.93,  p = 4.82e-53,  N = 1,424,452
코헨 w      = 0.01301
크라메르 V  = 0.00920
2014 비율 [0.32046 0.48718 0.19236]
2015 비율 [0.32427 0.48821 0.18752]
차이(%p)  [ 0.381  0.103 -0.484]

같은 비율을 유지하고 표본만 줄이면
  표본 1,424,452: χ² =  240.927,  p = 0.0000
  표본   142,445: χ² =   24.092,  p = 0.0000
  표본    14,244: χ² =    2.418,  p = 0.2985
  표본     7,123: χ² =    1.196,  p = 0.5500
  표본     1,424: χ² =    0.248,  p = 0.8835

유의해지는 최소 표본 ≈ 35,400 (χ² = 6.008, p = 0.0496)

import numpy as np
from scipy import stats

CA_JAN2014 = np.array([
    248.75, 248.59, 248.63, 248.37, 248.02, 247.68, 247.36,
    246.85, 246.44, 246.06, 245.81, 245.48, 245.18, 244.87,
    244.55, 244.23, 243.89, 243.60, 243.34, 243.08, 242.85,
    242.64, 242.36, 242.15, 241.88, 241.64, 241.40, 241.14,
    240.91, 240.65, 240.42,
])
CA_JAN2015 = np.array([
    245.02, 244.88, 244.76, 244.65, 244.53, 244.42, 244.30,
    244.18, 244.08, 243.97, 243.85, 243.74, 243.63, 243.52,
    243.40, 243.28, 243.17, 243.06, 242.95, 242.83, 242.72,
    242.61, 242.49, 242.38, 242.27, 242.15, 242.04, 241.93,
    241.81, 241.70, 241.59,
])

for label, d in [("2014", CA_JAN2014), ("2015", CA_JAN2015)]:
    n = len(d)
    r1 = np.corrcoef(d[:-1], d[1:])[0, 1]
    diffs = np.diff(d)
    slope, _, r, p, _ = stats.linregress(np.arange(n), d)
    print(f"{label}: 평균 {d.mean():.4f}  SD {d.std(ddof=1):.4f}")
    print(f"   1차 자기상관 {r1:.4f}")
    print(f"   일별 변화   평균 {diffs.mean():+.4f}  SD {diffs.std(ddof=1):.4f}  "
          f"하락 비율 {np.mean(diffs < 0):.3f}")
    print(f"   선형추세    기울기 {slope:+.4f}/일  R² = {r**2:.4f}  p = {p:.2e}")
    print(f"   유효 표본크기 n_eff ≈ n(1-r)/(1+r) = {n * (1 - r1) / (1 + r1):.2f}")

2014: 평균 244.4781  SD 2.6689
   1차 자기상관 0.9994
   일별 변화   평균 -0.2777  SD 0.0911  하락 비율 0.967
   선형추세    기울기 -0.2926/일  R² = 0.9934  p = 3.65e-33
   유효 표본크기 n_eff ≈ n(1-r)/(1+r) = 0.01
2015: 평균 243.2874  SD 1.0320
   1차 자기상관 1.0000
   일별 변화   평균 -0.1143  SD 0.0073  하락 비율 1.000
   선형추세    기울기 -0.1135/일  R² = 1.0000  p = 3.85e-64
   유효 표본크기 n_eff ≈ n(1-r)/(1+r) = 0.00

import numpy as np
from scipy import stats

rng = np.random.default_rng(3141)
M, n = 10_000, 31

gens = {"정규":     lambda k: rng.standard_normal(k),
        "t(3)":     lambda k: rng.standard_t(3, k),
        "지수":     lambda k: rng.exponential(1, k),
        "로그정규": lambda k: rng.lognormal(0, 1, k)}

print(f"{'분포':>10s} {'':>6s} {'Shapiro 통과':>13s} {'항상 t':>8s} "
      f"{'항상 MWU':>10s} {'2단계':>8s}")
for name, gen in gens.items():
    rows = []
    for shift in [0.0, 0.8]:
        passed = a = b = c = 0
        for _ in range(M):
            x, y = gen(n), gen(n) + shift
            ok = (stats.shapiro(x).pvalue > 0.05
                  and stats.shapiro(y).pvalue > 0.05)
            passed += ok
            p_t = stats.ttest_ind(x, y, equal_var=False).pvalue
            p_m = stats.mannwhitneyu(x, y).pvalue
            a += p_t < 0.05
            b += p_m < 0.05
            c += (p_t if ok else p_m) < 0.05
        rows.append((passed / M, a / M, b / M, c / M))
    print(f"{name:>10s} {'수준':>6s} {rows[0][0]:13.4f} {rows[0][1]:8.4f} "
          f"{rows[0][2]:10.4f} {rows[0][3]:8.4f}")
    print(f"{'':>10s} {'검정력':>6s} {rows[1][0]:13.4f} {rows[1][1]:8.4f} "
          f"{rows[1][2]:10.4f} {rows[1][3]:8.4f}")

        분포           Shapiro 통과     항상 t     항상 MWU      2단계
        정규     수준        0.9010   0.0507     0.0489   0.0517
              검정력        0.9013   0.8677     0.8510   0.8668
      t(3)     수준        0.2848   0.0470     0.0512   0.0521
              검정력        0.2861   0.5094     0.6725   0.6591
        지수     수준        0.0003   0.0484     0.0480   0.0481
              검정력        0.0005   0.8666     0.9750   0.9750
      로그정규     수준        0.0000   0.0429     0.0501   0.0501
              검정력        0.0000   0.4090     0.8787   0.8787

① 독립성 — 자료 수집 방식을 본다. 검정으로 확인할 수 없다.
     ↓ (통과해야 다음으로)
② 등분산 — 확인하지 말고 웰치를 쓴다.
     ↓
③ 정규성 — 그림으로 본다. 검정은 보조.