ADP 실기 — 통계 검정, 가설부터 쓰고 시작한다

ADP 실기의 통계 문항은 배점의 절반이다. 그런데 코드는 대개 한 줄이다. 점수는 그 앞뒤에 있다 — 가설을 세우고, 통계량을 구하고, 채택 여부를 문장으로 쓰는 것까지가 한 문제다.

23회, 24회, 25회, 26회 모두 문제 형식이 똑같았다.

n-1. 연구가설과 귀무가설을 설정하시오        (5점)
n-2. 검정통계량을 구하고 가설을 채택하시오    (5점)

그래서 어떤 검정이든 이 네 줄을 먼저 쓰고 코드를 짰다.

1. 가설 설정 — $H_0$ 와 $H_1$
2. 유의수준 설정 — 보통 $\alpha = 0.05$
3. 검정통계량과 유의확률 계산
4. 기각 여부 판정 — $p < \alpha$ 면 귀무가설을 기각한다

표본 크기

26회에 “불량률이 0.9일 때 오차의 한계가 5%가 되도록 하는 최소 표본 크기”가 나왔다.

$n = \frac{Z_{\alpha/2}^2 \; p(1-p)}{E^2}$

기호	뜻
$Z_{\alpha/2}$	신뢰수준에 해당하는 값. 95% 면 1.96
$p$	모비율(불량률)의 추정치
$E$	허용 오차의 한계

from scipy.stats import norm
import math

p, E, alpha = 0.9, 0.05, 0.05
z = norm.ppf(1 - alpha/2)              # 1.959963...
n = z**2 * p * (1-p) / E**2
print(n, math.ceil(n))                 # 138.29... → 139

올림한다. 138.3 명을 뽑을 수는 없고, 내림하면 요구한 정밀도에 못 미친다. $p$ 를 모르면 $p = 0.5$ 를 넣는다 — $p(1-p)$ 가 그때 최대라서 가장 보수적인 표본 크기가 된다.

모분산의 신뢰구간

25회 문제 — 표본 10개의 분산이 90일 때 신뢰도 95%로 모분산의 신뢰구간을 추정하시오.

$\frac{(n-1)s^2}{\chi^2_{\alpha/2,\,n-1}} \le \sigma^2 \le \frac{(n-1)s^2}{\chi^2_{1-\alpha/2,\,n-1}}$

from scipy.stats import chi2

n, s2, alpha = 10, 90, 0.05
lower = (n-1) * s2 / chi2.ppf(1 - alpha/2, n-1)   # 42.58
upper = (n-1) * s2 / chi2.ppf(alpha/2, n-1)       # 299.96
print(lower, upper)

분모가 뒤집혀 있다. 큰 카이제곱 값으로 나눈 쪽이 하한이다. 평균의 신뢰구간처럼 “추정치 ± 무엇” 꼴이 아니라서 여기서 자주 틀린다. 그리고 구간이 42.6~300.0 으로 어이없이 넓은데, 이게 정상이다 — 분산 추정은 표본이 10개면 이 정도로 못 미덥다는 게 이 문제가 보여 주려는 것이다.

t 검정 세 가지

from scipy import stats

# 일표본 — 모평균이 2.3인가
t = stats.ttest_1samp(data, 2.3, alternative='less')

# 대응표본 — 같은 대상의 전후 (혈압약 복용 전후)
pt = stats.ttest_rel(before, after)

# 독립표본 — 다른 두 집단 (남학생 vs 여학생)
it = stats.ttest_ind(group_a, group_b)

	언제
일표본	표본 평균이 특정 값과 다른가
대응표본	같은 대상을 두 번 측정했다. 전후·처치
독립표본	다른 대상 두 집단을 비교한다

25회의 혈압약 문제가 대응표본, 26회의 남녀 혈압 문제가 독립표본이었다. 같은 사람을 두 번 쟀는지만 보면 갈린다.

alternative 는 단측 검정에 쓴다. “복용 후 혈압이 더 낮아졌을 것”처럼 방향이 있는 가설이면 'less' 나 'greater' 를 명시한다. 안 쓰면 양측이고, 그러면 p 값이 두 배가 된다.

독립표본 t 검정은 등분산을 전제한다. 아니면 equal_var=False (Welch) 를 준다. 그래서 순서가 이렇게 된다 — 정규성 확인 → 등분산 확인 → t 검정.

요약값만 주어졌을 때

24회는 데이터가 아니라 평균·표준편차·Z 임계값만 줬다. 이럴 땐 손으로 쓴다.

$Z = \frac{\bar{x}_1 - \bar{x}_2}{\sqrt{\dfrac{s_1^2}{n_1} + \dfrac{s_2^2}{n_2}}}$

import numpy as np
z = (5.7 - 5.6) / np.sqrt(0.03**2/n1 + 0.04**2/n2)
print(z, abs(z) > 1.65)   # 문제에서 준 Z(0.05) = 1.65 와 비교

정규성과 등분산

검정을 고르기 전에 하는 검정들이다. 둘 다 $p > 0.05$ 여야 “가정을 만족한다”가 된다. 귀무가설이 “정규분포를 따른다”, “분산이 같다”이므로 방향이 반대다.

from scipy import stats

stats.shapiro(data)             # 샤피로-윌크 — 표본이 작을 때(n < 50) 기본
stats.ks_2samp(x, y)            # 콜모고로프-스미르노프 — 두 분포가 같은가
stats.normaltest(data)          # 왜도·첨도 기반 (D'Agostino)

stats.levene(g1, g2, g3)        # 레빈 — 등분산. 정규성이 의심스러워도 견딘다
stats.bartlett(g1, g2, g3)      # 바틀렛 — 정규성을 전제한다

카이제곱 — 독립성 검정

교차표를 만들고 넣는다. 26회의 지역별 지지율, 23회의 학과별 성적 문제가 전부 이것이다.

import pandas as pd
from scipy.stats import chi2_contingency, fisher_exact

ct = pd.crosstab(df['gender'], df['smoking'])
chi2_stat, p_value, dof, expected = chi2_contingency(ct)

print(chi2_stat, p_value, dof)
print(expected)     # 독립일 때의 기댓값

expected 를 반드시 출력한다. 23회는 “학과와 성적이 독립일 때 기댓값을 구하시오”가 따로 5점이었다. 손으로 구하는 식도 같이 알아 둔다.

$E_{ij} = \frac{(\text{행 합}) \times (\text{열 합})}{\text{전체 합}}$

기대빈도가 5 미만인 칸이 있으면 카이제곱을 쓰면 안 된다. 그때는 피셔의 정확검정이다.

odds_ratio, p_value = fisher_exact(ct)     # 2×2 에만 쓴다

분산분석

세 집단 이상의 평균 비교. 두 집단이면 t 검정이다.

from scipy import stats
fo = stats.f_oneway(g1, g2, g3)       # 일원배치 — 빠른 길

표를 만들어 제출해야 하면 statsmodels 쪽이다. 21회는 “이원분산분석을 수행하고 통계표를 작성”이었다.

import statsmodels.api as sm
from statsmodels.formula.api import ols

# 일원배치
model = ols('value ~ item', data=DF).fit()
print(sm.stats.anova_lm(model))

# 이원배치 — item:layer 가 교호작용이다
model = ols('value ~ item + layer + item:layer', data=DF).fit()
print(sm.stats.anova_lm(model, type=2))

item + layer 만 쓰면 주효과만 본다. item:layer 를 넣어야 교호작용이 들어간다. 이원분산분석 문제에서 이걸 빼면 표의 행 하나가 통째로 빈다.

교호작용은 그림으로 확인한다. 두 선이 평행하면 교호작용이 없다.

from statsmodels.graphics.factorplots import interaction_plot

interaction_plot(DF['item'], DF['layer'], DF['value'],
                 colors=['pink', 'green'], markers=['D', '^'], ms=10)
plt.show()

잔차의 정규성은 QQ 플롯으로 본다.

sm.qqplot(model.resid, stats.t, fit=True, line='45')
plt.show()

사후검정

분산분석은 “어딘가 다르다”까지만 말한다. 어디가 다른지는 사후검정이다.

from statsmodels.stats.multicomp import pairwise_tukeyhsd, MultiComparison

# Tukey's HSD
tukey = pairwise_tukeyhsd(DF['value'], DF['item'])
print(tukey)

# Scheffé
schf = MultiComparison(DF['value'], DF['item']).allpairtest(stats.ttest_ind, method='s')
print(schf[0])

정리하다 발견한 것 하나 — 예전 메모에 Duncan 검정을 pairwise_tukeyhsd(값, 그룹, 'Duncan') 으로 적어 두었는데 틀렸다. 세 번째 인자는 alpha 라서 문자열을 넣으면 유의수준 자리에 들어간다. statsmodels 에 Duncan 은 없다. 사후검정을 여러 개 비교해야 하면 scikit-posthocs 를 쓴다.

사후검정	성격
Tukey HSD	모든 쌍을 비교. 집단 크기가 같을 때 표준
Scheffé	가장 보수적. 어떤 대비도 가능하지만 잘 안 나온다
Bonferroni	유의수준을 비교 횟수로 나눈다. 단순하고 안전

상관 분석

DF.cov()                                              # 공분산
DF['value'].corr(DF['item'], method='pearson')        # 피어슨 — 선형, 연속
DF['value'].corr(DF['item'], method='spearman')       # 스피어만 — 순위, 단조
DF['value'].corr(DF['item'], method='kendall')        # 켄달 — 순위, 표본이 작을 때

corr, p = stats.pearsonr(DF['value'], DF['item'])     # 계수와 유의확률을 같이
corr, p = stats.kendalltau(DF['value'], DF['item'])

상관계수만 쓰고 끝내면 안 된다. 0.4 라는 숫자 자체는 아무 말도 하지 않는다. pearsonr 로 p 값까지 내야 “유의한 상관이 있다”고 쓸 수 있다.

정규성이 깨졌거나 순위 자료면 스피어만이다. 시각화는 세 가지를 준비했다.

plt.scatter(DF['item'], DF['value'])          # 산점도
sns.pairplot(DF)                               # 산점도 행렬
sns.heatmap(DF.corr(), cmap='coolwarm', annot=True, vmin=-1, vmax=1)

히트맵의 vmin=-1, vmax=1 을 빼면 색 눈금이 데이터 범위에 맞춰 늘어나서, 0.3 짜리 상관이 새빨갛게 보인다. 읽는 사람을 속이는 그림이 된다.

비모수 검정

정규성이 깨졌거나 순위·부호만 있을 때. 23회·25회에 연달아 나왔다.

from scipy import stats
from statsmodels.stats.descriptivestats import sign_test

# 부호 검정 — 23회. 진공관 수명이 1만 시간이라는 주장
stat, p = sign_test(data, mu0=10000)

# 크루스칼-왈리스 — 25회. 세 공장의 생산량 (일원분산분석의 비모수 버전)
stats.kruskal(x, y, z)

# 윌콕슨 부호순위 — 대응표본 t 의 비모수 버전
stats.wilcoxon(before, after)

# 만-휘트니 U — 독립표본 t 의 비모수 버전
stats.mannwhitneyu(a, b)

부호 검정에서 “유효한 샘플의 수를 계산”이 따로 배점이었다. 기준값과 정확히 같은 관측치는 버린다. 12개 중 하나가 정확히 10000이면 유효 표본은 11이다.

22회의 aabbaaaabbbbab 문제는 런 검정(연의 수 검정)이다. 두 값이 무작위로 섞였는지를 본다.

from statsmodels.sandbox.stats.runs import runstest_1samp
z, p = runstest_1samp(binary_sequence)

모수 검정	비모수 대응
일표본 t	부호 검정, 윌콕슨 부호순위
대응표본 t	윌콕슨 부호순위
독립표본 t	만-휘트니 U
일원배치 분산분석	크루스칼-왈리스

잡다한 계산 문제들

25회는 통계 문항 안에 계산 문제가 섞여 있었다.

from scipy.stats import hmean, gmean

# 갈 때 4km/h, 올 때 5km/h 의 왕복 평균 속도 → 조화평균
hmean([4, 5])            # 4.444...

# 매출이 3000 → 4000 → 5000 일 때 연평균 증가 배율 → 기하평균
gmean([4000/3000, 5000/4000])    # 1.1547

속도(비율)의 평균은 조화평균, 성장률의 평균은 기하평균이다. 산술평균을 쓰면 둘 다 틀린다.

베이즈 정리 문제(24회 — 양성 판정을 받은 사람이 실제 양성일 확률)는 공식 그대로다.

$P(A|B) = \frac{P(B|A)\,P(A)}{P(B)}$

NPV 문제(25회)는 라이브러리를 찾지 말고 직접 쓴다. np.npv 는 numpy 1.20 에서 빠졌고 numpy-financial 은 시험 환경에 없다.

npv = sum(cf / (1 + r)**t for t, cf in enumerate(cash_flows))