대응표본 비모수 검정¶

import numpy as np
from scipy.stats import wilcoxon

# 훈련 전후의 점수. 같은 사람을 두 번 잰 대응자료다.
before = np.array([70, 68, 75, 80, 72, 74, 69, 77, 73, 76])
after = np.array([72, 69, 78, 85, 75, 76, 70, 79, 74, 80])

# 부호순위검정은 대응 t 검정의 비모수 대응이다. 차이의 부호와 크기 순위만
# 쓰므로 정규성을 요구하지 않는다. 열 쌍 모두 점수가 올랐으므로 통계량이
# 0 이 되고, 그때 나올 수 있는 가장 작은 p-값이 (1/2)^10 의 두 배다.
stat, p_value = wilcoxon(after, before)

print(f"Test Statistic: {stat}")     # 0.0
print(f"P-value: {p_value}")         # 0.001953125

# 결과 해석
alpha = 0.05
if p_value < alpha:
    print("Reject the null hypothesis: Significant improvement in scores.")
else:
    print("Fail to reject the null hypothesis: No significant improvement in scores.")

Test Statistic: 0.0
P-value: 0.001953125
Reject the null hypothesis: Significant improvement in scores.

import itertools
from scipy import stats

d = after - before
r = stats.rankdata(np.abs(d))
n = len(d)
print(f"차이        {d.tolist()}")
print(f"|차이| 순위  {r.tolist()}   합 {r.sum()}  (n(n+1)/2 = {n * (n + 1) / 2})")

signs = np.array(list(itertools.product([0, 1], repeat=n)))

def null_T(ranks):
    W_plus = signs @ ranks
    return np.minimum(W_plus, ranks.sum() - W_plus)

T_plain = null_T(np.arange(1, n + 1, dtype=float))   # 동점이 없을 때
T_tied = null_T(r)                                   # 이 자료의 중간순위
print(f"\nT 가 가질 수 있는 값: 동점 없음 {len(set(np.round(T_plain, 6)))} 가지,"
      f"  이 자료 {len(set(np.round(T_tied, 6)))} 가지")
print("  c   동점없음 P(T<=c)   이 자료 P(T<=c)")
for c in (0, 6, 7, 8, 9):
    print(f" {c:>2} {float((T_plain <= c + 1e-9).mean()):>16.6f} "
          f"{float((T_tied <= c + 1e-9).mean()):>17.6f}")
for lab, T in (("동점 없음", T_plain), ("이 자료  ", T_tied)):
    c = max(k for k in range(28) if (T <= k + 1e-9).mean() <= 0.05)
    print(f"  {lab}: 임계값 c = {c},  실제 수준 {float((T <= c + 1e-9).mean()):.6f}")

차이        [2, 1, 3, 5, 3, 2, 1, 2, 1, 4]
|차이| 순위  [5.0, 2.0, 7.5, 10.0, 7.5, 5.0, 2.0, 5.0, 2.0, 9.0]   합 55.0  (n(n+1)/2 = 55.0)

T 가 가질 수 있는 값: 동점 없음 28 가지,  이 자료 43 가지
  c   동점없음 P(T<=c)   이 자료 P(T<=c)
  0         0.001953          0.001953
  6         0.027344          0.021484
  7         0.037109          0.039062
  8         0.048828          0.042969
  9         0.064453          0.062500
  동점 없음: 임계값 c = 8,  실제 수준 0.048828
  이 자료  : 임계값 c = 8,  실제 수준 0.042969

from scipy.stats import binom
import numpy as np

before = np.array([70, 68, 75, 80, 72, 74, 69, 77, 73, 76])
after = np.array([72, 69, 78, 85, 75, 76, 70, 79, 74, 80])

# 같은 자료에 부호검정을 돌린다. 크기는 버리고 부호만 세므로 계산이 더 쉽다.
differences = after - before
n_plus = np.sum(differences > 0)
n_minus = np.sum(differences < 0)
n = n_plus + n_minus
W = min(n_plus, n_minus)

# 귀무가설 아래에서 부호가 +일 확률이 0.5 이므로, 통계량이 이항분포를 따른다.
# 표본이 작을 때는 정규근사 대신 이렇게 정확분포를 쓰는 편이 낫다.
p_value = min(1.0, 2 * binom.cdf(W, n, 0.5))  # 양측

print(f"n+ = {n_plus}, n- = {n_minus}")   # n+ = 10, n- = 0
print(f"P-value: {p_value:.6f}")          # 0.001953

n+ = 10, n- = 0
P-value: 0.001953

d_sorted = np.sort(differences)
print(f"정렬한 차이 {d_sorted.tolist()}")
print(" k   구간          포함률")
for k in range(1, 5):
    cov = 1 - 2 * binom.cdf(k - 1, n, 0.5)
    print(f" {k}  [{d_sorted[k - 1]}, {d_sorted[n - k]}]  "
          f"  {cov:.6f}   (수준 {2 * binom.cdf(k - 1, n, 0.5):.6f})")

k = 2
lo, hi = d_sorted[k - 1], d_sorted[n - k]
level = 2 * binom.cdf(k - 1, n, 0.5)
print(f"\n95% 이상을 보장하는 가장 좁은 구간: k = {k},  [{lo}, {hi}],"
      f"  포함률 {1 - level:.6f}")
print(f"  0 을 포함하는가? {lo <= 0 <= hi}  ->  수준 {level:.6f} 에서 기각")
print(f"  부호검정 p-값 {p_value:.6f} < {level:.6f} ?  {p_value < level}")

# 포함률을 모의로 확인한다 (연속분포, 참 중앙값 0)
rng = np.random.default_rng(1)
B = 2000
for k in (2, 3):
    hit = sum(1 for _ in range(B)
              if (x := np.sort(rng.normal(0, 1, n)))[k - 1] <= 0 <= x[n - k])
    theory = 1 - 2 * binom.cdf(k - 1, n, 0.5)
    print(f"  k={k}: 모의 포함률 {hit / B:.4f},  이론 {theory:.4f},"
          f"  MC 오차 {np.sqrt(theory * (1 - theory) / B):.4f}")

정렬한 차이 [1, 1, 1, 2, 2, 2, 3, 3, 4, 5]
 k   구간          포함률
 1  [1, 5]    0.998047   (수준 0.001953)
 2  [1, 4]    0.978516   (수준 0.021484)
 3  [1, 3]    0.890625   (수준 0.109375)
 4  [2, 3]    0.656250   (수준 0.343750)

95% 이상을 보장하는 가장 좁은 구간: k = 2,  [1, 4],  포함률 0.978516
  0 을 포함하는가? False  ->  수준 0.021484 에서 기각
  부호검정 p-값 0.001953 < 0.021484 ?  True
  k=2: 모의 포함률 0.9810,  이론 0.9785,  MC 오차 0.0032
  k=3: 모의 포함률 0.8905,  이론 0.8906,  MC 오차 0.0070

import numpy as np
from scipy import stats
before = np.array([70, 68, 75, 80, 72, 74, 69, 77, 73, 76])
after  = np.array([72, 69, 78, 85, 75, 76, 70, 79, 74, 80])
d = after - before

print(stats.ttest_rel(after, before).pvalue)                    # 0.000325
print(stats.wilcoxon(d, method='exact').pvalue)                 # 0.001953
print(stats.wilcoxon(d, method='approx', correction=False).pvalue)  # 0.004809
print(stats.binomtest(10, 10).pvalue)                           # 0.001953

0.0003248947130212966
0.001953125
0.004808563078561356
0.001953125

import numpy as np, itertools
n = 10
ranks = np.arange(1, n + 1)
W = np.array([np.dot(s, ranks) for s in itertools.product([0, 1], repeat=n)])
T = np.minimum(W, n * (n + 1) // 2 - W)      # min(W+, W-)
for c in range(6, 12):
    print(c, (T <= c).mean())

6 0.02734375
7 0.037109375
8 0.048828125
9 0.064453125
10 0.083984375
11 0.10546875

from scipy.stats import binom, binomtest
n = 10
for W in range(6):
    naive = 2 * binom.cdf(W, n, 0.5)
    print(W, round(naive, 4), round(binomtest(W, n, 0.5).pvalue, 4))

0 0.002 0.002
1 0.0215 0.0215
2 0.1094 0.1094
3 0.3438 0.3438
4 0.7539 0.7539
5 1.2461 1.0

import numpy as np
from scipy import stats
d = np.array([2, 1, 3, 5, 3, 2, 1, 2, 1, 4])
print(d.mean(), d.std(ddof=1))                          # 2.4, 1.3499
p_neg = stats.norm.cdf(0, d.mean(), d.std(ddof=1))
print(p_neg, (1 - p_neg) ** 10)                         # 0.0377, 0.681

2.4 1.3498971154211057
0.037709049885978514 0.680869711192701