스튜던트 t 분포¶

(1) 해석적으로. 비의 로그를 잡는다. \(t_d\)의 밀도가 \(f_{t_d}(x) = C_d\,(1 + x^2/d)^{-(d+1)/2}\)이므로

이다. \(R\)은 \(x^2\)만 통해 \(x\)에 의존하므로 \(u = x^2\)로 바꿔 미분한다.

분모가 양수이므로 부호는 \(u - 1\)이 정한다. 비는 \(u < 1\)에서 줄고 \(u > 1\)에서 늘며, 최소가 정확히 \(|x| = 1\)에서 잡힌다. 자유도가 어떤 값이어도 그 자리는 \(1\)이다.

출발점은 \(1\)보다 작다. \(u = 0\)에서

인데, 가운데 등식은 다음 절 연습문제 10에 나오는 \(E[\sqrt{\chi^2_d}] = \sqrt2\,\Gamma((d+1)/2)/\Gamma(d/2)\)이고, 부등식은 \(\sqrt{\cdot}\)이 오목하다는 옌센 부등식이다. \(t_d\)의 봉우리가 정규분포보다 낮은 몫이 곧 그 절의 편향상수 \(c_4\)이며, 같은 옌센 부등식이 두 자리에서 일하고 있는 셈이다.

끝점은 \(1\)보다 크다. \(x \to \infty\)에서 \(t\)의 꼬리는 \(|x|^{-(d+1)}\)로 다항식처럼 줄고 정규는 \(e^{-x^2/2}\)로 줄므로 \(R(u) \to +\infty\)다.

그러므로 \(R\)은 음수에서 출발해 \(u = 1\)까지 더 내려갔다가 돌아서 \(+\infty\)로 간다. 영점은 \(u > 1\)에 정확히 하나뿐이다. 그 자리를 \(x^* > 1\)이라 하면 \(x > 0\)에서 교차가 한 번이고, 두 밀도가 모두 대칭이므로

다. 대칭인 두 밀도가 \(0\)에서 다르면 교차 횟수는 짝수여야 한다는 것과도 맞는다.

여기서 미리 못박아 둘 것이 있다. "꼬리가 두껍다"는 말은 "어디서나 높다"는 뜻이 아니다. \(|x| < x^*\)에서는 \(t_d\)가 정규보다 오히려 낮다. 밀도의 총합이 양쪽 다 1이니 어딘가가 높으면 어딘가는 낮아야 하고, 그 낮은 구간이 바로 눈에 잘 띄는 가운데다.

(2) 수치적으로.

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

fig, (ax_full, ax_tail) = plt.subplots(1, 2, figsize=(12, 3))
x = np.linspace(-4, 4, 200)

# 왼쪽: 전체 모습. 두 곡선이 거의 겹쳐 보인다.
ax_full.plot(x, stats.norm().pdf(x), label='Normal')
ax_full.plot(x, stats.t(df=10).pdf(x), label='t(10)')
ax_full.set_title('Full PDF')
ax_full.legend()

# 오른쪽: 꼬리만 확대. x[-50:] 은 x 배열의 마지막 50개, 즉 오른쪽 끝이다.
# 전체 그림에서는 밀도가 너무 작아 안 보이던 차이가 여기서 드러난다.
# **꼬리는 언제나 확대해서 봐야 한다.** 검정에서 문제가 되는 곳이 바로 꼬리다.
ax_tail.plot(x[-50:], stats.norm().pdf(x[-50:]), label='Normal')
ax_tail.plot(x[-50:], stats.t(df=10).pdf(x[-50:]), label='t(10)')
ax_tail.set_title('Right Tail (zoomed)')
ax_tail.legend()

plt.tight_layout()
plt.show()

왼쪽에서는 두 곡선이 겹쳐 보이지만 오른쪽 확대에서는 \(t(10)\)이 정규분포 위로 또렷이 올라앉아 있다. 자유도가 10이나 되는데도 그렇다.

(1)이 유도한 세 가지, 곧 \(|x| = 1\)의 최소점과 \(0\)에서의 비와 교차점을 하나씩 확인한다.

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

d = 10
ratio = lambda x: stats.t(d).pdf(x) / stats.norm.pdf(x)

# 0 에서의 비가 c4 상수와 같다는 등식을 먼저 확인한다.
c4 = np.sqrt(2 / d) * special.gamma((d + 1) / 2) / special.gamma(d / 2)
print(f"f_t(0)/phi(0) = {ratio(0):.9f},   E[sqrt(chi2_d)]/sqrt(d) = {c4:.9f}")

# 비가 |x| = 1 에서 최소라는 유도를 격자로 확인한다.
xs = np.linspace(0, 3, 30_001)
print(f"비가 최소인 곳 = {xs[ratio(xs).argmin()]:.4f},  그때의 비 = {ratio(xs).min():.4f}")

# 두 밀도가 만나는 곳. 유도에 따르면 |x| > 1 에 하나뿐이다.
x_star = optimize.brentq(lambda x: ratio(x) - 1, 1.01, 10)
print(f"교차점 x* = {x_star:.4f}  (그림 오른쪽 칸은 {np.linspace(-4, 4, 200)[-50]:.4f} 부터 시작한다)")

print(f"{'x':>8}{'t(10) 밀도':>13}{'정규 밀도':>13}{'비':>9}")
for x in (0.0, 1.0, 1.5928, 2.0302, 3.0, 4.0):
    print(f"{x:>8.4f}{stats.t(d).pdf(x):>13.6f}{stats.norm.pdf(x):>13.3e}{ratio(x):>9.3f}")

# 꼬리확률로 옮겨 보면
for x in (2.0302, 4.0):
    print(f"P(X > {x:.4f}):  t(10) {stats.t(d).sf(x):.3e}   정규 {stats.norm.sf(x):.3e}"
          f"   비 {stats.t(d).sf(x)/stats.norm.sf(x):.2f}")

출력:

f_t(0)/phi(0) = 0.975350077,   E[sqrt(chi2_d)]/sqrt(d) = 0.975350077
비가 최소인 곳 = 1.0000,  그때의 비 = 0.9520
교차점 x* = 1.5928  (그림 오른쪽 칸은 2.0302 부터 시작한다)
       x     t(10) 밀도        정규 밀도        비
  0.0000     0.389108    3.989e-01    0.975
  1.0000     0.230362    2.420e-01    0.952
  1.5928     0.112203    1.122e-01    1.000
  2.0302     0.058303    5.080e-02    1.148
  3.0000     0.011401    4.432e-03    2.572
  4.0000     0.002031    1.338e-04   15.176
P(X > 2.0302):  t(10) 3.489e-02   정규 2.117e-02   비 1.65
P(X > 4.0000):  t(10) 1.259e-03   정규 3.167e-05   비 39.76

셋 모두 유도와 맞는다. \(0\)에서의 비는 아홉째 자리까지 \(c_4\)와 같고, 격자가 찾은 최소점은 \(1.0000\)이며, 교차점은 \(x^* = 1.5928\) 하나다.

그림의 두 칸이 서로 다른 구간을 보여 준다. 오른쪽 칸은 x[-50:], 곧 \(x \ge 2.0302\)만 그리는데 이는 교차점 \(1.5928\)보다 바깥이므로 전부 \(t\)가 더 높은 구간이다. 거기서 비가 \(1.148\)에서 \(15.176\)까지 오르니 눈에 보이는 것이 당연하다. 반대로 왼쪽 칸에서 눈길이 가는 봉우리 근처는 \(|x| < 1.5928\), 곧 \(t\)가 더 낮은 구간이다. 그 구간에서 비는 \(0.952\)와 \(1\) 사이에만 머문다. 두 곡선이 "겹쳐 보인다"는 인상은 차이가 \(5\%\)를 넘지 않는다는 사실의 반영이다.

왼쪽 칸이 가리는 것이 둘 있다. 하나는 교차점이다. \(x^* = 1.5928\)에서 두 밀도가 정확히 같으므로 그 자리에는 볼 것이 아무것도 없고, 그래서 선이 교차하는 모습 자체가 보이지 않는다. 다른 하나는 꼬리의 차이다. \(x = 4\)에서 비가 \(15\)배인데도 밀도가 \(0.002\)와 \(0.00013\)이라 세로축에서 한 화소도 차지하지 못한다. 배율이 가장 큰 곳이 그림에서 가장 작게 보인다는 것이 선형 세로축의 성질이다.

검정에 쓰이는 양으로 옮기면 차이가 더 또렷해진다. 밀도가 아니라 꼬리확률로 재면 \(x = 2.0302\)에서 \(1.65\)배, \(x = 4\)에서 \(39.76\)배다. 임계값과 \(p\)값은 모두 이쪽 눈금에서 결정되며, 그래서 자유도 \(10\)이라는 "꽤 큰" 자유도에서도 \(t\)를 써야 한다.

(1) 이론값. 두 비가 따르는 분포는 둘 다 정확히 알려져 있다. 분모에 \(\sigma\)를 넣으면

이고, \(S\)로 갈아 끼우면 본문의 유도대로 \(T \sim t_{n-1} = t_9\)다. 4.2절의 \(\operatorname{Var}(t_d) = d/(d-2)\)에 \(d = 9\)를 넣으면

이다. 늘어난 몫 \(d/(d-2) - 1 = 2/(d-2) = 2/7 = 0.2857\)이 \(\sigma\)를 모른 채 추정하느라 치른 값이며, 이 식이 \(d > 2\)에서만 뜻을 갖는다는 것도 함께 기억해 둘 일이다(연습문제 8). \(d = 9\)는 넉넉히 넘으므로 분산이 유한하다.

꼬리확률은 같은 임계값에서 재어야 비교가 된다.

곧 \(1.96\)으로 양측 \(5\%\) 검정을 설계했다고 믿었는데 실제로는 \(8.2\%\)에서 기각한다. 명목의 \(1.63\)배다. \(t_9\)에서 양쪽 \(5\%\)를 지키는 올바른 임계값은 \(t_{9,\,0.975} = 2.2622\)로 \(1.96\)보다 \(15.4\%\) 크다(연습문제 9의 표에 같은 수치가 있다).

(2) 모의실험. 아래 코드는 쪽에 있던 것에 비교 줄만 덧붙인 것이다. \(Z\)는 이미 뽑아 둔 samples에서 분모만 바꿔 계산하므로 난수를 새로 쓰지 않고, 따라서 그림도 달라지지 않는다.

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

np.random.seed(0)
n, mu, sigma = 10, 0, 10
n_sim = 10_000

# (n, n_sim) 배열이므로 **열 하나가 표본 하나**다. axis=0 으로 집계한다.
samples = np.random.normal(mu, sigma, (n, n_sim))
x_bar = samples.mean(axis=0)
s = samples.std(axis=0, ddof=1)

# t 통계량. 분모에 참 sigma(=10)가 아니라 **표본표준편차 s** 를 넣는 것이 요점이다.
# sigma를 썼다면 이 값은 정확히 N(0,1)을 따랐을 것이다.
# s를 쓰면 분모 자체가 흔들리므로 분포가 넓어지고 꼬리가 두꺼워진다.
# 그 넓어진 정도를 정확히 기술하는 것이 t(n-1) 이다.
t_stats = (x_bar - mu) / (s / np.sqrt(n))

# **같은 표본**에서 분모만 참 sigma 로 되돌린 것. 비교 대상이다.
# 난수를 새로 뽑지 않으므로 아래 그림은 달라지지 않는다.
z_stats = (x_bar - mu) / (sigma / np.sqrt(n))

sd_t = np.sqrt((n - 1) / (n - 3))              # Var(t_d) = d/(d-2)
p_t = 2 * stats.t(n - 1).sf(1.959964)
print(f"표준편차     : z 이론 1.0000 모의 {z_stats.std(ddof=1):.4f}  |  t 이론 {sd_t:.4f} 모의 {t_stats.std(ddof=1):.4f}")
print(f"P(|.| > 1.96): z 이론 0.0500 모의 {np.mean(np.abs(z_stats) > 1.959964):.4f}  |  t 이론 {p_t:.4f} 모의 {np.mean(np.abs(t_stats) > 1.959964):.4f}")
print(f"KS 검정 p : t_9 과 맞는가 {stats.kstest(t_stats, 't', args=(n-1,)).pvalue:.4f}"
      f"   N(0,1) 과 맞는가 {stats.kstest(t_stats, 'norm').pvalue:.2e}")
print(f"올바른 임계값 t(9, 0.975) = {stats.t(n-1).ppf(0.975):.4f}")

fig, ax = plt.subplots(figsize=(12, 3))
bins = np.arange(-6, 6, 0.1)
ax.hist(t_stats, bins=bins, density=True, alpha=0.7, label=f'Simulated $t_{{{n-1}}}$')
ax.plot(bins, stats.t(n-1).pdf(bins), '--r', lw=2, label=f'$t_{{{n-1}}}$ PDF')
ax.legend()
ax.spines[['top', 'right']].set_visible(False)
plt.show()

출력:

표준편차     : z 이론 1.0000 모의 0.9904  |  t 이론 1.1339 모의 1.1300
P(|.| > 1.96): z 이론 0.0500 모의 0.0459  |  t 이론 0.0816 모의 0.0826
KS 검정 p : t_9 과 맞는가 0.8758   N(0,1) 과 맞는가 6.74e-04
올바른 임계값 t(9, 0.975) = 2.2622

네 수치가 모두 이론과 맞는다. \(Z\) 쪽은 \(0.9904\)와 \(0.0459\)로 \(1\)과 \(0.05\)에 붙고(되풀이 1만 번에서 표준편차의 몬테카를로 오차가 \(0.0071\), 비율의 오차가 \(0.0022\)이니 각각 \(1.4\)배와 \(1.9\)배 안이다), \(T\) 쪽은 \(1.1300\)과 \(0.0826\)으로 이론값 \(1.1339\)와 \(0.0816\)에 붙는다. 히스토그램에 겹친 빨간 곡선이 \(t_9\) 밀도인데, KS 검정이 \(t_9\)에 대해 \(p = 0.88\)을 주고 \(N(0,1)\)에 대해 \(p = 0.00067\)을 주어 그림으로는 닮아 보이는 둘을 가른다.

잃는 것은 유의수준이다. 분모 한 글자를 \(\sigma\)에서 \(S\)로 바꾸었을 뿐인데 표준편차가 \(13\%\) 넓어졌고, 그 넓어진 분포에 \(1.96\)을 그대로 대면 기각률이 \(5\%\)에서 \(8.2\%\)로 올라간다. 신뢰구간 쪽으로 옮기면 명목 \(95\%\) 구간이 실제로는 \(91.8\%\)만 담는다는 뜻이고, 연습문제 2의 표에서 \(n = 10\) 줄이 바로 이 수치다.

그런데 이것은 손실이 아니라 청구서다. 분포가 \(t_9\)로 바뀐 것을 알고 임계값을 \(2.2622\)로 고치면 명목 \(5\%\)가 정확히 지켜진다. 모집단이 정규이기만 하면 근사가 아니라 등식이기 때문이다. \(\sigma\)를 모른다는 문제는 이 자리에서 값이 매겨지고 끝난다. 뒤의 두 절에서 볼 \(\chi^2\)과 \(F\)에서는 이렇게 끝나지 않는다.

\(\mathrm{SE} = s/\sqrt n = 8000/3 \approx 2667\). 검정통계량: \(t = (45000 - 40000)/2667 \approx 1.875\).

\(t_8\) 아래에서 \(P(T \ge 1.875) \approx 0.048\)로 약 4.8%이다.

참고: \(\sigma\)를 모르고 표본에서 \(s\)를 추정하여 추가적인 불확실성이 들어오므로 \(z\) 대신 \(t\)를 사용한다.

\(\sigma\)를 모르면 \((\bar X - \mu)/(S/\sqrt n)\)의 정확한 분포는 \(t_{n-1}\)이다. 따라서 \(z = 1.96\)으로 만든 구간이 실제로 담는 확률은 \(P(|T_{n-1}| \le 1.96)\)이며, 이것은 95%가 아니다.

읽는 법. \(n = 5\)에서 95%라고 이름 붙인 구간이 실제로는 88%만 담는다. 7%포인트의 차이는 작지 않다. \(n = 30\)에서도 94.0%로, 검정 쪽으로 옮겨 말하면 명목 5% 양측검정의 실제 유의수준이 5.97%가 된다. 경험 법칙이 말하는 "거의 같다"는 이 정도를 뜻한다.

결론. 이 오차는 크지 않지만 치를 이유가 없는 오차다. 올바른 임계값을 얻는 데 드는 수고가 0이기 때문이다. \(\sigma\)를 모르면 표본크기와 무관하게 \(t\)를 쓰면 되고, \(n\)이 커지면 임계값이 저절로 1.96으로 수렴하므로 큰 표본에서 손해도 없다. 30이라는 문턱은 계산기가 없던 시절 분포표를 뒤지는 수고를 아끼려고 생긴 것이며, 지금은 그 수고가 사라졌다.

한 가지 주의할 점이 있다. 여기서 평가한 것은 \(t\) 대 \(z\)의 선택뿐이다. 같은 숫자 30이 "\(n \ge 30\)이면 중심극한정리가 듣는다"는 별개의 주장에도 쓰이는데, 그것은 위 표와 아무 상관이 없다. 위 계산은 모집단이 정규라고 가정하고 한 것이기 때문이다.

귀무가설 \(\mu = \mu_0\) 아래에서 \(\bar X \sim N(\mu_0, \sigma^2/n)\)이므로 \(Z = (\bar X - \mu_0)/(\sigma/\sqrt n) \sim N(0, 1)\)이다.

표본분산 \(s^2\)을 \(\sigma^2\)으로 척도조정하면 카이제곱 분포를 따른다: (정규 자료에 대해) \((n-1)s^2/\sigma^2 \sim \chi^2_{n-1}\).

정규 자료에서 \(\bar X\)와 \(s^2\)은 독립이다(정규성에만 특유한 자명하지 않은 사실이다).

따라서:

이며 \(V \sim \chi^2_{n-1}\)은 \(Z\)와 독립이다. \(t\)의 정의에 의해 이는 \(t_{n-1}\)이다.

\(\sigma\) 대신 \(s\)를 쓰면 분모에 카이제곱이 들어온다. \(t\) 분포는 정규분포보다 두꺼운 꼬리로 이 추가 불확실성을 반영한다.

\(t_3\): \(P(|T| > 2) = 2 \cdot P(T > 2)\). \(t_3\) 분포표에서 \(P(T > 2) \approx 0.07\)이므로 \(P(|T| > 2) \approx 0.14\)이다.

\(P(|T| > 4) \approx 2 \cdot 0.014 = 0.028\).

정규분포: \(P(|Z| > 2) \approx 0.046\), \(P(|Z| > 4) \approx 6 \times 10^{-5}\).

\(|x| = 4\)에서 비교하면 정규분포의 확률은 \(6 \times 10^{-5}\)로 사실상 0인데 \(t_3\)의 확률은 \(0.028\)로 400배가 넘는다. \(t_3\)은 정규분포보다 극단적인 결과에 훨씬 많은 확률을 부여한다.

실무적 귀결: 같은 유의수준에서 소표본 \(t\) 검정은 \(z\) 검정보다 검정력이 낮다. 두꺼운 꼬리를 상쇄하기 위해 임계값이 더 크기 때문이다. 이는 가정(\(\sigma\)를 안다는 것)과 꼬리에 대한 보수성 사이의 맞바꿈이다.

왜도에 훨씬 민감하다.

왜 그런가. \(t = \sqrt n(\bar X-\mu)/S\)에서 분자와 분모가 독립이 아니라는 것이 열쇠다. 정규모집단에서는 \(\bar X\)와 \(S^2\)이 독립이지만, 모집단이 치우쳐 있으면 그렇지 않다. 사실 이 독립성은 정규분포만의 성질이다. 둘이 독립이면 모집단이 정규라는 역이 성립한다는 것이 루카치의 정리(Lukacs 1942)다. 게리(Geary 1936)가 먼저 관련된 결과를 얻었으나, 이 특징화를 가리킬 때는 루카치로 부르는 것이 표준이다.

구체적으로 \(\operatorname{Cov}(\bar X, S^2) = \mu_3/n\)이다. 모집단의 3차 중심적률이 곧 공분산을 만든다. 오른쪽으로 치우친 모집단(\(\mu_3>0\))이면 \(\bar X\)가 우연히 클 때 \(S\)도 함께 크게 나온다. 그러면 분자가 커져도 분모가 함께 커져 \(t\)가 눌리고, 반대로 \(\bar X\)가 작을 때는 \(S\)도 작아 \(t\)가 과장된다.

그 결과 \(t\) 통계량의 분포 자체가 왼쪽으로 치우친다(모집단이 오른쪽으로 치우쳤을 때). 에지워스 전개로 보면 오차의 주항이

로 왜도 \(\gamma_1\)에 비례하고 \(n^{-1/2}\) 차수다.

첨도의 영향. 첨도는 \(n^{-1}\) 차수의 항에만 들어간다. 수렴이 한 차수 빠르므로 같은 \(n\)에서 영향이 훨씬 작다.

이 구조에서 세 가지가 따라 나온다. 첫째, 한쪽 꼬리가 양쪽보다 위험하다. 왜도가 만드는 오차는 한쪽 꼬리를 부풀리고 다른 쪽을 줄이므로 양쪽을 함께 보면 두 오차가 상당 부분 상쇄되지만, 한쪽만 보면 그대로 남는다. 다음 문제에서 그 비대칭을 실제로 재어 본다. 둘째, 그래서 \(t\)를 근거로 만든 구간도 좌우가 같은 폭이어서는 안 되는 상황이 생기며, 이때 쓰는 왜도 보정이나 부트스트랩 계열의 구간은 16장과 8장의 몫이다. 셋째, 꼬리가 두껍기만 한 대칭 모집단(왜도가 0이고 첨도만 큰 경우, 예를 들어 \(t_3\))에서는 \(t\) 통계량의 표본분포가 거의 어긋나지 않는다. 연습문제 7이 이것을 확인한다.

한 문장으로 줄이면 치우침은 표본분포의 중심과 꼬리를 비틀고, 두꺼운 꼬리는 그러지 않는다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
B = 200_000

print(f"{'n':>6}{'양쪽':>9}{'오른쪽':>9}{'왼쪽':>9}")
for n in (10, 30, 100, 1000):
    c2 = stats.t(n - 1).ppf(0.975)      # 양쪽 5%
    c1 = stats.t(n - 1).ppf(0.95)       # 한쪽 5%
    X = rng.exponential(size=(B, n))    # Exp(1): mu = sigma = 1
    T = (X.mean(axis=1) - 1) / (X.std(axis=1, ddof=1) / np.sqrt(n))
    print(f"{n:>6}{np.mean(np.abs(T) > c2):>9.3f}"
          f"{np.mean(T > c1):>9.3f}{np.mean(T < -c1):>9.3f}")

출력:

     n       양쪽      오른쪽       왼쪽
    10    0.099    0.014    0.132
    30    0.073    0.022    0.098
   100    0.057    0.032    0.075
  1000    0.051    0.043    0.058

줄어든다. 명목 5% 양쪽이 \(n = 10\)에서 0.099였다가 \(n = 1000\)에서 0.051로 내려온다. 한쪽 꼬리도 0.014와 0.132로 크게 갈려 있던 것이 0.043과 0.058로 모인다.

한쪽 꼬리가 훨씬 더 어긋난다는 점이 이 표의 요점이다. \(n = 10\)에서 양쪽은 0.099로 명목의 두 배지만 왼쪽 하나만 보면 0.132로 거의 세 배다. 오른쪽이 0.014로 명목에 크게 못 미치는 덕분에 둘을 합칠 때 오차가 일부 상쇄되기 때문이다. 연습문제 5가 말한 "한쪽 꼬리가 양쪽보다 위험하다"가 이 숫자다. 어긋나는 방향이 왼쪽인 것도 예고와 맞는다. 지수모집단은 오른쪽으로 치우쳐 \(\mu_3 > 0\)이므로 \(t\)의 분포가 왼쪽으로 끌린다.

\(S^2\)과 무엇이 다른가. 줄어든다는 것 자체가 다르다. 다음 절 연습문제 6에서 같은 지수모집단·같은 명목 5%로 \((n-1)S^2/\sigma^2\)을 쟀을 때는 \(n = 10\)의 0.233이 \(n = 1000\)에서 0.324로 올라간다. 여기서 고장 나는 것은 정리의 수렴이라 표본이 약이 되고, 거기서 고장 나는 것은 정리의 전제라 표본이 약이 되지 못한다. 같은 자료, 같은 모집단, 반대 방향이다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
B = 200_000

print(f"{'n':>6}{'양쪽':>9}{'오른쪽':>9}{'왼쪽':>9}")
for n in (10, 30, 100):
    c2 = stats.t(n - 1).ppf(0.975)
    c1 = stats.t(n - 1).ppf(0.95)
    X = stats.t(5).rvs(size=(B, n), random_state=rng)   # 대칭, beta2 = 9
    T = X.mean(axis=1) / (X.std(axis=1, ddof=1) / np.sqrt(n))
    print(f"{n:>6}{np.mean(np.abs(T) > c2):>9.3f}"
          f"{np.mean(T > c1):>9.3f}{np.mean(T < -c1):>9.3f}")

출력:

     n       양쪽      오른쪽       왼쪽
    10    0.045    0.048    0.048
    30    0.047    0.049    0.049
   100    0.049    0.050    0.049

\(n = 10\)에서 이미 맞는다. 첫 줄의 세 수치가 0.045, 0.048, 0.048로 명목값 언저리에 있고 좌우가 갈리지도 않는다. 첨도가 9로 지수분포와 똑같은데도 앞 문제의 0.099, 0.014, 0.132와는 딴판이다. 어긋남을 만든 것은 첨도가 아니라 왜도였다는 뜻이며, 에지워스 전개의 주항이 \(\gamma_1/\sqrt n\)이고 첨도가 \(n^{-1}\) 항에만 들어간다는 연습문제 5의 계산이 이렇게 나타난다.

오류율이 지켜진다고 안심할 일은 아니다. \(t_5\) 모집단에서 이따금 나오는 극단값이 \(S\)를 부풀리므로 같은 명목 수준에서 진짜 차이를 잡아내는 힘은 떨어진다. 표본분포의 어긋남은 없고 통계량의 효율이 나쁜 것이며, 두 가지는 별개의 문제다. 그 효율 비교는 순위 방법을 다루는 16장의 몫이다.

이 두 문제를 나란히 놓으면 이 쪽의 결론이 정리된다. \(t\)의 표본분포는 정규성 가운데 대칭성에 기대고 있으며, 그 기댐조차 표본을 키우면 느슨해진다.

평균. \(t_d\)의 밀도는 \(0\)에 대해 대칭이므로, 적분이 수렴하기만 하면 평균은 \(0\)이다. 꼬리가 \(|t|^{-(d+1)}\)로 줄므로 \(\int |t|\,f(t)\,dt\)가 수렴하려면 \(d > 1\)이어야 한다.

분산. \(T = Z/\sqrt{V/d}\)에서 \(Z \perp V\)이므로

이다. \(E[d/V]\)는 역카이제곱의 기댓값이며 \(d>2\)에서만 유한하다. 따라서

이다. \(\square\)

import numpy as np
from scipy import stats

print(f"{'d':>5}{'Var 이론':>12}{'모의':>14}")
for d in (1, 2, 3, 5, 10, 30):
    X = stats.t.rvs(d, size=500_000, random_state=2)
    th = 'inf' if d <= 2 else f"{d/(d-2):.4f}"
    print(f"{d:>5}{th:>12}{X.var():>14.4f}")

출력:

    d      Var 이론            모의
    1         inf   317384.7758
    2         inf       23.3971
    3      3.0000        2.9669
    5      1.6667        1.6499
   10      1.2500        1.2435
   30      1.0714        1.0644

\(d = 1\)의 표본분산이 \(317{,}384\)로 터무니없다. 무한한 양을 유한 표본으로 재려 했으니 당연하다. 난수 씨앗을 바꾸면 이 값이 완전히 달라지는데, 수렴하지 않는다는 것이 바로 그런 모습이다. \(d = 2\)의 \(23.4\)도 마찬가지다.

\(d = 1\)이면 코시분포다. \(t_1\)의 밀도는

로 표준 코시분포와 같다. 평균조차 존재하지 않으므로 큰수의 법칙이 성립하지 않는다. 코시 표본의 평균은 아무리 표본을 키워도 한 점으로 수렴하지 않고, 다시 코시분포를 따른다(7장에서 다룬다).

\(d \le 2\)에서 분산이 무한한 것은 \(t\)의 정의에서 곧바로 읽힌다. 분모 \(\sqrt{V/d}\)가 \(0\) 근처의 값을 충분히 자주 가지면 \(T\)가 폭발하는데, \(\chi^2_d\)는 \(d\)가 작을수록 \(0\) 근처에 질량이 많다. \(d=1\)이면 밀도가 \(0\)에서 발산하기까지 한다.

실무적 함의. \(t\) 분포를 잡음 모형으로 쓸 때(강건 회귀, 베이즈 모형) 자유도를 \(2\) 이하로 두면 분산이 존재하지 않아 평균제곱오차 같은 기준이 무의미해진다. 보통 \(d \ge 3\), 흔히 \(4\sim7\)을 쓰는 이유이며, 꼬리는 두껍게 두되 적률은 남겨 두려는 절충이다.

from scipy import stats

print(f"{'n':>6}{'t_{n-1,.975}':>14}{'z 대비':>12}")
for n in (3, 5, 10, 20, 30, 60, 120):
    tq = stats.t.ppf(0.975, n - 1)
    print(f"{n:>6}{tq:>14.4f}{(tq/1.959964 - 1) * 100:>11.2f}%")

출력:

     n  t_{n-1,.975}        z 대비
     3        4.3027     119.53%
     5        2.7764      41.66%
    10        2.2622      15.42%
    20        2.0930       6.79%
    30        2.0452       4.35%
    60        2.0010       2.09%
   120        1.9801       1.03%

\(n = 3\)이면 구간이 두 배가 넘는다. 자유도 \(2\)에서 임계값이 \(4.30\)이라, \(z\)를 쓴 구간보다 \(120\%\) 넓다. \(\sigma\)를 모른다는 사실 하나가 이만한 대가를 요구한다.

\(n = 30\)에서 \(4.35\%\). 흔히 말하는 "\(n \ge 30\)이면 \(z\)를 써도 된다"는 관행이 여기서 나온다. \(4\%\) 차이는 실무에서 대개 무시할 만하다. 다만 연습문제 2에서 본 대로 포함확률로 보면 이야기가 조금 다르므로, 굳이 \(z\)를 쓸 이유가 없다는 점도 함께 기억하라. 계산 비용이 같으니 \(t\)를 쓰면 된다.

감소 속도가 빠르다는 점이 요점이다.

\(t_{d} \to N(0,1)\)의 수렴은 \(d\)에 대해 \(O(1/d)\) 수준이라, 자유도가 조금만 늘어도 차이가 빠르게 줄어든다. 문제가 되는 구간은 \(n < 10\)이며, 거기서는 \(z\)와 \(t\)의 차이가 결론을 바꿀 만큼 크다.

작은 표본에서 진짜 문제는 따로 있다. \(n\)이 \(5\) 이하면 \(t\)를 제대로 쓰더라도 정규성 가정 자체를 확인할 방법이 없다. 2장 연습문제 10에서 본 대로 \(n=5\)의 표본왜도는 \(1.5\)를 넘을 수 없어 치우침을 탐지하지도 못한다. 임계값을 정확히 고르는 것보다 가정이 맞는지가 더 큰 문제다.

조건부로 보면 정규다. \(Z \sim N(0,1)\)과 \(V \sim \chi^2_d\)가 독립일 때 \(V = v\)로 고정하면

이다. 따라서 \(t_d\)는 분산이 \(d/V\)로 무작위한 정규분포들의 혼합이다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
d = 5
Z = rng.normal(0, 1, 500_000)
V = rng.chisquare(d, 500_000)
T = Z / np.sqrt(V / d)

print(f"  Z/sqrt(V/d) 가 t_{d} 인가?  KS p = "
      f"{stats.kstest(T, 't', args=(d,)).pvalue:.4f}")
print(f"  분산 모의 {T.var():.4f}   이론 {d/(d-2):.4f}")

출력:

  Z/sqrt(V/d) 가 t_5 인가?  KS p = 0.1685
  분산 모의 1.6672   이론 1.6667

꼬리가 두꺼워지는 메커니즘이 이 표현에 다 들어 있다. 분산이 하나로 고정된 정규분포는 꼬리가 \(e^{-t^2/2}\)로 아주 빠르게 죽는다. 그런데 분산 자체가 무작위라면, 가끔 \(V\)가 작게 나와 분산이 매우 큰 정규분포에서 값을 뽑게 된다. 그런 회차의 관측이 멀리 나가고, 그것들이 모여 두꺼운 꼬리를 만든다.

"이질적인 분산의 혼합이 꼬리를 두껍게 한다" 는 것이 일반 원리다. 4장 음이항분포 연습문제 2에서 포아송의 비율을 감마로 섞어 과대산포를 얻은 것과 정확히 같은 구조다.

세 번째 줄이 대비를 분명히 한다. 평균을 섞으면 정규가 유지되지만(앞 절 정규분포 쪽 연습문제 7의 선형결합), 분산을 섞으면 정규를 벗어난다.

어디에 쓰이는가.

• 강건 회귀. 잔차에 정규 대신 \(t\)를 가정하면 이상치에 덜 휘둘린다. 척도혼합 표현 덕분에 이것을 "관측마다 가중치가 다른 정규 회귀"로 바꿀 수 있고, 그 가중치를 반복 갱신하는 것이 EM 알고리즘이 된다.
• 금융 수익률 모형. 변동성이 시기마다 달라지는 자료를 "정규인데 분산이 변한다"로 보면 주변분포가 두꺼운 꼬리를 갖는다. 14장에서 수익률의 비정규성을 다룰 때 이 관점이 핵심이 된다.
• 베이즈 추론. 분산에 역감마 사전분포를 두고 적분해 없애면 평균의 사후분포가 정확히 \(t\)가 된다. \(t\) 분포가 베이즈와 빈도주의 양쪽에서 같은 자리에 나타나는 이유다. \(\square\)