이표본 검정¶

import numpy as np
from scipy import stats

n_f, n_s = 100, 100
x_bar_f, x_bar_s = 1.85, 1.65
s_f, s_s = 1.3, 1.2

# n이 각각 100이라 크므로 sigma 자리에 표본표준편차를 넣고 z를 쓴다.
statistic = (x_bar_f - x_bar_s) / np.sqrt(s_f**2/n_f + s_s**2/n_s)
p_value = stats.norm().sf(abs(statistic)) * 2

print(f"statistic : {statistic:.4f}")
print(f"p value   : {p_value:.4f}")

statistic : 1.1305
p value   : 0.2583

import numpy as np
from scipy import stats

se = np.sqrt(1.3**2 / 100 + 1.2**2 / 100)
z_crit = stats.norm.ppf(0.975)
print(f"SE             = {se:.6f}")
print(f"기각 경계 차이 = {z_crit * se:.6f}")

# 0.2 의 차를 양측 5% 검정으로 80% 확률로 잡아내려면 집단당 몇 명인가.
n_formula = (z_crit + stats.norm.ppf(0.8)) ** 2 * (1.3**2 + 1.2**2) / 0.2**2
print(f"공식이 주는 n  = {n_formula:.2f}")
for n in (613, 614, 615):
    s = np.sqrt(3.13 / n)
    power = stats.norm.sf(z_crit - 0.2 / s) + stats.norm.cdf(-z_crit - 0.2 / s)
    print(f"  n = {n}: 검정력 {power:.4f}")

SE             = 0.176918
기각 경계 차이 = 0.346753
공식이 주는 n  = 614.17
  n = 613: 검정력 0.7992
  n = 614: 검정력 0.7999
  n = 615: 검정력 0.8005

import numpy as np
from scipy import stats

n_1, n_2 = 30, 25
X_1_bar, X_2_bar = 62.5, 58.9
s_1, s_2 = 7.2, 6.4
df = n_1 + n_2 - 2

s_p_square = ((n_1 - 1) * s_1**2 + (n_2 - 1) * s_2**2) / df
statistic = (X_1_bar - X_2_bar) / np.sqrt(s_p_square * (1 / n_1 + 1 / n_2))

print(f"s_p       = {np.sqrt(s_p_square):.4f}")
print(f"statistic = {statistic:.4f}   (df = {df})")
print(f"단측 p    = {stats.t(df).sf(statistic):.4f}   임계값 {stats.t(df).ppf(0.95):.4f}")
print(f"양측 p    = {2 * stats.t(df).sf(statistic):.4f}   임계값 {stats.t(df).ppf(0.975):.4f}")

# 방향을 자료를 보고 고르면 실제 수준은 2alpha 가 된다.
rng = np.random.default_rng(7)
B = 200_000
x = rng.normal(0, 1, (B, n_1))
y = rng.normal(0, 1, (B, n_2))
v = ((n_1 - 1) * x.var(1, ddof=1) + (n_2 - 1) * y.var(1, ddof=1)) / df
t = (x.mean(1) - y.mean(1)) / np.sqrt(v * (1 / n_1 + 1 / n_2))
c95 = stats.t(df).ppf(0.95)
print(f"\nH0 가 참일 때 {B}번 중 기각 비율")
print(f"  방향을 미리 고정한 단측 : {np.mean(t > c95):.4f}")
print(f"  자료를 보고 방향을 고름 : {np.mean(np.abs(t) > c95):.4f}")
print(f"  양측                    : {np.mean(np.abs(t) > stats.t(df).ppf(0.975)):.4f}")

s_p       = 6.8493
statistic = 1.9409   (df = 53)
단측 p    = 0.0288   임계값 1.6741
양측 p    = 0.0576   임계값 2.0057

H0 가 참일 때 200000번 중 기각 비율
  방향을 미리 고정한 단측 : 0.0500
  자료를 보고 방향을 고름 : 0.1008
  양측                    : 0.0503

import numpy as np
from scipy import stats

X_1_bar, X_2_bar = 1.3, 1.6
s_1, s_2 = 0.5, 0.3
n_1, n_2 = 22, 24

# 표준편차가 0.5와 0.3으로 다르므로 합동하지 않고 Welch를 쓴다.
statistic = (X_1_bar - X_2_bar) / np.sqrt(s_1**2 / n_1 + s_2**2 / n_2)

# Welch-Satterthwaite 자유도.
# 분모에 n_i가 아니라 **n_i - 1**이 들어간다는 점에 주의하라.
# n으로 잘못 쓰면 자유도가 부풀어 기각하기 쉬워진다.
top = (s_1**2 / n_1 + s_2**2 / n_2)**2
bottom = (s_1**2 / n_1)**2 / (n_1 - 1) + (s_2**2 / n_2)**2 / (n_2 - 1)
df = top / bottom

p_value = 2 * stats.t(df).cdf(-abs(statistic))
print(f"{df = :.4f}")
print(f"{statistic = :.4f}")
print(f"{p_value   = :.4f}")

alpha = 0.05
if p_value <= alpha:
    print("Reject H_0")
else:
    print("Fail to reject H_0")

df = 33.7874
statistic = -2.4403
p_value   = 0.0201
Reject H_0

import numpy as np
from scipy import stats

s_1, s_2 = 0.5, 0.3
n_1, n_2 = 22, 24

# 같은 자료에 합동 t 를 걸어 본다.
s_p_square = ((n_1 - 1) * s_1**2 + (n_2 - 1) * s_2**2) / (n_1 + n_2 - 2)
se_pool = np.sqrt(s_p_square * (1 / n_1 + 1 / n_2))
se_welch = np.sqrt(s_1**2 / n_1 + s_2**2 / n_2)
t_pool = (1.3 - 1.6) / se_pool

print(f"웰치 SE   = {se_welch:.6f}")
print(f"합동 SE   = {se_pool:.6f}   비 = {se_pool / se_welch:.6f}")
R = s_p_square * (1 / n_1 + 1 / n_2) / (s_1**2 / n_1 + s_2**2 / n_2)
print(f"R = {R:.6f}   sqrt(R) = {np.sqrt(R):.6f}")
print(f"합동 t    = {t_pool:.4f}  (df = {n_1 + n_2 - 2})"
      f"  p = {2 * stats.t(n_1 + n_2 - 2).cdf(-abs(t_pool)):.4f}")

웰치 SE   = 0.122938
합동 SE   = 0.120390   비 = 0.979280
R = 0.958988   sqrt(R) = 0.979280
합동 t    = -2.4919  (df = 44)  p = 0.0165

X_1_bar, X_2_bar = 1.85, 1.65
s_1, s_2 = 1.3, 1.2
n_1, n_2 = 100, 100

# 합동분산은 두 표본분산을 자유도로 가중평균한 것이다.
# 여기서는 n이 같아 단순 평균과 같아진다.
s_p_square = ((n_1 - 1) * s_1**2 + (n_2 - 1) * s_2**2) / (n_1 + n_2 - 2)
statistic = (X_1_bar - X_2_bar) / np.sqrt(s_p_square / n_1 + s_p_square / n_2)
df = n_1 + n_2 - 2
p_value = 2 * stats.t(df).cdf(-abs(statistic))

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

df = 198.0000
statistic = 1.1305
p_value   = 0.2596

import numpy as np
from scipy import stats

s_1, s_2, n = 1.3, 1.2, 100
s_p_square = ((n - 1) * s_1**2 + (n - 1) * s_2**2) / (2 * n - 2)

print(f"s_p^2            = {s_p_square:.6f}")
print(f"(s_1^2+s_2^2)/2  = {(s_1**2 + s_2**2) / 2:.6f}")
print(f"합동 SE          = {np.sqrt(s_p_square * 2 / n):.10f}")
print(f"z 검정의 분모    = {np.sqrt(s_1**2 / n + s_2**2 / n):.10f}")
print(f"t_.975,198 = {stats.t(2 * n - 2).ppf(0.975):.6f}"
      f"    z_.975 = {stats.norm.ppf(0.975):.6f}")

s_p^2            = 1.565000
(s_1^2+s_2^2)/2  = 1.565000
합동 SE          = 0.1769180601
z 검정의 분모    = 0.1769180601
t_.975,198 = 1.972017    z_.975 = 1.959964

import numpy as np
from scipy import stats

X_1_bar, X_2_bar = 120, 116
s_1, s_2 = 15, 13
n_1 = n_2 = 65
diff = X_1_bar - X_2_bar

# n_1 = n_2 이므로 합동 SE 와 웰치 SE 가 정확히 같다.
se = np.sqrt(s_1**2 / n_1 + s_2**2 / n_2)
df_pool = n_1 + n_2 - 2
a, b = s_1**2 / n_1, s_2**2 / n_2
df_welch = (a + b) ** 2 / (a**2 / (n_1 - 1) + b**2 / (n_2 - 1))

print(f"SE = {se:.6f}")
for name, df in (("합동", df_pool), ("웰치", df_welch)):
    m = stats.t(df).ppf(0.995) * se
    print(f"{name}  df = {df:7.2f}  t_.995 = {stats.t(df).ppf(0.995):.4f}"
          f"  여유폭 = {m:.4f}  99% CI = ({diff - m:.4f}, {diff + m:.4f})")

t = diff / se
print(f"\nt = {t:.4f}   양측 p = {2 * stats.t(df_pool).sf(t):.4f}")
m95 = stats.t(df_pool).ppf(0.975) * se
print(f"95% CI = ({diff - m95:.4f}, {diff + m95:.4f})")

SE = 2.462019
합동  df =  128.00  t_.995 = 2.6148  여유폭 = 6.4377  99% CI = (-2.4377, 10.4377)
웰치  df =  125.47  t_.995 = 2.6156  여유폭 = 6.4396  99% CI = (-2.4396, 10.4396)

t = 1.6247   양측 p = 0.1067
95% CI = (-0.8715, 8.8715)

import numpy as np
from scipy.stats import ttest_ind

team_a = [120, 118, 125, 130, 115, 122, 121, 119, 117, 123, 124, 126, 127, 118, 116]
team_b = [135, 132, 137, 140, 136, 130, 134, 138, 139, 133, 131, 142, 141,
           129, 128, 135, 137, 136, 134, 132]

# 표본크기가 15와 20으로 다르다. 이런 상황이 Welch를 쓸 이유다.
stat, p_value = ttest_ind(team_a, team_b, equal_var=False)
print(f"Test Statistic: {stat:.4f}")
print(f"P-value: {p_value:.4f}")

alpha = 0.05
if p_value < alpha:
    print("Reject H0: The means are significantly different.")
else:
    print("Fail to reject H0.")

Test Statistic: -9.4407
P-value: 0.0000
Reject H0: The means are significantly different.

from scipy import stats

a = np.array(team_a, dtype=float)
b = np.array(team_b, dtype=float)
n1, n2 = len(a), len(b)
N = n1 + n2
v1, v2 = a.var(ddof=1), b.var(ddof=1)
print(f"n = {n1}, {n2}   평균 {a.mean():.4f}, {b.mean():.4f}   "
      f"표준편차 {np.sqrt(v1):.4f}, {np.sqrt(v2):.4f}")
print(f"분산비 s1^2/s2^2 = {v1 / v2:.6f}")

se_w = np.sqrt(v1 / n1 + v2 / n2)
sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (N - 2)
se_p = np.sqrt(sp2 * (1 / n1 + 1 / n2))
ident = (N - 1) * (n1 - n2) * (v1 - v2) / ((N - 2) * n1 * n2)
print(f"\nSE_welch = {se_w:.9f}   SE_pool = {se_p:.9f}")
print(f"SE_p^2 - SE_w^2 = {se_p**2 - se_w**2:.9f}   닫힌 꼴 = {ident:.9f}")
print(f"(n1-n2)(s1^2-s2^2) = ({n1 - n2}) x ({v1 - v2:+.4f}) = "
      f"{(n1 - n2) * (v1 - v2):+.4f}  →  합동 SE 가 더 작다")

d = a.mean() - b.mean()
nu = (v1 / n1 + v2 / n2) ** 2 / ((v1 / n1) ** 2 / (n1 - 1) + (v2 / n2) ** 2 / (n2 - 1))
print(f"\n{'검정':>8}{'t':>12}{'자유도':>10}{'p':>14}")
print(f"{'Welch':>8}{d / se_w:12.6f}{nu:10.4f}"
      f"{2 * stats.t.sf(abs(d / se_w), nu):14.3e}")
print(f"{'합동':>8}{d / se_p:12.6f}{N - 2:10d}"
      f"{2 * stats.t.sf(abs(d / se_p), N - 2):14.3e}")
print(f"자유도 범위 검사: {min(n1, n2) - 1} <= {nu:.4f} <= {N - 2}  "
      f"→ {min(n1, n2) - 1 <= nu <= N - 2}")
print(f"코헨의 d = {d / np.sqrt(sp2):.6f}")

n = 15, 20   평균 121.4000, 134.9500   표준편차 4.3883, 3.9400
분산비 s1^2/s2^2 = 1.240501

SE_welch = 1.435267827   SE_pool = 1.412757530
SE_p^2 - SE_w^2 = -0.064109896   닫힌 꼴 = -0.064109896
(n1-n2)(s1^2-s2^2) = (-5) x (+3.7335) = -18.6673  →  합동 SE 가 더 작다

      검정           t       자유도             p
   Welch   -9.440747   28.3975     2.942e-10
      합동   -9.591172        33     4.556e-11
자유도 범위 검사: 14 <= 28.3975 <= 33  → True
코헨의 d = -3.276009

import numpy as np
from scipy import stats

positive_A, positive_B = 58, 52
n_A, n_B = 100, 100
p_hat_A, p_hat_B = positive_A / n_A, positive_B / n_B

# H0가 "두 비율이 같다"이므로 그 공통값을 전체를 합쳐 추정한다.
# 신뢰구간을 만들 때는 이렇게 합동하지 않는다. 목적이 다르기 때문이다.
p_pooled = (positive_A + positive_B) / (n_A + n_B)
statistic = (p_hat_A - p_hat_B) / (np.sqrt(p_pooled * (1 - p_pooled)) * np.sqrt(1/n_A + 1/n_B))
p_value = stats.norm().sf(abs(statistic)) * 2

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

alpha = 0.05
if p_value <= alpha:
    print("Reject H_0")
else:
    print("Fail to reject H_0")

statistic = 0.8528
p_value = 0.3938
Fail to reject H_0

p1, p2 = 0.58, 0.52
pbar = (p1 + p2) / 2
za, zb = stats.norm.ppf(0.975), stats.norm.ppf(0.80)
n_formula = ((za * np.sqrt(2 * pbar * (1 - pbar))
              + zb * np.sqrt(p1 * (1 - p1) + p2 * (1 - p2))) / (p1 - p2)) ** 2
print(f"공식이 주는 지구당 n = {n_formula:.4f}  →  {int(np.ceil(n_formula))}")


def power_two_prop(m, p1, p2, alpha=0.05, one_sided=False):
    """각 집단 m 명인 합동 z 검정의 검정력 (정규근사)."""
    pbar = (p1 + p2) / 2
    se0 = np.sqrt(2 * pbar * (1 - pbar) / m)      # H0 아래의 표준오차
    se1 = np.sqrt(p1 * (1 - p1) / m + p2 * (1 - p2) / m)
    zc = stats.norm.ppf(1 - alpha) if one_sided else stats.norm.ppf(1 - alpha / 2)
    return stats.norm.sf((zc * se0 - abs(p1 - p2)) / se1)


print("\n   n(지구당)   검정력")
for m in (100, 500, 1000, 1078, 1079, 1200):
    print(f"{m:9d}   {power_two_prop(m, p1, p2):.5f}")
ms = np.arange(10, 5001)
pw = np.array([power_two_prop(int(m), p1, p2) for m in ms])
ok = ms[pw >= 0.80]
print(f"\n검정력 >= 0.80 인 가장 작은 n = {ok.min()},  연속인가 = "
      f"{bool(np.all(np.diff(ok) == 1))}")
print(f"지구당 100명 설계의 검정력 = {power_two_prop(100, p1, p2):.4f}")

공식이 주는 지구당 n = 1078.0413  →  1079

   n(지구당)   검정력
      100   0.13368
      500   0.47881
     1000   0.76980
     1078   0.79998
     1079   0.80035
     1200   0.84039

검정력 >= 0.80 인 가장 작은 n = 1079,  연속인가 = True
지구당 100명 설계의 검정력 = 0.1337

# 11월과 12월에 각 m 명을 새로 조사하는 설계. 지지율이 0.50 에서 0.45 로
# 떨어진 것을 단측 5% 에서 80% 검정력으로 잡으려면 m 이 얼마여야 하는가.
p_nov, p_dec = 0.50, 0.45
pbar = (p_nov + p_dec) / 2
z_a, z_b = stats.norm.ppf(0.95), stats.norm.ppf(0.80)
m_formula = ((z_a * np.sqrt(2 * pbar * (1 - pbar))
              + z_b * np.sqrt(p_nov * (1 - p_nov) + p_dec * (1 - p_dec)))
             / (p_nov - p_dec)) ** 2
print(f"단측 공식이 주는 매달 n = {m_formula:.4f}  →  {int(np.ceil(m_formula))}")
ms = np.arange(10, 5001)
pw = np.array([power_two_prop(int(m), p_nov, p_dec, one_sided=True) for m in ms])
ok = ms[pw >= 0.80]
print(f"전수 탐색: 가장 작은 n = {ok.min()},  연속인가 = "
      f"{bool(np.all(np.diff(ok) == 1))}")
print(f"\n{'n(매달)':>8}{'단측 검정력':>12}{'양측 검정력':>12}")
for m in (500, 1000, 1233, 1500):
    print(f"{m:>8}{power_two_prop(m, p_nov, p_dec, one_sided=True):>12.4f}"
          f"{power_two_prop(m, p_nov, p_dec):>12.4f}")

# 같은 사람을 두 번 조사했다면: 대응 설계의 표준오차
print("\n같은 1000명을 두 번 조사했다면 (대응 설계)")
print(f"{'rho':>6}{'독립 SE':>10}{'대응 SE':>10}{'비':>8}")
se_ind = np.sqrt(p_nov * (1 - p_nov) / 1000 + p_dec * (1 - p_dec) / 1000)
for rho in (0.0, 0.3, 0.6, 0.9):
    se_pair = np.sqrt(p_nov * (1 - p_nov) / 1000 + p_dec * (1 - p_dec) / 1000
                      - 2 * rho * np.sqrt(p_nov * (1 - p_nov) * p_dec * (1 - p_dec)) / 1000)
    print(f"{rho:>6.1f}{se_ind:>10.6f}{se_pair:>10.6f}{se_pair / se_ind:>8.4f}")

단측 공식이 주는 매달 n = 1232.3734  →  1233
전수 탐색: 가장 작은 n = 1233,  연속인가 = True

   n(매달)      단측 검정력      양측 검정력
     500      0.4754      0.3530
    1000      0.7240      0.6100
    1233      0.8002      0.7008
    1500      0.8640      0.7832

같은 1000명을 두 번 조사했다면 (대응 설계)
   rho     독립 SE     대응 SE       비
   0.0  0.022305  0.022305  1.0000
   0.3  0.022305  0.018662  0.8367
   0.6  0.022305  0.014107  0.6325
   0.9  0.022305  0.007054  0.3162

# 동전 두 개를 각 n 번 던져 비율 차이를 재는 설계.
# H0 아래에서 p1 = p2 = 1/2 이므로 차이의 표준오차가 닫힌 꼴로 나온다.
print(f"{'n(동전마다)':>12}{'SE(차이)':>12}{'95% 구간의 반너비':>18}")
for n in (25, 100, 200, 768, 769, 1000):
    se = np.sqrt(0.25 / n + 0.25 / n)
    print(f"{n:>12}{se:>12.6f}{1.959964 * se:>18.6f}")
need = 1.959964 ** 2 * 0.5 / 0.05 ** 2
print(f"\n반너비 <= 0.05 에 필요한 n = 1.96^2 * 0.5 / 0.05^2 = {need:.4f}"
      f"  →  {int(np.ceil(need))}")

# "두 비율이 같은가" 와 "둘 다 1/2 인가" 는 다른 물음이다.
print("\n두 동전이 모두 0.6 으로 치우쳐 있다면 (n = 400)")
k1, k2 = 240, 240                      # 둘 다 60%
p1, p2 = k1 / 400, k2 / 400
P = (k1 + k2) / 800
z_two = (p1 - p2) / np.sqrt(P * (1 - P) * (1 / 400 + 1 / 400))
print(f"  이표본 검정 (H0: p1 = p2):  z = {z_two:.4f}, "
      f"p = {2 * stats.norm.sf(abs(z_two)):.4f}  → 기각 못 함")
z_one = (p1 - 0.5) / np.sqrt(0.25 / 400)
print(f"  일표본 검정 (H0: p1 = 1/2): z = {z_one:.4f}, "
      f"p = {2 * stats.norm.sf(abs(z_one)):.3e}  → 강하게 기각")

     n(동전마다)      SE(차이)       95% 구간의 반너비
          25    0.141421          0.277181
         100    0.070711          0.138590
         200    0.050000          0.097998
         768    0.025516          0.050009
         769    0.025499          0.049977
        1000    0.022361          0.043826

반너비 <= 0.05 에 필요한 n = 1.96^2 * 0.5 / 0.05^2 = 768.2918  →  769

두 동전이 모두 0.6 으로 치우쳐 있다면 (n = 400)
  이표본 검정 (H0: p1 = p2):  z = 0.0000, p = 1.0000  → 기각 못 함
  일표본 검정 (H0: p1 = 1/2): z = 4.0000, p = 6.334e-05  → 강하게 기각

n_2000, n_2015 = 400, 600
positive_2000, positive_2015 = 132, 228
p_hat_2000, p_hat_2015 = positive_2000 / n_2000, positive_2015 / n_2015

p_pooled = (positive_2000 + positive_2015) / (n_2000 + n_2015)
statistic = (p_hat_2000 - p_hat_2015) / (np.sqrt(p_pooled * (1 - p_pooled)) * np.sqrt(1/n_2000 + 1/n_2015))
# H1이 p_2000 < p_2015 이므로 왼쪽 꼬리를 센다.
p_value = stats.norm().cdf(statistic)

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

alpha = 0.05
if p_value <= alpha:
    print("Reject H_0: significant increase in myopia")
else:
    print("Fail to reject H_0")

statistic = -1.6137
p_value = 0.0533
Fail to reject H_0

d = p_hat_2000 - p_hat_2015
se = np.sqrt(p_pooled * (1 - p_pooled) * (1 / n_2000 + 1 / n_2015))
print(f"p-hat: {p_hat_2000:.4f}, {p_hat_2015:.4f}   차이 {d:+.4f}")
print(f"합동 p = {p_pooled:.4f},  SE = {se:.7f},  z = {d / se:.6f}")
print(f"단측 p = {stats.norm.cdf(d / se):.6f},  양측 p = "
      f"{2 * stats.norm.cdf(d / se):.6f}")

# 연속성 보정을 넣으면
num = abs(d) - 0.5 * (1 / n_2000 + 1 / n_2015)
print(f"연속성 보정: |d| - (1/n1+1/n2)/2 = {num:.6f},  z = {num / se:.6f},  "
      f"단측 p = {stats.norm.sf(num / se):.6f}")

# 비율을 고정한 채 표본만 k 배로 늘리면
k_need = (stats.norm.ppf(0.95) / abs(d / se)) ** 2
print(f"\n단측 0.05 를 넘으려면 표본을 k = {k_need:.4f} 배로 늘려야 한다")
print(f"{'배율':>6}{'n2000':>7}{'n2015':>7}{'z':>11}{'단측 p':>11}")
for f in (1.00, 1.04, 1.50, 2.00):
    m1, m2 = round(n_2000 * f), round(n_2015 * f)
    P = (p_hat_2000 * m1 + p_hat_2015 * m2) / (m1 + m2)
    s = np.sqrt(P * (1 - P) * (1 / m1 + 1 / m2))
    print(f"{f:>6.2f}{m1:>7}{m2:>7}{d / s:>11.6f}{stats.norm.cdf(d / s):>11.6f}")

# 신뢰구간 (합동하지 않은 표준오차)
se_u = np.sqrt(p_hat_2000 * (1 - p_hat_2000) / n_2000
               + p_hat_2015 * (1 - p_hat_2015) / n_2015)
print(f"\n비합동 SE = {se_u:.7f},  95% CI = "
      f"({d - 1.959964 * se_u:+.6f}, {d + 1.959964 * se_u:+.6f})")

p-hat: 0.3300, 0.3800   차이 -0.0500
합동 p = 0.3600,  SE = 0.0309839,  z = -1.613743
단측 p = 0.053292,  양측 p = 0.106583
연속성 보정: |d| - (1/n1+1/n2)/2 = 0.047917,  z = 1.546504,  단측 p = 0.060991

단측 0.05 를 넘으려면 표본을 k = 1.0389 배로 늘려야 한다
    배율  n2000  n2015          z       단측 p
  1.00    400    600  -1.613743   0.053292
  1.04    416    624  -1.645701   0.049913
  1.50    600    900  -1.976424   0.024053
  2.00    800   1200  -2.282177   0.011239

비합동 SE = 0.0307476,  95% CI = (-0.110264, +0.010264)

# 두 비율이 같다는 귀무가설 아래에서는 둘을 합쳐 하나의 비율로 보는 것이
# 맞다. 아래 p_pooled 가 그것이며, 표준오차를 이 값으로 만든다.
positive_male, positive_female = 24, 14
n_male, n_female = 259, 241
p_hat_male, p_hat_female = positive_male / n_male, positive_female / n_female

p_pooled = (positive_male + positive_female) / (n_male + n_female)
statistic = (p_hat_male - p_hat_female) / (np.sqrt(p_pooled * (1 - p_pooled)) * np.sqrt(1/n_male + 1/n_female))
p_value = stats.norm().sf(statistic)      # H1: p_male > p_female 이므로 오른쪽 꼬리

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

statistic = 1.4577
p_value = 0.0725

a, c = positive_male, positive_female
RR = p_hat_male / p_hat_female
log_rr = np.log(RR)
var_rare = 1 / a + 1 / c                      # 희귀사건 근사
var_exact = 1 / a - 1 / n_male + 1 / c - 1 / n_female
print(f"이환율 {p_hat_male:.6f} 대 {p_hat_female:.6f},  위험비 RR = {RR:.6f}")
print(f"log RR = {log_rr:.6f}")
print(f"{'분산식':>22}{'분산':>12}{'SE':>10}{'z':>10}{'단측 p':>10}")
for lab, v in (("1/a + 1/c (희귀근사)", var_rare),
               ("1/a-1/n1+1/c-1/n2", var_exact)):
    print(f"{lab:>22}{v:>12.6f}{np.sqrt(v):>10.6f}"
          f"{log_rr / np.sqrt(v):>10.6f}{stats.norm.sf(log_rr / np.sqrt(v)):>10.6f}")
print(f"{'비율 차 검정 (보기)':>22}{'':>12}{'':>10}{statistic:>10.6f}"
      f"{p_value:>10.6f}")

print(f"\n전체 개체 수는 {n_male + n_female}마리인데 사건은 {a + c}건뿐이다.")
print(f"  1/a + 1/c = 1/{a} + 1/{c} = {var_rare:.6f}")
print(f"  개체 수를 두 배로 늘려도 사건 수가 그대로면 분산은 그대로다.")

need_se = log_rr / stats.norm.ppf(0.95)
scale = var_rare / need_se ** 2
print(f"\n단측 0.05 를 넘기려면 SE <= {need_se:.6f},  곧 1/a+1/c <= "
      f"{need_se**2:.6f}")
print(f"  사건 수를 {scale:.4f} 배로 늘려야 한다  →  "
      f"{a * scale:.1f}건 대 {c * scale:.1f}건")

이환율 0.092664 대 0.058091,  위험비 RR = 1.595146
log RR = 0.466965
                   분산식          분산        SE         z      단측 p
      1/a + 1/c (희귀근사)    0.113095  0.336296  1.388553  0.082484
     1/a-1/n1+1/c-1/n2    0.105085  0.324168  1.440504  0.074862
          비율 차 검정 (보기)                        1.457690  0.072463

전체 개체 수는 500마리인데 사건은 38건뿐이다.
  1/a + 1/c = 1/24 + 1/14 = 0.113095
  개체 수를 두 배로 늘려도 사건 수가 그대로면 분산은 그대로다.

단측 0.05 를 넘기려면 SE <= 0.283895,  곧 1/a+1/c <= 0.080596
  사건 수를 1.4032 배로 늘려야 한다  →  33.7건 대 19.6건

import numpy as np
from scipy import stats

lo, hi = -0.04, 0.14
center, half = (lo + hi) / 2, (hi - lo) / 2
se_ci = half / stats.norm.ppf(0.975)
print(f"구간 ({lo}, {hi})  →  추정값 {center:.4f},  반너비 {half:.4f}")
print(f"비합동 SE = {half:.4f}/1.959964 = {se_ci:.6f}")
print(f"그 SE 로 읽은 z = {center / se_ci:.6f},  양측 p = "
      f"{2 * stats.norm.sf(abs(center / se_ci)):.6f}")

# 분모가 다르면 쌍대성이 깨질 수 있다
k1, n1, k2, n2 = 499, 500, 47, 50
p1, p2 = k1 / n1, k2 / n2
d = p1 - p2
P = (k1 + k2) / (n1 + n2)
se_pool = np.sqrt(P * (1 - P) * (1 / n1 + 1 / n2))
se_wald = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
z95 = stats.norm.ppf(0.975)
print(f"\n쌍대성이 깨지는 보기: {k1}/{n1} 대 {k2}/{n2}")
print(f"  p-hat {p1:.4f} 대 {p2:.4f},  차이 {d:+.4f}")
print(f"  합동 SE  {se_pool:.6f}  →  z = {d / se_pool:.4f},  "
      f"양측 p = {2 * stats.norm.sf(d / se_pool):.3e}  → 기각")
print(f"  Wald SE  {se_wald:.6f}  →  95% CI = "
      f"({d - z95 * se_wald:+.6f}, {d + z95 * se_wald:+.6f})  → 0 을 담는다")
print(f"  두 SE 의 비 = {se_wald / se_pool:.4f}")


def wilson(k, n, z=z95):
    p = k / n
    den = 1 + z * z / n
    c = (p + z * z / (2 * n)) / den
    h = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / den
    return c - h, c + h


l1, u1 = wilson(k1, n1)
l2, u2 = wilson(k2, n2)
nc_lo = d - np.sqrt((p1 - l1) ** 2 + (u2 - p2) ** 2)
nc_hi = d + np.sqrt((u1 - p1) ** 2 + (p2 - l2) ** 2)
print(f"\n  윌슨 구간: p1 ({l1:.4f}, {u1:.4f}),  p2 ({l2:.4f}, {u2:.4f})")
print(f"  뉴컴 구간: ({nc_lo:+.6f}, {nc_hi:+.6f})  → 0 을 담지 않는다")
print(f"  피셔 정확검정 양측 p = "
      f"{stats.fisher_exact([[k1, n1 - k1], [k2, n2 - k2]]).pvalue:.6f}")

구간 (-0.04, 0.14)  →  추정값 0.0500,  반너비 0.0900
비합동 SE = 0.0900/1.959964 = 0.045919
그 SE 로 읽은 z = 1.088869,  양측 p = 0.276212

쌍대성이 깨지는 보기: 499/500 대 47/50
  p-hat 0.9980 대 0.9400,  차이 +0.0580
  합동 SE  0.012603  →  z = 4.6021,  양측 p = 4.183e-06  → 기각
  Wald SE  0.033645  →  95% CI = (-0.007943, +0.123943)  → 0 을 담는다
  두 SE 의 비 = 2.6696

  윌슨 구간: p1 (0.9888, 0.9996),  p2 (0.8378, 0.9794)
  뉴컴 구간: (+0.017545, +0.160184)  → 0 을 담지 않는다
  피셔 정확검정 양측 p = 0.002660

import numpy as np
from scipy import stats

A = np.array([4.2, 5.1, 3.8, 6.0, 4.7, 5.5, 3.9, 4.4, 5.8,
              4.1, 6.3, 4.9, 5.2, 3.6, 4.6, 12.4, 5.0, 4.3])
B = np.array([5.6, 6.2, 5.9, 7.1, 6.5, 5.4, 6.8, 5.1, 6.0,
              7.4, 5.7, 6.1, 6.6, 5.3, 6.9, 5.8, 6.4, 7.0])
n1, n2 = len(A), len(B)

print(f"A: n={n1}  평균 {A.mean():.4f}  SD {A.std(ddof=1):.4f}  "
      f"중앙값 {np.median(A):.4f}")
print(f"B: n={n2}  평균 {B.mean():.4f}  SD {B.std(ddof=1):.4f}  "
      f"중앙값 {np.median(B):.4f}\n")

z = (A.mean() - B.mean()) / np.sqrt(A.var(ddof=1) / n1 + B.var(ddof=1) / n2)
print(f"z 검정       z={z:8.4f}                 p={2 * stats.norm.sf(abs(z)):.4f}")

r = stats.ttest_ind(A, B, equal_var=True)
print(f"합동 t       t={r.statistic:8.4f}  df={n1 + n2 - 2:6d}   p={r.pvalue:.4f}")

r = stats.ttest_ind(A, B, equal_var=False)
print(f"Welch t      t={r.statistic:8.4f}  df={r.df:6.2f}   p={r.pvalue:.4f}")

u = stats.mannwhitneyu(A, B)
print(f"만·휘트니    U={u.statistic:8.1f}                 p={u.pvalue:.4f}")

res = stats.permutation_test(
    (A, B), lambda x, y, axis=0: x.mean(axis) - y.mean(axis),
    permutation_type='independent', n_resamples=99_999,
    random_state=1, alternative='two-sided')
print(f"순열검정                                p={res.pvalue:.4f}")

A: n=18  평균 5.2111  SD 1.9517  중앙값 4.8000
B: n=18  평균 6.2111  SD 0.6685  중앙값 6.1500

z 검정       z= -2.0565                 p=0.0397
합동 t       t= -2.0565  df=    34   p=0.0475
Welch t      t= -2.0565  df= 20.93   p=0.0524
만·휘트니    U=    45.5                 p=0.0002
순열검정                                p=0.0337

A2 = A[A < 10]                                   # 12.4 를 빼면
print(f"A(이상점 제외): n={len(A2)}  평균 {A2.mean():.4f}  "
      f"SD {A2.std(ddof=1):.4f}")
print(f"  Welch p = {stats.ttest_ind(A2, B, equal_var=False).pvalue:.6f}")

A(이상점 제외): n=17  평균 4.7882  SD 0.7921
  Welch p = 0.000003

import numpy as np
from scipy import stats

x = np.array([3.1, 4.5, 2.8, 5.0, 3.9, 4.2])
y = np.array([5.5, 6.1, 4.8, 7.0, 5.9, 6.3])
for method in ["exact", "asymptotic"]:
    r = stats.mannwhitneyu(x, y, method=method)
    print(f"{method:12s} U={r.statistic:5.1f}  p={r.pvalue:.6f}")

print("\n동점이 있는 자료 (정수로 반올림)")
xt = np.array([3, 4, 3, 5, 4, 4])
yt = np.array([5, 6, 5, 7, 6, 6])
r = stats.mannwhitneyu(xt, yt, method="asymptotic")
print(f"  정규근사(동점 보정)  U={r.statistic:5.1f}  p={r.pvalue:.6f}")
r = stats.mannwhitneyu(xt, yt, method="exact")
print(f"  정확법(동점 무시)    U={r.statistic:5.1f}  p={r.pvalue:.6f}")

print("\n두 집단이 완전히 분리됐을 때 도달 가능한 최소 양측 p")
for k in [3, 4, 5, 6, 8]:
    a, b = np.arange(k), np.arange(k) + 100
    print(f"  n1=n2={k}: p = {stats.mannwhitneyu(a, b, method='exact').pvalue:.6f}")

exact        U=  1.0  p=0.004329
asymptotic   U=  1.0  p=0.008239

동점이 있는 자료 (정수로 반올림)
  정규근사(동점 보정)  U=  1.0  p=0.006845
  정확법(동점 무시)    U=  1.0  p=0.004329

두 집단이 완전히 분리됐을 때 도달 가능한 최소 양측 p
  n1=n2=3: p = 0.100000
  n1=n2=4: p = 0.028571
  n1=n2=5: p = 0.007937
  n1=n2=6: p = 0.002165
  n1=n2=8: p = 0.000155

import numpy as np
from scipy import stats

A = np.array([4.2, 5.1, 3.8, 6.0, 4.7, 5.5, 3.9, 4.4, 5.8,
              4.1, 6.3, 4.9, 5.2, 3.6, 4.6, 12.4, 5.0, 4.3])
B = np.array([5.6, 6.2, 5.9, 7.1, 6.5, 5.4, 6.8, 5.1, 6.0,
              7.4, 5.7, 6.1, 6.6, 5.3, 6.9, 5.8, 6.4, 7.0])

def perm_p(x, y, stat, B_=19_999, seed=0):
    """같은 순열 틀에 통계량만 갈아 끼운다."""
    rng = np.random.default_rng(seed)
    pooled = np.r_[x, y]
    n1 = len(x)
    obs = abs(stat(pooled[:n1], pooled[n1:]))
    cnt = 1
    for _ in range(B_):
        z = rng.permutation(pooled)
        cnt += abs(stat(z[:n1], z[n1:])) >= obs
    return cnt / (B_ + 1)

stats_to_try = {
    "평균 차이":      lambda a, b: a.mean() - b.mean(),
    "합동 t 통계량":  lambda a, b: stats.ttest_ind(a, b, equal_var=True).statistic,
    "Welch t 통계량": lambda a, b: stats.ttest_ind(a, b, equal_var=False).statistic,
    # 순위합은 양수이므로 abs() 가 통하도록 기댓값을 빼서 중심화한다
    "순위합(만·휘트니)": lambda a, b: (stats.rankdata(np.r_[a, b])[:len(a)].sum()
                                    - len(a) * (len(a) + len(b) + 1) / 2),
    "중앙값 차이":    lambda a, b: np.median(a) - np.median(b),
    "20% 절사평균 차이": lambda a, b: (stats.trim_mean(a, 0.2)
                                      - stats.trim_mean(b, 0.2)),
}
for name, f in stats_to_try.items():
    print(f"{name:20s} 순열 p = {perm_p(A, B, f):.4f}")

print(f"\n참고: 합동 t 검정      p = "
      f"{stats.ttest_ind(A, B, equal_var=True).pvalue:.4f}")
print(f"      만·휘트니 검정   p = {stats.mannwhitneyu(A, B).pvalue:.4f}")