런 검정¶

import numpy as np
import scipy.stats as stats


def runs_test(data):
    """+1/-1 수열에 대한 Wald-Wolfowitz 런 검정.

    같은 부호가 이어지는 덩어리를 런이라 한다. 런이 너무 적으면 뭉쳐
    있다는 뜻이고 너무 많으면 번갈아 난다는 뜻이며, 둘 다 무작위가 아니다.
    """
    data = np.asarray(data)
    N = data.shape[0]
    N_plus = (data == 1).sum()
    N_minus = N - N_plus

    # 무작위라는 가정 아래에서 런 개수의 평균과 표준편차
    mu = 2 * N_plus * N_minus / N + 1
    sigma = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    # 이웃한 원소의 곱은 부호가 바뀌는 자리에서만 -1 이다. 그 개수를 세면
    # 런의 경계 수가 나온다.
    R = (N_plus + N_minus + 1 - np.sum(data[1:] * data[:-1])) / 2

    statistic = (R - mu) / sigma
    p_value = 2 * stats.norm.cdf(-abs(statistic))
    return statistic, p_value

def runs_detail(data):
    """R 과 중간값을 모두 돌려주어 추적할 수 있게 한다."""
    data = np.asarray(data)
    N = data.shape[0]
    N_plus = int((data == 1).sum())
    N_minus = N - N_plus
    mu = 2 * N_plus * N_minus / N + 1
    sigma = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    prod = int(np.sum(data[1:] * data[:-1]))
    R = (N + 1 - prod) / 2
    z = (R - mu) / sigma
    return N_plus, N_minus, prod, R, mu, sigma, z, 2 * stats.norm.cdf(-abs(z))

seqs = {
    "뭉침": np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
    "과잉 교대": np.array([1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0]),
}
for name, s in seqs.items():
    for tag, x in (("±1 로 변환", s * 2 - 1), ("{0,1} 그대로", s)):
        Np, Nm, prod, R, mu, sigma, z, p = runs_detail(x)
        print(f"{name} / {tag}: N+={Np} N-={Nm} sum(x_i x_(i+1))={prod:>3} "
              f"R={R:>5} mu={mu:.4f} z={z:>9.4f} p={p:.4f}")

뭉침 / ±1 로 변환: N+=6 N-=11 sum(x_i x_(i+1))= 14 R=  2.0 mu=8.7647 z=  -3.7335 p=0.0002
뭉침 / {0,1} 그대로: N+=6 N-=11 sum(x_i x_(i+1))=  5 R=  6.5 mu=8.7647 z=  -1.2499 p=0.2113
과잉 교대 / ±1 로 변환: N+=9 N-=8 sum(x_i x_(i+1))=-10 R= 14.0 mu=9.4706 z=   2.2775 p=0.0228
과잉 교대 / {0,1} 그대로: N+=9 N-=8 sum(x_i x_(i+1))=  2 R=  8.0 mu=9.4706 z=  -0.7395 p=0.4596

import itertools
from collections import Counter
from math import comb

def exact_runs_pmf(N_plus, N_minus):
    """가능한 모든 배치를 열거해 R 의 정확 귀무분포를 얻는다."""
    N = N_plus + N_minus
    tally = Counter()
    for cb in itertools.combinations(range(N), N_plus):
        x = -np.ones(N, int)
        x[list(cb)] = 1
        tally[int((N + 1 - np.sum(x[1:] * x[:-1])) // 2)] += 1
    return tally, comb(N, N_plus)

for name, N_plus, N_minus, R_obs in (("뭉침", 6, 11, 2), ("과잉 교대", 9, 8, 14)):
    N = N_plus + N_minus
    mu = 2 * N_plus * N_minus / N + 1
    sigma = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    tally, total = exact_runs_pmf(N_plus, N_minus)
    dev = abs(R_obs - mu)
    hit = sum(c for r, c in tally.items() if abs(r - mu) >= dev - 1e-12)
    lower = sum(c for r, c in tally.items() if r <= R_obs)
    upper = sum(c for r, c in tally.items() if r >= R_obs)
    p_two = min(1.0, 2 * min(lower, upper) / total)
    print(f"{name}: N+={N_plus} N-={N_minus}  배치 {total} 가지,  "
          f"R 의 범위 {min(tally)} ~ {max(tally)}")
    print(f"   정규근사      p = {2 * stats.norm.cdf(-abs((R_obs - mu) / sigma)):.6f}")
    print(f"   편차 기준     p = {hit}/{total} = {hit / total:.6f}")
    print(f"   작은 꼬리 x2  p = {p_two:.6f}")

뭉침: N+=6 N-=11  배치 12376 가지,  R 의 범위 2 ~ 13
   정규근사      p = 0.000189
   편차 기준     p = 2/12376 = 0.000162
   작은 꼬리 x2  p = 0.000323
과잉 교대: N+=9 N-=8  배치 24310 가지,  R 의 범위 2 ~ 17
   정규근사      p = 0.022753
   편차 기준     p = 622/24310 = 0.025586
   작은 꼬리 x2  p = 0.040559

data = np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
z, p = runs_test(data * 2 - 1)
print(f"z = {z:.4f}, p = {p:.4f}")
# R = 2, mu = 8.76, sigma = 1.81
# z = -3.7335, p = 0.0002 → 무작위성 기각

z = -3.7335, p = 0.0002

tally, total = exact_runs_pmf(6, 11)
print(f"배치 {total} 가지,  R 의 범위 {min(tally)} ~ {max(tally)}"
      f"  (2*min(N+,N-)+1 = {2 * min(6, 11) + 1})")
print("R :", [r for r in sorted(tally)])
print("개수:", [tally[r] for r in sorted(tally)])
print(f"R = 2 인 배치 {tally[2]} 가지  ->  (+^6 -^11) 과 (-^11 +^6)")
print(f"R = 13 인 배치 {tally[13]} 가지  (C(10, 6) = {comb(10, 6)})")
print(f"합계 {sum(tally.values())} = {total}")
print(f"\n편차 기준 정확 p = {tally[2]}/{total} = {tally[2] / total:.6f}")
print(f"정규근사       p = 0.000189,  비 = {0.000189 / (tally[2] / total):.2f} 배")

배치 12376 가지,  R 의 범위 2 ~ 13  (2*min(N+,N-)+1 = 13)
R : [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
개수: [2, 15, 100, 325, 900, 1650, 2400, 2700, 2100, 1470, 504, 210]
R = 2 인 배치 2 가지  ->  (+^6 -^11) 과 (-^11 +^6)
R = 13 인 배치 210 가지  (C(10, 6) = 210)
합계 12376 = 12376

편차 기준 정확 p = 2/12376 = 0.000162
정규근사       p = 0.000189,  비 = 1.17 배

data = np.array([1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0])
z, p = runs_test(data * 2 - 1)
# R = 14, mu = 9.47, sigma = 1.99
# z = +2.2775, p = 0.0228 → 무작위성 기각

import numpy as np
import scipy.stats as stats

data = np.array([1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0])
seq = data * 2 - 1  # convert to +1/-1

N = len(seq)
N_plus = (seq == 1).sum()   # 9
N_minus = N - N_plus        # 8

mu = 2 * N_plus * N_minus / N + 1               # 9.4706
sigma = np.sqrt((mu - 1) * (mu - 2) / (N - 1))  # 1.9887
R = (N - np.sum(seq[1:] * seq[:-1]) + 1) / 2    # 14

z = (R - mu) / sigma
p = 2 * stats.norm.cdf(-abs(z))
print(f"R = {R}, Z = {z:.4f}, p = {p:.4f}")
# R = 14.0, Z = 2.2775, p = 0.0228  →  H0 기각

R = 14.0, Z = 2.2775, p = 0.0228

import numpy as np
from math import comb
import scipy.stats as stats

def runs_pmf(N_plus, N_minus):
    N = N_plus + N_minus
    tot = comb(N, N_plus)
    pmf = {}
    for k in range(1, min(N_plus, N_minus) + 1):
        pmf[2 * k] = 2 * comb(N_plus - 1, k - 1) * comb(N_minus - 1, k - 1) / tot
        pmf[2 * k + 1] = (comb(N_plus - 1, k - 1) * comb(N_minus - 1, k)
                          + comb(N_plus - 1, k) * comb(N_minus - 1, k - 1)) / tot
    return pmf

crit = stats.norm.isf(0.025)
for n in (5, 10, 15, 20, 50):
    N = 2 * n
    mu = 2 * n * n / N + 1
    sd = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    pmf = runs_pmf(n, n)
    size = sum(p for R, p in pmf.items() if abs((R - mu) / sd) > crit)
    print(N, round(size, 4))

10 0.0794
20 0.037
30 0.0398
40 0.0363
100 0.0555