일표본 검정¶

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def plot_z_statistic(statistic, ax, alternative='two-sided'):
    """검정통계량 위치에서 잘린 꼬리를 칠해 p-값을 눈으로 보여준다.

    칠해진 넓이가 곧 p-값이다. 양측이면 좌우 두 조각의 합이다.
    """
    x = np.linspace(-4, 4, 100)
    y = stats.norm().pdf(x)
    ax.plot(x, y, '-k')

    if alternative == 'less':
        x_fill = np.linspace(-4, statistic, 100)
        y_fill = stats.norm().pdf(x_fill)
        ax.fill_between(x_fill, y_fill, color='r', alpha=0.2)
    elif alternative == 'greater':
        x_fill = np.linspace(statistic, 4, 100)
        y_fill = stats.norm().pdf(x_fill)
        ax.fill_between(x_fill, y_fill, color='r', alpha=0.2)
    elif alternative == 'two-sided':
        x_fill_left = np.linspace(-4, -abs(statistic), 100)
        y_fill_left = stats.norm().pdf(x_fill_left)
        ax.fill_between(x_fill_left, y_fill_left, color='r', alpha=0.2)
        x_fill_right = np.linspace(abs(statistic), 4, 100)
        y_fill_right = stats.norm().pdf(x_fill_right)
        ax.fill_between(x_fill_right, y_fill_right, color='r', alpha=0.2)

    ax.spines['left'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.set_yticks([])

mu = 50
n = 500
x_bar = 48
s = 20.3

# n = 500 >= 30 이므로 sigma를 몰라도 s를 넣고 z를 쓴다.
statistic = (x_bar - mu) / (s / np.sqrt(n))
# sf(x) = 1 - cdf(x) 이고 꼬리확률에서는 sf가 수치적으로 더 정확하다.
# abs를 씌우고 2를 곱해 양측으로 만든다.
p_value = stats.norm().sf(abs(statistic)) * 2

print(f"Statistic: {statistic:.4f}")
print(f"P-value  : {p_value:.4f}")

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

fig, ax = plt.subplots(figsize=(12, 3))
plot_z_statistic(statistic, ax=ax, alternative='two-sided')
plt.show()

Statistic: -2.2030
P-value  : 0.0276
Reject H0 (Choose H1)

mu = 50
n = 500
x_bar = 48
s = 20.3

statistic = (x_bar - mu) / (s / np.sqrt(n))
# H1이 "작다" 쪽이므로 왼쪽 꼬리만 센다. abs도, 2를 곱하는 것도 없다.
p_value = stats.norm().cdf(statistic)

print(f"Statistic : {statistic:.4f}")
print(f"P-value   : {p_value:.4f}")

alpha = 0.05
if p_value <= alpha:
    print("We choose H1, or using statistician's jargon, reject H0")
else:
    print("We choose H0, or using statistician's jargon, fail to reject H0")

fig, ax = plt.subplots(figsize=(12, 3))
plot_z_statistic(statistic, ax=ax, alternative='less')
plt.show()

Statistic : -2.2030
P-value   : 0.0138
We choose H1, or using statistician's jargon, reject H0

mu = 50
n = 500
x_bar = 52
s = 20.3

statistic = (x_bar - mu) / (s / np.sqrt(n))
# 이번에는 H1이 "크다" 쪽이므로 오른쪽 꼬리만 센다.
p_value = stats.norm().sf(statistic)

print(f"Statistic : {statistic:.4f}")
print(f"P-value   : {p_value:.4f}")

alpha = 0.05
if p_value <= alpha:
    print("We choose H1, or using statistician's jargon, reject H0")
else:
    print("We choose H0, or using statistician's jargon, fail to reject H0")

fig, ax = plt.subplots(figsize=(12, 3))
plot_z_statistic(statistic, ax=ax, alternative='greater')
plt.show()

Statistic : 2.2030
P-value   : 0.0138
We choose H1, or using statistician's jargon, reject H0

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

def plot_t_statistic(statistic, df, ax, alternative='two-sided'):
    """t 분포를 그리고 p-값에 해당하는 꼬리를 칠한다.

    대립가설의 방향에 따라 칠하는 쪽이 달라진다. 양측이면 양쪽을 모두
    칠하므로 넓이가 두 배가 되고, 그래서 양측 p-값이 단측의 두 배다.
    """
    x = np.linspace(-4, 4, 100)
    y = stats.t(df).pdf(x)
    ax.plot(x, y, '-k')

    if alternative == 'less':
        x_fill = np.linspace(-4, statistic, 100)
        y_fill = stats.t(df).pdf(x_fill)
        ax.fill_between(x_fill, y_fill, color='k', alpha=0.2)
    elif alternative == 'greater':
        x_fill = np.linspace(statistic, 4, 100)
        y_fill = stats.t(df).pdf(x_fill)
        ax.fill_between(x_fill, y_fill, color='k', alpha=0.2)
    elif alternative == 'two-sided':
        x_fill_left = np.linspace(-4, -abs(statistic), 100)
        y_fill_left = stats.t(df).pdf(x_fill_left)
        ax.fill_between(x_fill_left, y_fill_left, color='k', alpha=0.2)
        x_fill_right = np.linspace(abs(statistic), 4, 100)
        y_fill_right = stats.t(df).pdf(x_fill_right)
        ax.fill_between(x_fill_right, y_fill_right, color='k', alpha=0.2)

    ax.spines['left'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.set_yticks(())

# sigma 를 모르므로 s 를 쓴다. 그래서 z 가 아니라 t 로 간다.
mu_0 = 5
x_bar = 4
s = 2
n = 25

statistic = (x_bar - mu_0) / (s / np.sqrt(n))
df = n - 1                      # s를 자료에서 추정했으므로 하나를 잃는다
p_value = stats.t(df).cdf(statistic)

print(f"T-statistic: {statistic:.4f}")
print(f"P-value    : {p_value:.4f}")

fig, ax = plt.subplots(figsize=(12, 3))
plot_t_statistic(statistic, df=df, ax=ax, alternative='less')
plt.show()

T-statistic: -2.5000
P-value    : 0.0098

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
M, n = 400_000, 25
c_left = stats.t(n - 1).ppf(0.05)        # 좌측 임계값 (우측은 대칭이라 -c_left)
c_two = stats.t(n - 1).ppf(0.975)        # 양측 임계값

print(f"n = {n},  되풀이 {M:,}회,  명목 수준 0.05")
print(f"{'모집단':>10s} {'왜도':>7s} {'좌측':>8s} {'우측':>8s} {'양측':>8s}")
for label, draw, mu, g1 in [
    ("정규", lambda s: rng.normal(0, 1, s), 0.0, 0.0),
    ("지수", lambda s: rng.exponential(1, s), 1.0, 2.0),
    ("로그정규", lambda s: rng.lognormal(0, 1, s), np.exp(0.5), 6.1849),
]:
    x = draw((M, n))
    # 참 평균 mu 를 귀무가설값으로 넣는다. 곧 H0 가 참인 세상이다.
    t = (x.mean(1) - mu) / (x.std(1, ddof=1) / np.sqrt(n))
    print(f"{label:>10s} {g1:7.2f} {(t < c_left).mean():8.4f} "
          f"{(t > -c_left).mean():8.4f} {(np.abs(t) > c_two).mean():8.4f}")
print(f"\n몬테카를로 오차 = sqrt(0.05*0.95/{M}) "
      f"= {np.sqrt(0.05 * 0.95 / M):.5f}")

n = 25,  되풀이 400,000회,  명목 수준 0.05
       모집단      왜도       좌측       우측       양측
        정규    0.00   0.0500   0.0501   0.0502
        지수    2.00   0.1033   0.0214   0.0769
      로그정규    6.18   0.1626   0.0089   0.1245

몬테카를로 오차 = sqrt(0.05*0.95/400000) = 0.00034

n = 7
df = n - 1
# 원자료 없이 t 통계량만 있어도 p-값을 구할 수 있다.
# 필요한 것은 통계량과 자유도뿐이다.
statistic = -1.9
p_value = stats.t(df).cdf(statistic)
print(f"{statistic = :.4f}")
print(f"{p_value = :.4f}")

fig, ax = plt.subplots(figsize=(12, 3))
plot_t_statistic(statistic, df=df, ax=ax, alternative='less')
plt.show()

statistic = -1.9000
p_value = 0.0531

import numpy as np
from scipy import stats

statistic = -1.9
print(f"{'df':>6s} {'임계값':>10s} {'p-값':>9s}")
for df in [6, 10, 20, 30, 60, 120, 1000]:
    print(f"{df:6d} {stats.t(df).ppf(0.05):10.4f} {stats.t(df).cdf(statistic):9.4f}")
# 자유도를 무한히 키운 극한이 표준정규다.
print(f"{'∞':>6s} {stats.norm.ppf(0.05):10.4f} {stats.norm.cdf(statistic):9.4f}")

    df        임계값       p-값
     6    -1.9432    0.0531
    10    -1.8125    0.0433
    20    -1.7247    0.0360
    30    -1.6973    0.0335
    60    -1.6706    0.0311
   120    -1.6577    0.0299
  1000    -1.6464    0.0289
     ∞    -1.6449    0.0287

n = 6
df = n - 1
statistic = 2.75
# 양측이므로 한쪽 꼬리를 두 배 한다.
p_value = stats.t(df).sf(statistic) * 2
print(f"{statistic = :.4f}")
print(f"{p_value = :.4f}")

fig, ax = plt.subplots(figsize=(12, 3))
plot_t_statistic(statistic, df=df, ax=ax, alternative='two-sided')
plt.show()

statistic = 2.7500
p_value = 0.0403

import numpy as np
from scipy import stats

mu_0, x_bar, s, n = 530, 528, 4, 20
se = s / np.sqrt(n)
df = n - 1
statistic = (x_bar - mu_0) / se
p_value = 2 * stats.t(df).cdf(statistic)
t_crit = stats.t(df).ppf(0.975)
lo, hi = x_bar - t_crit * se, x_bar + t_crit * se

print(f"표준오차 = {se:.6f}")
print(f"t = {statistic:.5f}   (-sqrt(5) = {-np.sqrt(5):.5f})")
print(f"양측 p = {p_value:.5f}")
print(f"임계값 t(0.025, {df}) = {t_crit:.5f}")
print(f"95% 신뢰구간 = ({lo:.4f}, {hi:.4f})")
print(f"구간이 530 을 담는가: {lo <= mu_0 <= hi}")
print(f"부족량 530 - mu 의 95% 구간 = "
      f"({mu_0 - hi:.4f}, {mu_0 - lo:.4f}) mL")
print(f"표시량에 대한 비율 = ({(mu_0 - hi) / mu_0 * 100:.3f}%, "
      f"{(mu_0 - lo) / mu_0 * 100:.3f}%)")

표준오차 = 0.894427
t = -2.23607   (-sqrt(5) = -2.23607)
양측 p = 0.03754
임계값 t(0.025, 19) = 2.09302
95% 신뢰구간 = (526.1279, 529.8721)
구간이 530 을 담는가: False
부족량 530 - mu 의 95% 구간 = (0.1279, 3.8721) mL
표시량에 대한 비율 = (0.024%, 0.731%)

samples = np.array([78, 83, 68, 72, 88])

n = samples.shape[0]
df = n - 1
x_bar = samples.mean()
s = samples.std(ddof=1)
mu = 70

confidence_level = 0.95
alpha = 1 - confidence_level
t_score = (x_bar - mu) / (s / np.sqrt(n))
# H1이 mu > 70인 단측이므로 오른쪽 꼬리만 센다.
# 여기에 양측 공식 sf(abs(t))*2 를 쓰면 p-값이 두 배가 되어
# 이 자료에서는 결론이 뒤집힌다(0.0484 대 0.0969).
p_value = stats.t(df=df).sf(t_score)

print(f"Test statistic (t-score): {t_score:.4f}")
print(f"p-value                 : {p_value:.4f}")

if p_value <= alpha:
    print("Reject H_0: Sufficient evidence to support the alternative hypothesis.")
else:
    print("Fail to reject H_0: Insufficient evidence to support the alternative hypothesis.")

fig, ax = plt.subplots(figsize=(12, 3))
plot_t_statistic(t_score, df=df, ax=ax, alternative='greater')
ax.legend(["t-distribution", f"t statistic = {t_score:.4f}"])
plt.show()

Test statistic (t-score): 2.1600
p-value                 : 0.0484
Reject H_0: Sufficient evidence to support the alternative hypothesis.

# 표시량 128 에 못 미치는지를 묻는 단측검정이다. 그래서 아래에서 cdf 를 쓴다.
mu_0 = 128
x_bar = 127.2
s = 2.1
n = 12
statistic = (x_bar - mu_0) / (s / np.sqrt(n))
df = n - 1
p_value = stats.t(df).cdf(statistic)

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

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

fig, ax = plt.subplots(figsize=(12, 3))
plot_t_statistic(statistic, df=df, ax=ax, alternative='less')
plt.show()

statistic = -1.3197
p_value = 0.1069
Fail to reject H_0

import numpy as np
from scipy import stats

mu_0, mu_true, sigma, n = 128, 127.2, 2.1, 12
df = n - 1
ncp = (mu_true - mu_0) / (sigma / np.sqrt(n))     # 비중심 모수
crit = stats.t(df).ppf(0.05)                      # 좌측 임계값

print(f"비중심 모수 delta = {ncp:.6f}")
print(f"좌측 임계값       = {crit:.6f}")
print(f"비중심 t 로 구한 정확한 검정력 = "
      f"{stats.nct(df, ncp).cdf(crit):.4f}")
print(f"정규근사 (z 임계값을 쓴 어림)  = "
      f"{stats.norm.cdf(-stats.norm.ppf(0.95) - ncp):.4f}")
print(f"정규근사 (t 임계값을 쓴 어림)  = "
      f"{stats.norm.cdf(crit - ncp):.4f}")

비중심 모수 delta = -1.319658
좌측 임계값       = -1.795885
비중심 t 로 구한 정확한 검정력 = 0.3422
정규근사 (z 임계값을 쓴 어림)  = 0.3725
정규근사 (t 임계값을 쓴 어림)  = 0.3170

import numpy as np
from scipy import stats

p_0, p_1, alpha, beta = 0.49, 0.55, 0.05, 0.20
z_a = stats.norm.ppf(1 - alpha / 2)
z_b = stats.norm.ppf(1 - beta)

n_formula = ((z_a * np.sqrt(p_0 * (1 - p_0))
              + z_b * np.sqrt(p_1 * (1 - p_1))) / (p_1 - p_0)) ** 2
print(f"정규근사 공식이 주는 n = {n_formula:.2f}  ->  {int(np.ceil(n_formula))}")


def exact_power(n, p_true, p_0=0.49, alpha=0.05):
    """전수 열거로 구한 양측 점수검정의 정확한 기각확률."""
    k = np.arange(n + 1)
    z = (k / n - p_0) / np.sqrt(p_0 * (1 - p_0) / n)
    reject = np.abs(z) > stats.norm.ppf(1 - alpha / 2)
    return stats.binom.pmf(k, n, p_true)[reject].sum()


print(f"그 n 에서의 정확 검정력 = {exact_power(544, p_1):.4f}  (목표 0.80)")
print(f"그 n 에서의 정확 크기   = {exact_power(544, p_0):.4f}  (명목 0.05)")

# 격자에서 "처음 0.80 을 넘는 n" 을 집으면 안 된다. 검정력이 n 의
# 증가함수가 아니기 때문이다. 격자 전체를 보고 되돌아오는 칸을 센다.
ns = np.arange(300, 901)
pw = np.array([exact_power(n, p_1) for n in ns])
first = ns[pw >= 0.80].min()
dips = ns[(ns > first) & (pw < 0.80)]
print(f"\n처음 0.80 에 닿는 n = {first}  (검정력 {exact_power(first, p_1):.4f})")
print(f"그보다 큰데도 0.80 에 못 미치는 n = {list(dips)}")
print(f"0.80 이 무너지지 않는 가장 작은 n = {dips.max() + 1}")

정규근사 공식이 주는 n = 543.27  ->  544
그 n 에서의 정확 검정력 = 0.7985  (목표 0.80)
그 n 에서의 정확 크기   = 0.0484  (명목 0.05)

처음 0.80 에 닿는 n = 535  (검정력 0.8017)
그보다 큰데도 0.80 에 못 미치는 n = [536, 538, 540, 542, 544, 553, 555, 557]
0.80 이 무너지지 않는 가장 작은 n = 558

import numpy as np
from scipy import stats

n, k, p_0 = 1000, 920, 0.90
sd = np.sqrt(n * p_0 * (1 - p_0))          # X 의 표준편차
se_0 = np.sqrt(p_0 * (1 - p_0) / n)        # p_hat 의 표준오차

z = (k / n - p_0) / se_0
exact = stats.binomtest(k, n, p_0, alternative="greater").pvalue
z_cc = (k - 0.5 - n * p_0) / sd            # 연속성 보정
g1 = (1 - 2 * p_0) / sd                    # 이항분포의 왜도
# 연속성 보정 위에 왜도 항을 하나 더 얹은 에지워스 근사
edge = stats.norm.sf(z_cc) + (g1 / 6) * (z_cc**2 - 1) * stats.norm.pdf(z_cc)

print(f"z (점수)  = {z:.6f}")
print(f"z (연속성 보정) = {z_cc:.6f}")
print(f"이항분포의 왜도 = {g1:.6f}\n")
for label, value in [("정확 이항", exact),
                     ("정규근사", stats.norm.sf(z)),
                     ("정규근사 + 연속성 보정", stats.norm.sf(z_cc)),
                     ("정규근사 + 연속성 + 왜도", edge)]:
    print(f"{label:24s} p = {value:.6f}   오차 {abs(value - exact):.6f}")

z (점수)  = 2.108185
z (연속성 보정) = 2.055480
이항분포의 왜도 = -0.084327

정확 이항                    p = 0.017612   오차 0.000000
정규근사                     p = 0.017507   오차 0.000104
정규근사 + 연속성 보정            p = 0.019916   오차 0.002305
정규근사 + 연속성 + 왜도          p = 0.017730   오차 0.000118

p_hat = 22 / 200
p = 0.08
n = 200

# 표준오차의 분모에 p_hat이 아니라 가설값 p를 넣는다.
# H0가 참이라는 가정 아래의 확률을 재는 것이 검정이기 때문이다.
statistic = (p_hat - p) / np.sqrt(p * (1 - p) / n)
p_value = stats.norm().sf(abs(statistic)) * 2

print(f"Statistic: {statistic:.4f}")
print(f"P-value : {p_value:.4f}")

Statistic: 1.5639
P-value : 0.1179

from scipy.stats import binomtest

ph, p0, n = 22 / 200, 0.08, 200
se_score = np.sqrt(p0 * (1 - p0) / n)          # H0 의 p0 를 쓴다 (점수검정)
se_wald = np.sqrt(ph * (1 - ph) / n)           # 관측된 p-hat 을 쓴다 (발트)
print(f"p-hat = {ph},  p0 = {p0}")
print(f"SE(점수) = {se_score:.7f}   SE(발트) = {se_wald:.7f}"
      f"   비 = {se_wald / se_score:.4f}")
for lab, s in (("점수 (p0)", se_score), ("발트 (p-hat)", se_wald)):
    z = (ph - p0) / s
    print(f"  {lab:>12}: z = {z:.6f},  양측 p = {2 * stats.norm.sf(abs(z)):.6f}")

z_cc = (abs(ph - p0) - 0.5 / n) / se_score
print(f"  {'연속성 보정':>12}: z = {z_cc:.6f},  양측 p = "
      f"{2 * stats.norm.sf(z_cc):.6f}")
print(f"  {'정확 이항':>12}: p = {binomtest(22, 200, p0).pvalue:.6f}")

p-hat = 0.11,  p0 = 0.08
SE(점수) = 0.0191833   SE(발트) = 0.0221246   비 = 1.1533
       점수 (p0): z = 1.563858,  양측 p = 0.117851
    발트 (p-hat): z = 1.355954,  양측 p = 0.175114
        연속성 보정: z = 1.433537,  양측 p = 0.151705
         정확 이항: p = 0.117778

# 비율의 검정에서 표준오차는 관측된 p_hat 이 아니라 귀무가설의 p 로 만든다.
# 귀무가설이 참이라는 전제 아래에서의 흩어짐을 재야 하기 때문이다.
p_hat = 40 / 120
p = 0.26
n = 120

statistic = (p_hat - p) / np.sqrt(p * (1 - p) / n)
p_value = stats.norm().sf(statistic)      # 단측이므로 오른쪽 꼬리만

print(f"Statistic: {statistic:.4f}")
print(f"P-value: {p_value:.4f}")

Statistic: 1.8314
P-value: 0.0335

p0, n = 0.26, 120
se = np.sqrt(p0 * (1 - p0) / n)
k_thresh = (p0 + stats.norm.ppf(0.95) * se) * n
print(f"SE = {se:.7f}")
print(f"기각에 필요한 p-hat > {p0 + stats.norm.ppf(0.95) * se:.6f}"
      f"  →  k > {k_thresh:.4f}")
print(f"\n{'k':>4}{'p-hat':>9}{'z':>9}{'근사 단측 p':>13}{'정확 단측 p':>13}")
for k in (38, 39, 40, 41):
    z = (k / n - p0) / se
    print(f"{k:>4}{k / n:>9.4f}{z:>9.4f}{stats.norm.sf(z):>13.5f}"
          f"{binomtest(k, n, p0, alternative='greater').pvalue:>13.5f}")
print(f"\n관측값 40 에서 양측으로 하면 p = "
      f"{2 * stats.norm.sf((40 / n - p0) / se):.5f}")

SE = 0.0400416
기각에 필요한 p-hat > 0.325863  →  k > 39.1035

   k    p-hat        z      근사 단측 p      정확 단측 p
  38   0.3167   1.4152      0.07851      0.09665
  39   0.3250   1.6233      0.05226      0.06674
  40   0.3333   1.8314      0.03352      0.04464
  41   0.3417   2.0395      0.02070      0.02892

관측값 40 에서 양측으로 하면 p = 0.06704

k = 113
n = 200
p_0 = 0.5

# 정확검정: 이항분포에서 P(X >= 113)을 그대로 더한다. 근사가 없다.
result = stats.binomtest(k, n, p=p_0, alternative="greater")
print(f"Exact P-value: {result.pvalue:.4f}")

# 근사검정: 이항분포를 정규분포로 바꿔 계산한다.
p_hat = k / n
approx_statistic = (p_hat - p_0) / np.sqrt(p_0 * (1 - p_0) / n)
approx_p_value = stats.norm().sf(approx_statistic)
print(f"Approximate Statistic: {approx_statistic:.4f}")
print(f"Approximate P-value: {approx_p_value:.4f}")

Exact P-value: 0.0384
Approximate Statistic: 1.8385
Approximate P-value: 0.0330

k, n, p_0 = 113, 200, 0.5
p_hat = k / n
se = np.sqrt(p_0 * (1 - p_0) / n)
z_plain = (p_hat - p_0) / se
z_cc = ((k - 0.5) / n - p_0) / se          # 막대의 절반을 되돌려 준다
print(f"{'방법':>14}{'z':>11}{'단측 p':>12}{'정확값과의 차':>14}")
exact = stats.binomtest(k, n, p=p_0, alternative="greater").pvalue
print(f"{'정확 이항':>14}{'':>11}{exact:>12.6f}{0.0:>14.6f}")
for lab, z in (("정규근사", z_plain), ("연속성 보정", z_cc)):
    p = stats.norm.sf(z)
    print(f"{lab:>14}{z:>11.6f}{p:>12.6f}{p - exact:>+14.6f}")
print(f"\n보정의 크기: 0.5/n = {0.5 / n} (비율 척도),  "
      f"z 로는 {z_plain - z_cc:.6f}")
print(f"이항분포의 왜도 (1-2p0)/sqrt(n p0 q0) = "
      f"{(1 - 2 * p_0) / np.sqrt(n * p_0 * (1 - p_0)):.6f}  (p0 = 0.5 이므로 0)")

            방법          z        단측 p       정확값과의 차
         정확 이항               0.038419      0.000000
          정규근사   1.838478    0.032996     -0.005423
        연속성 보정   1.767767    0.038550     +0.000131

보정의 크기: 0.5/n = 0.0025 (비율 척도),  z 로는 0.070711
이항분포의 왜도 (1-2p0)/sqrt(n p0 q0) = 0.000000  (p0 = 0.5 이므로 0)

# n*p_0 = 13 으로 작은 편이라 정규근사가 미덥지 않다.
# 그래서 이항분포로 정확히 구한 값과 근사값을 나란히 놓고 견준다.
k = 11
n = 65
p_0 = 0.2

result = stats.binomtest(k, n, p=p_0, alternative="less")
print(f"Exact P-value: {result.pvalue:.4f}")

p_hat = k / n
approx_statistic = (p_hat - p_0) / np.sqrt(p_0 * (1 - p_0) / n)
approx_p_value = stats.norm().cdf(approx_statistic)
print(f"Approximate Statistic: {approx_statistic:.4f}")
print(f"Approximate P-value: {approx_p_value:.4f}")

Exact P-value: 0.3301
Approximate Statistic: -0.6202
Approximate P-value: 0.2676

k, n, p_0 = 11, 65, 0.2
p_hat = k / n
se = np.sqrt(p_0 * (1 - p_0) / n)
exact = stats.binomtest(k, n, p=p_0, alternative="less").pvalue
z_plain = (p_hat - p_0) / se
z_cc = ((k + 0.5) / n - p_0) / se          # 왼쪽 꼬리이므로 +0.5
print(f"{'방법':>14}{'z':>11}{'단측 p':>12}{'정확값과의 차':>14}")
print(f"{'정확 이항':>14}{'':>11}{exact:>12.6f}{0.0:>14.6f}")
for lab, z in (("정규근사", z_plain), ("연속성 보정", z_cc)):
    p = stats.norm.cdf(z)
    print(f"{lab:>14}{z:>11.6f}{p:>12.6f}{p - exact:>+14.6f}")
print(f"\nn p0 = {n * p_0}, n(1-p0) = {n * (1 - p_0)}")
print(f"이항분포의 왜도 (1-2p0)/sqrt(n p0 q0) = "
      f"{(1 - 2 * p_0) / np.sqrt(n * p_0 * (1 - p_0)):.6f}")
print("보기 14 (p0 = 0.5) 의 왜도는 0 이었다. 그 차이가 남은 어긋남의 정체다.")

            방법          z        단측 p       정확값과의 차
         정확 이항               0.330066      0.000000
          정규근사  -0.620174    0.267572     -0.062494
        연속성 보정  -0.465130    0.320919     -0.009146

n p0 = 13.0, n(1-p0) = 52.0
이항분포의 왜도 (1-2p0)/sqrt(n p0 q0) = 0.186052
보기 14 (p0 = 0.5) 의 왜도는 0 이었다. 그 차이가 남은 어긋남의 정체다.

import numpy as np
from scipy.stats import wilcoxon

task_times = np.array([16, 14, 15, 17, 13, 18, 14, 16, 15, 19])
hypothetical_median = 15

# 차이의 **크기 순위**에 부호를 붙여 더한다.
# 부호검정이 방향만 보는 것과 달리 크기 정보를 일부 살리므로 검정력이 높다.
differences = task_times - hypothetical_median
stat, p_value = wilcoxon(differences)

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

alpha = 0.05
if p_value < alpha:
    print("Reject the null hypothesis: The median is significantly different.")
else:
    print("Fail to reject the null hypothesis: No significant difference.")

Test Statistic: 10.5
P-value: 0.28702974430361805
Fail to reject the null hypothesis: No significant difference.

from scipy.stats import rankdata
from collections import Counter

d = task_times - hypothetical_median
nz = d[d != 0]
n = len(nz)
r = rankdata(np.abs(nz))
W_plus = r[nz > 0].sum()
W_minus = r[nz < 0].sum()
print(f"차이          {list(d)}")
print(f"0 을 뺀 차이  {list(nz)}   n = {n}  (버린 개수 {len(d) - n})")
print(f"|차이|의 순위 {list(r)}")
print(f"W+ = {W_plus},  W- = {W_minus},  합 = {W_plus + W_minus}"
      f"  (n(n+1)/2 = {n * (n + 1) / 2})")

E = n * (n + 1) / 4
V = n * (n + 1) * (2 * n + 1) / 24
cnt = Counter(np.abs(nz))
corr = sum(c ** 3 - c for c in cnt.values()) / 2
V_tie = (n * (n + 1) * (2 * n + 1) - corr) / 24
print(f"\nE[W] = n(n+1)/4 = {E},  Var = n(n+1)(2n+1)/24 = {V}")
print(f"동점 보정량 sum(t^3-t)/2 = {corr},  보정한 Var = {V_tie}")
for lab, v in (("보정 없음", V), ("동점 보정", V_tie)):
    z = (W_minus - E) / np.sqrt(v)
    print(f"  {lab}: z = {z:.6f},  양측 p = {2 * stats.norm.cdf(z):.11f}")
print(f"  scipy 가 준 값: {p_value:.11f}")

차이          [1, -1, 0, 2, -2, 3, -1, 1, 0, 4]
0 을 뺀 차이  [1, -1, 2, -2, 3, -1, 1, 4]   n = 8  (버린 개수 2)
|차이|의 순위 [2.5, 2.5, 5.5, 5.5, 7.0, 2.5, 2.5, 8.0]
W+ = 25.5,  W- = 10.5,  합 = 36.0  (n(n+1)/2 = 36.0)

E[W] = n(n+1)/4 = 18.0,  Var = n(n+1)(2n+1)/24 = 51.0
동점 보정량 sum(t^3-t)/2 = 33.0,  보정한 Var = 49.625
  보정 없음: z = -1.050210,  양측 p = 0.29362154393
  동점 보정: z = -1.064660,  양측 p = 0.28702974430
  scipy 가 준 값: 0.28702974430

from scipy.stats import binom
import numpy as np

def sign_test(data, median_hypothesis):
    """차이의 부호만 세는 검정. 크기는 완전히 버린다.

    H0가 "중앙값 = median_hypothesis"이면 각 관측값이 그보다 클 확률이 1/2이므로
    양수의 개수가 Bin(n, 0.5)를 따른다. 그래서 이항분포로 p-값을 구한다.
    """
    differences = np.array(data) - median_hypothesis
    n_plus = np.sum(differences > 0)
    n_minus = np.sum(differences < 0)
    ties = np.sum(differences == 0)      # 0인 것은 세지 않고 버린다
    W = min(n_plus, n_minus)
    n = n_plus + n_minus                 # 동점을 뺀 개수
    p_value = 2 * binom.cdf(W, n, 0.5)
    return p_value, n_plus, n_minus, ties

data = [9.8, 10.1, 9.9, 10.2, 10.4, 10.3, 10.0, 9.7, 10.5, 9.6, 10.0, 9.9, 10.2, 9.8, 10.1]
p_value, n_plus, n_minus, ties = sign_test(data, 10)
print(f"P-value: {p_value}, n+: {n_plus}, n-: {n_minus}, Ties: {ties}")

P-value: 1.0, n+: 7, n-: 6, Ties: 2

from scipy.stats import binom

n_pm = n_plus + n_minus
W = min(n_plus, n_minus)
print(f"n+ = {n_plus}, n- = {n_minus}, 동점 {ties}  →  유효 n = {n_pm}")
print(f"W = min = {W}")
print(f"binom.cdf({W}, {n_pm}, 0.5) = {binom.cdf(W, n_pm, 0.5):.12f}")
print(f"  n 이 홀수이면 cdf((n-1)/2) 가 정확히 1/2 이다 → p = "
      f"{2 * binom.cdf(W, n_pm, 0.5):.12f}")

print("\n'작은 쪽 꼬리의 두 배' 관례가 1 을 넘는 경우")
for m, w in ((13, 6), (14, 7), (20, 10)):
    print(f"  n = {m}, W = {w}:  2*cdf = {2 * binom.cdf(w, m, 0.5):.6f}"
          f"{'   ← 1 을 넘는다' if 2 * binom.cdf(w, m, 0.5) > 1 else ''}")

print("\n부호검정의 실제 수준 (이산성 때문에 명목 0.05 를 채우지 못한다)")
print(f"{'n':>4}{'기각 임계 k':>12}{'실제 수준':>12}")
for m in (10, 13, 14, 15, 20, 31, 32, 40):
    kk, lev = None, 0.0
    for k in range(0, m // 2 + 1):
        if 2 * binom.cdf(k, m, 0.5) < 0.05:
            kk, lev = k, 2 * binom.cdf(k, m, 0.5)
    print(f"{m:>4}{str(kk):>12}{lev:>12.5f}")

n+ = 7, n- = 6, 동점 2  →  유효 n = 13
W = min = 6
binom.cdf(6, 13, 0.5) = 0.500000000000
  n 이 홀수이면 cdf((n-1)/2) 가 정확히 1/2 이다 → p = 1.000000000000

'작은 쪽 꼬리의 두 배' 관례가 1 을 넘는 경우
  n = 13, W = 6:  2*cdf = 1.000000
  n = 14, W = 7:  2*cdf = 1.209473   ← 1 을 넘는다
  n = 20, W = 10:  2*cdf = 1.176197   ← 1 을 넘는다

부호검정의 실제 수준 (이산성 때문에 명목 0.05 를 채우지 못한다)
   n     기각 임계 k       실제 수준
  10           1     0.02148
  13           2     0.02246
  14           2     0.01294
  15           3     0.03516
  20           5     0.04139
  31           9     0.02945
  32           9     0.02006
  40          13     0.03848

rng = np.random.default_rng(7)
B = 40_000
print("정규자료에서 t 검정과 부호검정의 검정력 (양측 0.05)")
print(f"{'n':>4}{'참 mu':>8}{'t 검정':>10}{'부호검정':>11}{'비':>8}")
for m, mu in ((20, 0.5), (40, 0.5), (20, 0.3), (40, 0.3)):
    x = rng.normal(mu, 1, (B, m))
    t = x.mean(1) / (x.std(1, ddof=1) / np.sqrt(m))
    p_t = np.mean(2 * stats.t.sf(np.abs(t), m - 1) < 0.05)
    npl = (x > 0).sum(1)
    kk = np.minimum(npl, m - npl)
    p_s = np.mean(np.minimum(1, 2 * binom.cdf(kk, m, 0.5)) < 0.05)
    print(f"{m:>4}{mu:>8.2f}{p_t:>10.4f}{p_s:>11.4f}{p_s / p_t:>8.4f}")
print(f"\n점근 상대효율 2/pi = {2 / np.pi:.6f}")
print("이것은 '같은 검정력에 필요한 표본의 비'의 극한값이고,")
print("위 표의 '검정력의 비'와는 다른 양이다.")

정규자료에서 t 검정과 부호검정의 검정력 (양측 0.05)
   n    참 mu      t 검정       부호검정       비
  20    0.50    0.5638     0.3820  0.6775
  40    0.50    0.8674     0.6611  0.7622
  20    0.30    0.2446     0.1636  0.6689
  40    0.30    0.4560     0.2872  0.6298

점근 상대효율 2/pi = 0.636620
이것은 '같은 검정력에 필요한 표본의 비'의 극한값이고,
위 표의 '검정력의 비'와는 다른 양이다.

import numpy as np

np.random.seed(42)      # 아래 출력을 재현하려면 고정한다

def bootstrap_confidence_interval(data, statistic=np.mean, n_resamples=10000, ci=95):
    """원자료에서 **복원**추출로 재표본을 만들어 통계량의 분포를 흉내 낸다.

    size=len(data)가 핵심이다. 원자료와 같은 크기로 뽑아야
    "같은 크기의 표본을 다시 얻었다면"이라는 물음에 답하는 것이 된다.
    """
    bootstrap_distribution = np.array([
        statistic(np.random.choice(data, size=len(data), replace=True))
        for _ in range(n_resamples)
    ])
    lower_bound = np.percentile(bootstrap_distribution, (100 - ci) / 2)
    upper_bound = np.percentile(bootstrap_distribution, 100 - (100 - ci) / 2)
    return lower_bound, upper_bound, bootstrap_distribution

data = [5, 7, 9, 12, 15]
lower, upper, _ = bootstrap_confidence_interval(data)
print(f"95% CI for the Mean: ({lower:.2f}, {upper:.2f})")

95% CI for the Mean: (6.60, 12.80)

x = np.array(data, dtype=float)
m = len(x)
s1 = x.std(ddof=1)
s0 = x.std(ddof=0)
tc = stats.t.ppf(0.975, m - 1)
hw_t = tc * s1 / np.sqrt(m)
print(f"자료 {list(x)}   평균 {x.mean()}")
print(f"s(ddof=1) = {s1:.6f},  s(ddof=0) = {s0:.6f}")
print(f"t 구간: ({x.mean() - hw_t:.4f}, {x.mean() + hw_t:.4f})  반너비 {hw_t:.6f}")
se_boot = s0 / np.sqrt(m)
print(f"붓스트랩 평균의 참 표준편차 = s(ddof=0)/sqrt(n) = {se_boot:.6f}")
print(f"  정규근사 반너비 1.96*SE = {1.959964 * se_boot:.6f}")
print(f"  실제 붓스트랩 구간 (6.60, 12.80) 의 반너비 = {(12.80 - 6.60) / 2:.6f}")
print(f"\n예측한 너비의 비 = (1.96/t_4)*sqrt((n-1)/n) = "
      f"{(1.959964 / tc) * np.sqrt((m - 1) / m):.6f}")
print(f"실제 비 = {(12.80 - 6.60) / 2 / hw_t:.6f}")

rng = np.random.default_rng(11)
B, nn = 20_000, 5
tc5 = stats.t.ppf(0.975, nn - 1)
cov_t = cov_b = 0
for _ in range(B):
    z = rng.normal(0, 1, nn)
    cov_t += abs(z.mean()) <= tc5 * z.std(ddof=1) / np.sqrt(nn)
    bs = z[rng.integers(0, nn, (2000, nn))].mean(1)
    lo, hi = np.percentile(bs, [2.5, 97.5])
    cov_b += lo <= 0 <= hi
print(f"\nn = 5 정규자료, 명목 95% (반복 {B:,})")
print(f"  t 구간의 포함률          {cov_t / B:.4f}")
print(f"  백분위 붓스트랩의 포함률 {cov_b / B:.4f}")

자료 [5.0, 7.0, 9.0, 12.0, 15.0]   평균 9.6
s(ddof=1) = 3.974921,  s(ddof=0) = 3.555278
t 구간: (4.6645, 14.5355)  반너비 4.935517
붓스트랩 평균의 참 표준편차 = s(ddof=0)/sqrt(n) = 1.589969
  정규근사 반너비 1.96*SE = 3.116281
  실제 붓스트랩 구간 (6.60, 12.80) 의 반너비 = 3.100000

예측한 너비의 비 = (1.96/t_4)*sqrt((n-1)/n) = 0.631399
실제 비 = 0.628100

n = 5 정규자료, 명목 95% (반복 20,000)
  t 구간의 포함률          0.9487
  백분위 붓스트랩의 포함률 0.8366

import numpy as np

np.random.seed(42)      # 아래 출력을 재현하려면 고정한다

def permutation_test(group_a, group_b, n_permutations=10000):
    """집단 표시를 무작위로 뒤섞어 "차이가 없다"는 세상을 흉내 낸다.

    H0가 참이면 어느 관측값이 어느 집단에 속하는지가 아무 상관이 없다.
    그러니 표시를 뒤섞은 자료들이 곧 귀무분포다.
    """
    combined = np.concatenate([group_a, group_b])
    observed_diff = np.mean(group_a) - np.mean(group_b)
    perm_differences = []
    for _ in range(n_permutations):
        # 비복원이다. 값 자체는 그대로 두고 순서만 바꾼다.
        # 붓스트랩이 복원추출인 것과 여기서 갈린다.
        np.random.shuffle(combined)
        perm_diff = np.mean(combined[:len(group_a)]) - np.mean(combined[len(group_a):])
        perm_differences.append(perm_diff)
    perm_differences = np.array(perm_differences)
    p_value = np.mean(np.abs(perm_differences) >= np.abs(observed_diff))
    return p_value, observed_diff, perm_differences

group_a = [8, 7, 9, 10, 6]
group_b = [5, 6, 4, 3, 7]
p_value, observed_diff, _ = permutation_test(group_a, group_b)
print(f"Observed Difference: {observed_diff:.2f}, P-value: {p_value:.4f}")