16.2 일표본 비모수 검정¶

import numpy as np
import scipy.stats as stats

def runs_test(data):
    """Wald-Wolfowitz 런 검정. 관측 순서가 무작위인지 본다.

    +1 과 -1 이 늘어선 수열에서 같은 부호가 이어지는 덩어리를 런이라 한다.
    런이 너무 적으면 같은 값끼리 뭉쳐 있다는 뜻이고, 너무 많으면 번갈아
    나온다는 뜻이다. 둘 다 무작위가 아니다.

    매개변수
    --------
    data : +1 과 -1 로 이루어진 수열

    돌려주는 값
    ----------
    statistic : Z 통계량
    p_value : 양측 p-값
    """
    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

import itertools

def count_runs(x):
    x = np.asarray(x)
    return (len(x) + 1 - np.sum(x[1:] * x[:-1])) / 2

print(" N+  N-   배치수    E[R]열거    mu_R공식     Var열거  (mu-1)(mu-2)/(N-1)  교과서꼴")
for N_plus, N_minus in ((6, 11), (9, 8), (4, 4), (5, 3)):
    N = N_plus + N_minus
    mu = 2 * N_plus * N_minus / N + 1
    var_a = (mu - 1) * (mu - 2) / (N - 1)
    var_b = 2 * N_plus * N_minus * (2 * N_plus * N_minus - N) / (N ** 2 * (N - 1))
    Rs = []
    for comb in itertools.combinations(range(N), N_plus):
        x = -np.ones(N, int)
        x[list(comb)] = 1
        Rs.append(count_runs(x))
    Rs = np.array(Rs)
    print(f"{N_plus:>3} {N_minus:>3} {len(Rs):>8} "
          f"{Rs.mean():>11.6f} {mu:>11.6f} {Rs.var():>11.6f} "
          f"{var_a:>19.6f} {var_b:>9.6f}")

 N+  N-   배치수    E[R]열거    mu_R공식     Var열거  (mu-1)(mu-2)/(N-1)  교과서꼴
  6  11    12376    8.764706    8.764706    3.282872            3.282872  3.282872
  9   8    24310    9.470588    9.470588    3.955017            3.955017  3.955017
  4   4       70    5.000000    5.000000    1.714286            1.714286  1.714286
  5   3       56    4.750000    4.750000    1.473214            1.473214  1.473214

data = np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
statistic, p_value = runs_test(data * 2 - 1)  # {0,1} → {-1,+1} 변환
print(f"{statistic = :.4f}")   # -3.7335
print(f"{p_value   = :.4f}")   # 0.0002 → 무작위성 기각

statistic = -3.7335
p_value   = 0.0002

data = np.array([1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0])
statistic, p_value = runs_test(data * 2 - 1)
print(f"{statistic = :.4f}")   # +2.2775
print(f"{p_value   = :.4f}")   # 0.0228 → 무작위성 기각

statistic = 2.2775
p_value   = 0.0228

np.random.seed(1)
data = np.random.binomial(n=1, p=0.4, size=200)
statistic, p_value = runs_test(data * 2 - 1)
print(f"{statistic = :.4f}")   # +1.3294
print(f"{p_value   = :.4f}")   # 0.1837 → 무작위성을 기각하지 못함

statistic = 1.3294
p_value   = 0.1837

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]),
}
np.random.seed(1)
seqs["이항 n=200"] = np.random.binomial(n=1, p=0.4, size=200)

for name, s in seqs.items():
    x = s * 2 - 1
    N = len(x)
    N_plus = int((x == 1).sum())
    N_minus = N - N_plus
    R = (N + 1 - np.sum(x[1:] * x[:-1])) / 2
    mu = 2 * N_plus * N_minus / N + 1
    sigma = np.sqrt((mu - 1) * (mu - 2) / (N - 1))
    z, p = runs_test(x)
    print(f"{name}: N+={N_plus:>3} N-={N_minus:>3} R={R:>5.0f} "
          f"mu={mu:>9.6f} sigma={sigma:>8.6f} "
          f"z={(R - mu) / sigma:>9.6f}  (runs_test: {z:>9.6f}, p={p:.4f})")
    # 가장 긴 런
    best = cur = 1
    for i in range(1, N):
        cur = cur + 1 if s[i] == s[i - 1] else 1
        best = max(best, cur)
    print(f"  가장 긴 런의 길이 {best}")

뭉침: N+=  6 N-= 11 R=    2 mu= 8.764706 sigma=1.811870 z=-3.733550  (runs_test: -3.733550, p=0.0002)
  가장 긴 런의 길이 11
과잉 교대: N+=  9 N-=  8 R=   14 mu= 9.470588 sigma=1.988723 z= 2.277548  (runs_test:  2.277548, p=0.0228)
  가장 긴 런의 길이 2
이항 n=200: N+= 80 N-=120 R=  106 mu=97.000000 sigma=6.769723 z= 1.329449  (runs_test:  1.329449, p=0.1837)
  가장 긴 런의 길이 9

import numpy as np
import scipy.stats as stats

def sign_test(paired_data, test_type="two-sided"):
    """대응표본에 대한 부호검정.

    차이의 크기는 버리고 부호만 센다. 그래서 자료가 순서척도이기만 하면
    쓸 수 있고 이상치에도 끄떡없다. 대신 크기 정보를 버린 만큼 검정력이
    낮다 — 그 중간이 부호순위검정이다.

    매개변수
    --------
    paired_data : 모양 (n, 2) 인 배열. 0열이 처리 후, 1열이 처리 전이다.
    test_type : "less", "two-sided", "greater" 중 하나

    돌려주는 값
    ----------
    z : Z 통계량
    p_value : p-값
    """
    p_0, q_0 = 0.5, 0.5

    # 동점(차이가 0)은 아예 세지 않는다. 그래서 실제로 쓰이는 n 이
    # 원래 표본크기보다 작아진다.
    n_plus = np.sum(paired_data[:, 0] > paired_data[:, 1])
    n_minus = np.sum(paired_data[:, 0] < paired_data[:, 1])
    n = n_plus + n_minus
    p_hat = n_plus / n

    z = (p_hat - p_0) / np.sqrt(p_0 * q_0 / n)

    if test_type == "less":
        p_value = stats.norm.cdf(z)
    elif test_type == "two-sided":
        p_value = 2 * stats.norm.cdf(-abs(z))
    elif test_type == "greater":
        p_value = stats.norm.sf(z)

    return z, p_value

from math import comb

n, n_plus = 12, 10
T = 2 * n_plus - n
print(f"T = n+ - n- = {T},  sqrt(n) = {np.sqrt(n):.6f},  T/sqrt(n) = {T / np.sqrt(n):.6f}")
p_hat = n_plus / n
print(f"(p_hat - 0.5)/sqrt(0.25/n) = {(p_hat - 0.5) / np.sqrt(0.25 / n):.6f}   <- 같다")

exact = 2 * sum(comb(n, k) for k in range(n_plus, n + 1)) / 2 ** n
print(f"\n정확     : 2*({' + '.join(str(comb(n, k)) for k in range(n_plus, n + 1))})"
      f"/{2 ** n} = {exact:.6f}")
print(f"binomtest: {stats.binomtest(n_plus, n).pvalue:.6f}")

z0 = (p_hat - 0.5) / np.sqrt(0.25 / n)
zc = (abs(n_plus - n / 2) - 0.5) / np.sqrt(n / 4)
p0 = 2 * stats.norm.sf(z0)
pc = 2 * stats.norm.sf(zc)
print(f"\n보정 없음: z = {z0:.6f},  p = {p0:.6f},  차이 {p0 - exact:+.6f}")
print(f"연속성보정: z = {zc:.6f},  p = {pc:.6f},  차이 {pc - exact:+.6f}")
print(f"오차가 {abs(p0 - exact) / abs(pc - exact):.1f} 배 줄었다")

T = n+ - n- = 8,  sqrt(n) = 3.464102,  T/sqrt(n) = 2.309401
(p_hat - 0.5)/sqrt(0.25/n) = 2.309401   <- 같다

정확     : 2*(66 + 12 + 1)/4096 = 0.038574
binomtest: 0.038574

보정 없음: z = 2.309401,  p = 0.020921,  차이 -0.017653
연속성보정: z = 2.020726,  p = 0.043308,  차이 +0.004734
오차가 3.7 배 줄었다

# 처리 전후 점수 15쌍. 아래 부호검정과 부호순위검정을 같은 자료에 돌려
# p-값을 견줄 수 있다.
data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]
])
z, p_value = sign_test(data)
print(f"{z       = :.4f}")   # 2.3094
print(f"{p_value = :.4f}")   # 0.0209

z       = 2.3094
p_value = 0.0209

d = data[:, 0] - data[:, 1]
n_plus = int((d > 0).sum())
n_minus = int((d < 0).sum())
n_zero = int((d == 0).sum())
print(f"차이 {list(d)}")
print(f"n+ = {n_plus},  n- = {n_minus},  동점 = {n_zero},  n = {n_plus + n_minus}")

print("\n동점 처리별 p-값")
for label, np_eff, n_eff in (("버린다       ", n_plus, n_plus + n_minus),
                             ("반반씩 센다  ", n_plus + n_zero / 2, 15),
                             ("음의 방향으로", n_plus, 15)):
    p_hat = np_eff / n_eff
    z = (p_hat - 0.5) / np.sqrt(0.25 / n_eff)
    line = (f"  {label}: n+={np_eff:>5} n={n_eff:>3} p_hat={p_hat:.4f} "
            f"z={z:.4f} 근사 p={2 * stats.norm.sf(abs(z)):.4f}")
    if float(np_eff).is_integer():
        line += f" 정확 p={stats.binomtest(int(np_eff), n_eff).pvalue:.4f}"
    print(line)

차이 [17, -2, 6, -3, 14, 0, 6, 10, 12, 8, 9, 2, 0, 3, 0]
n+ = 10,  n- = 2,  동점 = 3,  n = 12

동점 처리별 p-값
  버린다       : n+=   10 n= 12 p_hat=0.8333 z=2.3094 근사 p=0.0209 정확 p=0.0386
  반반씩 센다  : n+= 11.5 n= 15 p_hat=0.7667 z=2.0656 근사 p=0.0389
  음의 방향으로: n+=   10 n= 15 p_hat=0.6667 z=1.2910 근사 p=0.1967 정확 p=0.3018

import numpy as np
import scipy.stats as stats

data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]
])

# 부호순위검정은 부호뿐 아니라 차이의 크기 순위까지 쓴다. 그래서 같은
# 자료에서 부호검정보다 작은 p-값이 나온다.
statistic, p_value = stats.wilcoxon(
    data[:, 0], data[:, 1],
    alternative="two-sided",
    method="approx",      # 정규근사
    zero_method="pratt"   # 0을 순위 매기기에 포함
)
print(f"{statistic = }")   # 11.0  (= W^-)
print(f"{p_value   = :.4f}")   # 0.0086

statistic = 11.0
p_value   = 0.0086

d_orig = data[:, 0] - data[:, 1]
# 부호는 그대로, 크기만 뒤바꾼다. |d| 의 다중집합이 같다.
d_swap = np.array([-17, -14, 12, 10, 9, 8, 6, 6, 3, 3, 2, 2, 0, 0, 0])
print(f"원래  |d| 정렬 {sorted(np.abs(d_orig))}")
print(f"뒤바꾼 |d| 정렬 {sorted(np.abs(d_swap))}   <- 같다")

for label, d in (("원래      ", d_orig), ("크기 뒤바꿈", d_swap)):
    r = stats.rankdata(np.abs(d))
    n_plus, n_minus = int((d > 0).sum()), int((d < 0).sum())
    n = n_plus + n_minus
    z_sign = (n_plus / n - 0.5) / np.sqrt(0.25 / n)
    w = stats.wilcoxon(d, method="approx", zero_method="pratt")
    print(f"\n{label}: n+={n_plus} n-={n_minus}  "
          f"부호검정 z={z_sign:.4f} p={2 * stats.norm.sf(abs(z_sign)):.4f}")
    print(f"            W-={r[d < 0].sum():.1f}  W+={r[d > 0].sum():.1f}  "
          f"부호순위 W={w.statistic:.1f} p={w.pvalue:.6f}")

원래  |d| 정렬 [0, 0, 0, 2, 2, 3, 3, 6, 6, 8, 9, 10, 12, 14, 17]
뒤바꾼 |d| 정렬 [0, 0, 0, 2, 2, 3, 3, 6, 6, 8, 9, 10, 12, 14, 17]   <- 같다

원래      : n+=10 n-=2  부호검정 z=2.3094 p=0.0209
            W-=11.0  W+=103.0  부호순위 W=11.0 p=0.008561

크기 뒤바꿈: n+=10 n-=2  부호검정 z=2.3094 p=0.0209
            W-=29.0  W+=85.0  부호순위 W=29.0 p=0.109526

import numpy as np, scipy.stats as stats
data = np.array([
    [93, 76], [70, 72], [81, 75], [65, 68], [79, 65],
    [54, 54], [94, 88], [91, 81], [77, 65], [65, 57],
    [95, 86], [89, 87], [78, 78], [80, 77], [76, 76]])
d = data[:, 0] - data[:, 1]

print(stats.binomtest(int((d > 0).sum()), int((d != 0).sum())).pvalue)
print(stats.wilcoxon(d, zero_method='pratt', method='approx').pvalue)
print(stats.wilcoxon(d, zero_method='wilcox', method='exact').pvalue)
print(stats.ttest_1samp(d, 0).pvalue)

0.03857421875
0.008560915878575636
0.007579201614253502
0.003815602209766083

import numpy as np
rng = np.random.default_rng(0)

def max_run(s):
    best = cur = 1
    for i in range(1, len(s)):
        cur = cur + 1 if s[i] == s[i - 1] else 1
        best = max(best, cur)
    return best

for n in (17, 50, 100, 200):
    m = [max_run(rng.integers(0, 2, n)) for _ in range(20000)]
    print(n, np.mean(m), np.percentile(m, [5, 95]))

17 4.41955 [3. 7.]
50 5.9724 [4. 9.]
100 6.97655 [ 5. 10.]
200 7.98065 [ 6. 11.]

import numpy as np, scipy.stats as stats
nonzero = [3, 5, 8, 12, 15, 2, 6, -4, -9, -1]
for nz in (0, 10, 20, 30, 40):
    d = np.array([0] * nz + nonzero)
    print(nz,
          round(stats.wilcoxon(d, zero_method='wilcox', method='approx').pvalue, 4),
          round(stats.wilcoxon(d, zero_method='pratt',  method='approx').pvalue, 4),
          round(stats.wilcoxon(d, zero_method='zsplit', method='approx').pvalue, 4))

0 0.1394 0.1394 0.1394
10 0.1394 0.1663 0.1912
20 0.1394 0.1792 0.245
30 0.1394 0.1859 0.2908
40 0.1394 0.1899 0.3297