무작위성에 대한 런 검정¶

1단계. 원소를 센다. \(n_1 = 8\)(양수), \(n_2 = 7\)(음수), \(N = 15\).

2단계. 런을 센다. \(+\!+\!+\), \(-\!-\), \(+\!+\), \(-\), \(+\), \(-\!-\!-\), \(+\!+\), \(-\)이므로 \(R = 8\).

3단계. 귀무가설 아래 평균과 표준편차를 계산한다.

4단계. \(Z\) 통계량을 계산한다.

5단계. \(p\)값을 계산한다. \(p = 2\,\Phi(-0.251) \approx 0.802\).

\(p = 0.802 \gg 0.05\)이므로 \(H_0\)을 기각하지 못한다. 자료는 무작위 수열과 일관된다.

(a) 런(같은 원소가 이어지는 최대 덩어리)을 찾으면

• 런 1: HH
• 런 2: TTT
• 런 3: HHHH
• 런 4: TT
• 런 5: H
• 런 6: T
• 런 7: HH
• 런 8: TTT
• 런 9: HH

런의 개수는 \(R = 9\)이다.

(b) \(n_H = 11\), \(n_T = 9\)이므로

(c)

양측 \(p\)값은

\(p = 0.377 \gg 0.05\)이므로 \(H_0\)을 기각하지 못한다. 이 수열은 비무작위성의 통계적으로 유의한 증거를 보이지 않는다. 관측된 런 개수 \(9\)가 기댓값 \(10.9\)에 가까워 무작위 수열과 일관된다.

두 수열 모두 \(n_1 = n_2 = 10\), \(N = 20\)이므로

import numpy as np
from scipy import stats

def runs_z(s):
    s = np.asarray(s); N = len(s)
    R = 1 + np.sum(s[1:] != s[:-1])
    n1 = s.sum(); n2 = N - n1
    mu = 2 * n1 * n2 / N + 1
    sd = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    z = (R - mu) / sd
    return R, z, 2 * stats.norm.sf(abs(z))

print(runs_z([1, 0] * 10))       # (20, 4.135, 3.55e-05)
print(runs_z([1, 1, 0, 0] * 5))  # (10, -0.459, 0.646)

출력:

(20, 4.135214625627066, 3.546230491468838e-05)
(10, -0.45946829173634074, 0.6458979170265631)

\(S_1\)은 잡아낸다. 완전한 교대는 런을 최대치인 \(20\)까지 밀어 올려 \(Z = +4.14\)가 된다.

\(S_2\)는 전혀 잡아내지 못한다. 주기가 4인 완벽하게 결정론적인 수열인데도 \(p = 0.646\)으로 무작위성과 완전히 일관되어 보인다. 런이 정확히 10개인데 이는 기댓값 11에 거의 딱 맞기 때문이다.

교훈: 런 검정은 지연 1의 자기상관에만 민감하다. 뭉침(지연 1 양의 상관)과 교대(지연 1 음의 상관)를 본다. 더 긴 주기의 구조는 런 개수를 정상 범위에 남겨 두므로 보이지 않는다. 주기적 패턴을 찾으려면 자기상관함수 전체나 주기도(periodogram)를 봐야 한다.

\(n_1 + n_2 = N\)을 고정하고 \(\mu_R = 2n_1n_2/N + 1\)을 보면, \(n_1 n_2\)는 \(n_1 = n_2 = N/2\)에서 최대가 된다. 즉 중앙값 부호화는 런 개수의 기댓값을 최대로 만든다(\(\mu_R \approx N/2 + 1\)).

이것이 좋은 이유는 검정의 해상도 때문이다. \(\mu_R\)이 크면 관측된 \(R\)이 위아래 양쪽으로 움직일 여지가 넓어져 뭉침과 교대를 모두 탐지할 수 있다. 극단적으로 \(n_1 = 1\)이면 \(\mu_R = 2(N-1)/N + 1 \approx 3\)이고 \(R\)이 가질 수 있는 값은 \(2\) 또는 \(3\)뿐이라 검정이 사실상 무력해진다.

평균을 쓰면 자료가 치우쳤을 때 \(n_1\)과 \(n_2\)가 크게 불균형해진다. 예를 들어 지수분포 자료에서는 관측값의 약 63%가 평균보다 작다.

import numpy as np
rng = np.random.default_rng(3)
x = rng.exponential(1, 2000)
print((x < x.mean()).mean())      # 0.6275  (이론값 1 - e^-1 = 0.632)
print((x < np.median(x)).mean())  # 0.5

출력:

0.6275
0.5

\(N = 2000\)에서 평균 부호화는 \(n_1 = 745\), \(n_2 = 1255\)를 주어 \(\mu_R = 936.0\)인 반면, 중앙값 부호화는 \(\mu_R = 1001.0\)이다. 손실이 크지는 않지만 치우침이 심해질수록 나빠진다.

더 중요한 문제는 해석이다. 중앙값 부호화에서 \(H_0\)은 "수열이 무작위이다"이고, 이는 자료의 분포와 무관하게 성립한다. 평균 부호화에서는 \(+\)의 비율 자체가 분포 모양에 의존하므로 검정이 순수하게 분포무관이 아니게 된다. 실무에서는 언제나 중앙값을 쓴다.

import numpy as np
rng = np.random.default_rng(5)
n = 40
t = np.arange(n)
y = 2 + 0.5 * t + rng.normal(0, 3, n)

def runs_z(res):
    s = (res > np.median(res)).astype(int)
    R = 1 + np.sum(s[1:] != s[:-1])
    n1 = s.sum(); n2 = n - n1; N = n
    mu = 2 * n1 * n2 / N + 1
    sd = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    return R, (R - mu) / sd

print("상수만:", runs_z(y - y.mean()))
b1, b0 = np.polyfit(t, y, 1)
print("추세 포함:", runs_z(y - (b0 + b1 * t)))

출력:

상수만: (8, -4.164701290861719)
추세 포함: (23, 0.6407232755171874)

상수만 적합하면 잔차의 부호가 앞쪽에 음수, 뒤쪽에 양수로 몰려 런이 \(8\)개에 그친다(기댓값 \(21\)). 잡음의 표준편차가 \(3\)으로 작지 않아 런이 \(2\)까지 내려가지는 않지만, 런 검정은 여전히 압도적으로 기각한다.

이것이 런 검정의 잔차 진단으로서의 가치이다. 검정이 잡아낸 것은 오차항의 자기상관이 아니라 모형의 설정 오류이다. 빠진 추세가 잔차에 체계적인 부호 패턴을 남겼고, 런 검정은 그 패턴을 본다.

반대로 옳은 모형을 적합하면 \(R = 23\)으로 기댓값 \(21\)에 가깝고 \(p = 0.52\)이다. 다만 주의할 점이 있다. 여기서 잔차는 자료에서 추정한 두 모수에 의존하므로 엄밀히 독립이 아니다. 순위 기반 귀무분포는 근사일 뿐이며, 잔차의 개수가 적으면 이 근사가 다소 보수적으로 작동한다.