mu_D에 대한 대응 t-검정¶

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

# 짝마다 차이를 만들고 나면 그 뒤로는 완전히 일표본 문제다.
# 사람마다 점수 수준이 54에서 95까지 흩어져 있어도 그 개인차가 여기서 사라진다.
difference = data[:,0] - data[:,1]
x_bar = difference.mean()
mu = 0
s = difference.std(ddof=1)
n = difference.shape[0]      # 관측값 30개가 아니라 짝 15개

t = (x_bar - mu) / (s / np.sqrt(n))
p_value = 2 * stats.t(n-1).cdf(-abs(t))

print(f"{t       = :.4f}")
print(f"{p_value = :.4f}")

# scipy로 확인. ttest_ind가 아니라 ttest_**rel**이다.
t2, p2 = stats.ttest_rel(data[:,0], data[:,1], alternative="two-sided")
print(f"\nscipy: t = {t2:.4f}, p = {p2:.4f}")

# 짝을 무시하고 독립 이표본으로 다루면 어떻게 되는지 비교한다.
t3, p3 = stats.ttest_ind(data[:,0], data[:,1])
print(f"ttest_ind (틀린 분석): t = {t3:.4f}, p = {p3:.4f}")

t       = 3.4616
p_value = 0.0038

scipy: t = 3.4616, p = 0.0038
ttest_ind (틀린 분석): t = 1.3354, p = 0.1925

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

# 차이가 전부 양수라 음의 순위가 하나도 없다. W = min(W+, W-) = 0이 된다.
stat, p_value = wilcoxon(after, before)

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

alpha = 0.05
if p_value < alpha:
    print("Reject H0: Significant improvement in scores.")
else:
    print("Fail to reject H0: No significant improvement.")

Test Statistic: 0.0
P-value: 0.001953125
Reject H0: Significant improvement in scores.

import numpy as np
from scipy import stats

rng = np.random.default_rng(45)
n, M = 20, 10_000

cases = [("정규", lambda s: rng.normal(0, 1, s), 0.5),
         ("t(3)", lambda s: rng.standard_t(3, s), 0.5),
         ("이중지수", lambda s: rng.laplace(0, 1 / np.sqrt(2), s), 0.5),
         ("로그정규(중심화)",
          lambda s: rng.lognormal(0, 1, s) - np.exp(0.5), 0.8)]
print(f"{'차이의 분포':>16s} {'t':>8s} {'윌콕슨':>9s} {'부호':>8s}")
for name, gen, shift in cases:
    c = np.zeros(3)
    for _ in range(M):
        d = gen(n) + shift
        c[0] += stats.ttest_1samp(d, 0).pvalue < 0.05
        c[1] += stats.wilcoxon(d).pvalue < 0.05
        k = int((d > 0).sum())
        c[2] += stats.binomtest(k, n, 0.5).pvalue < 0.05
    print(f"{name:>16s} {c[0] / M:8.4f} {c[1] / M:9.4f} {c[2] / M:8.4f}")

          차이의 분포        t       윌콕슨       부호
              정규   0.5623    0.5454   0.3826
            t(3)   0.3182    0.3850   0.3297
            이중지수   0.5927    0.6669   0.6290
       로그정규(중심화)   0.3875    0.2819   0.0753

import numpy as np
from scipy import stats

rng = np.random.default_rng(19)
n, M, rho = 20, 5_000, 0.7
L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
hit_cc = hit_ind = 0

for _ in range(M):
    z = rng.standard_normal((n, 2)) @ L.T
    x, y = 10 + z[:, 0], 10.5 + z[:, 1]      # 참 차이 0.5
    miss = rng.random(n) < 0.25              # y 가 25% 결측

    d = (y - x)[~miss]                       # 완전사례 대응 검정
    hit_cc += stats.ttest_1samp(d, 0).pvalue < 0.05

    yy = y[~miss]                            # 있는 자료를 다 쓰되 독립으로
    hit_ind += stats.ttest_ind(yy, x, equal_var=False).pvalue < 0.05

print(f"완전사례 대응 검정   검정력 {hit_cc / M:.4f}")
print(f"짝 깨고 독립 검정    검정력 {hit_ind / M:.4f}")

완전사례 대응 검정   검정력 0.6380
짝 깨고 독립 검정    검정력 0.2096

import numpy as np
from scipy import stats

rng = np.random.default_rng(72)
n, rho = 30, 0.7
L = np.linalg.cholesky(np.array([[1, rho], [rho, 1]]))
z = rng.standard_normal((n, 2)) @ L.T
pre = 50 + 10 * z[:, 0]
post = 53 + 10 * z[:, 1]
d = post - pre

sd_d = d.std(ddof=1)
sd_pre = pre.std(ddof=1)
sd_avg = np.sqrt((pre.var(ddof=1) + post.var(ddof=1)) / 2)
r = np.corrcoef(pre, post)[0, 1]

print(f"평균 차이 {d.mean():.4f}")
print(f"SD(차이) {sd_d:.4f},  SD(사전) {sd_pre:.4f},  "
      f"SD(평균) {sd_avg:.4f},  r = {r:.4f}")
print()
print(f"d_z   = 평균차/SD(차이)  = {d.mean() / sd_d:.4f}")
print(f"d_rm  = d_z × √(2(1-r))  = "
      f"{d.mean() / sd_d * np.sqrt(2 * (1 - r)):.4f}")
print(f"d_av  = 평균차/SD(평균)  = {d.mean() / sd_avg:.4f}")
print(f"d_pre = 평균차/SD(사전)  = {d.mean() / sd_pre:.4f}")

평균 차이 3.1253
SD(차이) 7.0163,  SD(사전) 8.5041,  SD(평균) 9.0700,  r = 0.7060

d_z   = 평균차/SD(차이)  = 0.4454
d_rm  = d_z × √(2(1-r))  = 0.3416
d_av  = 평균차/SD(평균)  = 0.3446
d_pre = 평균차/SD(사전)  = 0.3675

for rho_v in [0.0, 0.5, 0.7, 0.9, 0.95]:
    print(f"ρ = {rho_v:.2f}: d_z / d_av = "
          f"{1 / np.sqrt(2 * (1 - rho_v)):.3f}")

ρ = 0.00: d_z / d_av = 0.707
ρ = 0.50: d_z / d_av = 1.000
ρ = 0.70: d_z / d_av = 1.291
ρ = 0.90: d_z / d_av = 2.236
ρ = 0.95: d_z / d_av = 3.162