부호검정¶

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

    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, sqrt

crit = stats.norm.isf(0.025)
print(f"임계값 z_0.025 = {crit:.6f}")

# n = 21 을 손으로 따라가 본다.
n = 21
print(f"\nn = {n}: 문턱 = {crit:.6f} * sqrt({n}) = {crit * sqrt(n):.3f}")
ks = [k for k in range(n + 1) if abs(2 * k - n) > crit * sqrt(n)]
print(f"  기각역 k in {ks}")
terms = [comb(n, k) for k in range(7)]
print(f"  2*({' + '.join(map(str, terms))})/{2 ** n} = {2 * sum(terms)}/{2 ** n}"
      f" = {2 * sum(terms) / 2 ** n:.6f}")

print("\n  n   보정없음    보정    정확")
worst = {"plain": (0, 0), "corr": (0, 0), "exact": (0, 0)}
over = []
for n in range(5, 41):
    pmf = np.array([comb(n, k) for k in range(n + 1)]) / 2 ** n
    k = np.arange(n + 1)
    rej_plain = np.abs(2 * k - n) / np.sqrt(n) > crit
    rej_corr = (np.abs(k - n / 2) - 0.5) / np.sqrt(n / 4) > crit
    rej_exact = np.array([stats.binomtest(j, n).pvalue < 0.05 for j in k])
    s = (pmf[rej_plain].sum(), pmf[rej_corr].sum(), pmf[rej_exact].sum())
    for key, val in zip(("plain", "corr", "exact"), s):
        if val > worst[key][0]:
            worst[key] = (val, n)
    if s[0] > 0.05:
        over.append(n)
    if n in (5, 12, 17, 21, 25, 30, 39, 40):
        print(f"{n:>3} {s[0]:>10.4f} {s[1]:>8.4f} {s[2]:>8.4f}")

print(f"\n보정 없음: 명목 0.05 를 넘는 n {len(over)} 개 / 36 개  -> {over}")
print(f"  최대 실제 수준 {worst['plain'][0]:.6f} (n = {worst['plain'][1]})")
print(f"연속성 보정: 최대 {worst['corr'][0]:.6f} (n = {worst['corr'][1]})")
print(f"정확 이항  : 최대 {worst['exact'][0]:.6f} (n = {worst['exact'][1]})")

임계값 z_0.025 = 1.959964

n = 21: 문턱 = 1.959964 * sqrt(21) = 8.982
  기각역 k in [0, 1, 2, 3, 4, 5, 6, 15, 16, 17, 18, 19, 20, 21]
  2*(1 + 21 + 210 + 1330 + 5985 + 20349 + 54264)/2097152 = 164320/2097152 = 0.078354

  n   보정없음    보정    정확
  5     0.0625   0.0000   0.0000
 12     0.0386   0.0386   0.0386
 17     0.0490   0.0127   0.0490
 21     0.0784   0.0266   0.0266
 25     0.0433   0.0433   0.0433
 30     0.0428   0.0428   0.0428
 39     0.0533   0.0237   0.0237
 40     0.0385   0.0385   0.0385

보정 없음: 명목 0.05 를 넘는 n 17 개 / 36 개  -> [5, 8, 11, 14, 16, 19, 21, 22, 24, 26, 27, 29, 31, 32, 34, 36, 39]
  최대 실제 수준 0.078354 (n = 21)
연속성 보정: 최대 0.047031 (n = 37)
정확 이항  : 최대 0.049042 (n = 17)

paired_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]
])

diffs = paired_data[:, 0] - paired_data[:, 1]
print(diffs)
# [17 -2  6 -3 14  0  6 10 12  8  9  2  0  3  0]

z, p_value = sign_test(paired_data, test_type="two-sided")
print(f"Z = {z:.4f}, p = {p_value:.4f}")
# Z = 2.3094, p = 0.0209

[17 -2  6 -3 14  0  6 10 12  8  9  2  0  3  0]
Z = 2.3094, p = 0.0209

from scipy import stats
print(stats.binomtest(10, 12, alternative='greater').pvalue)   # 0.019287

0.019287109375

import numpy as np
from scipy.stats import binom

def sign_test_exact(paired_data):
    """대응관측에 대한 정확 양측 부호검정.

        정규근사 대신 이항분포로 p-값을 정확히 구한다. 표본이 작을 때는
        이쪽이 맞다.
        """
    diffs = paired_data[:, 0] - paired_data[:, 1]
    nonzero = diffs[diffs != 0]
    n = len(nonzero)
    n_plus = (nonzero > 0).sum()

    # 양측 p-값: 2 * min(P(X <= n_+), P(X >= n_+))
    p_left = binom.cdf(n_plus, n, 0.5)
    p_right = binom.sf(n_plus - 1, n, 0.5)  # P(X >= n_+)
    p_value = 2 * min(p_left, p_right)
    p_value = min(p_value, 1.0)  # cap at 1

    return n_plus, n, p_value

paired_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]
])

n_plus, n, p = sign_test_exact(paired_data)
print(f"n_+ = {n_plus}, n = {n}, exact p = {p:.4f}")
# n_+ = 10, n = 12, exact p = 0.0386

n_+ = 10, n = 12, exact p = 0.0386

from scipy.stats import binomtest
print(binomtest(10, 12).pvalue)   # 0.03857421875

0.03857421875

import numpy as np
from scipy import stats
d = np.array([17, -2, 6, -3, 14, 0, 6, 10, 12, 8, 9, 2, 0, 3, 0])

print(stats.binomtest(10, 12).pvalue)                        # 0.038574
print(stats.wilcoxon(d, zero_method='wilcox',
                     method='exact').pvalue)                 # 0.007579
print(stats.ttest_1samp(d, 0).pvalue)                        # 0.003816

0.03857421875
0.007579201614253502
0.003815602209766083