일원배치 분산분석: F-검정 절차¶

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
def load_data():
    """
    Load and preprocess plant growth data from the given URL.
    Returns:
        df (pd.DataFrame): The full DataFrame of plant growth data.
        data (tuple): A tuple containing weights for each group ('ctrl', 'trt1', 'trt2').
        df1 (int): Degrees of freedom between groups.
        df2 (int): Degrees of freedom within groups.
    """
    url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv'
    df = pd.read_csv(url, usecols=[1, 2])
    group_data = df.groupby('group')
    data_ctrl = group_data.get_group('ctrl').weight
    data_trt1 = group_data.get_group('trt1').weight
    data_trt2 = group_data.get_group('trt2').weight
    data = (data_ctrl, data_trt1, data_trt2)
    total_samples = data_ctrl.shape[0] + data_trt1.shape[0] + data_trt2.shape[0]
    num_groups = len(data)
    df1 = num_groups - 1
    df2 = total_samples - num_groups
    return df, data, df1, df2
def perform_anova(data_ctrl, data_trt1, data_trt2):
    """
    Perform one-way ANOVA on the given data.
    Returns:
        statistic (float): F-statistic of the ANOVA test.
        p_value (float): P-value of the ANOVA test.
    """
    statistic, p_value = stats.f_oneway(data_ctrl, data_trt1, data_trt2)
    print(f"\nANOVA Results:\nF-Statistic = {statistic:.4f}\nP-Value = {p_value:.4f}")
    return statistic, p_value
def plot_data(data_ctrl, data_trt1, data_trt2, df1, df2, statistic, p_value):
    """
    Plot boxplot of weights for each group and F-distribution with critical region.
    """
    fig, (ax_box, ax_pdf) = plt.subplots(1, 2, figsize=(12, 4))
    ax_box.boxplot([data_ctrl, data_trt1, data_trt2], labels=['ctrl', 'trt1', 'trt2'])
    ax_box.set_ylim(3, 7)
    ax_box.set_xlabel('Group')
    ax_box.set_ylabel('Weight')
    x_vals = np.linspace(0, 6, 100)
    pdf_vals = stats.f(df1, df2).pdf(x_vals)
    ax_pdf.plot(x_vals, pdf_vals, label='F-distribution PDF')
    ax_pdf.fill_between(x_vals[x_vals >= statistic], pdf_vals[x_vals >= statistic], color='red', alpha=0.3)
    ax_pdf.spines[['top','right']].set_visible(False)
    ax_pdf.spines[['bottom','left']].set_position("zero")
    ax_pdf.set_title("F-distribution and Critical Region")
    ax_pdf.legend()
    ax_pdf.annotate(f'P-Value = {p_value:.2%}', xy=(5.0, 0.1), xytext=(5.0, 0.8),
                    arrowprops=dict(color='k', width=0.2, headwidth=8), fontsize=12)
    plt.tight_layout()
    plt.show()
# 자료 읽기 → 분산분석 → 그림 순으로 돌린다
_, (data_ctrl, data_trt1, data_trt2), df1, df2 = load_data()
statistic, p_value = perform_anova(data_ctrl, data_trt1, data_trt2)
plot_data(data_ctrl, data_trt1, data_trt2, df1, df2, statistic, p_value)

ANOVA Results:
F-Statistic = 4.8461
P-Value = 0.0159

groups = [data_ctrl, data_trt1, data_trt2]
y = np.concatenate([g.values for g in groups])
n_i = np.array([len(g) for g in groups])
ybar_i = np.array([g.mean() for g in groups])
ybar = y.mean()

SST = ((y - ybar) ** 2).sum()
SSW = sum(((g.values - g.mean()) ** 2).sum() for g in groups)
SSB = (n_i * (ybar_i - ybar) ** 2).sum()
MSB, MSW = SSB / df1, SSW / df2
crit = stats.f(df1, df2).ppf(0.95)
print(f"SST = {SST:.5f} = SSW {SSW:.5f} + SSB {SSB:.5f}   (차이 {abs(SST - SSW - SSB):.1e})")
print(f"MSB = {MSB:.6f},  MSW = {MSW:.6f},  F = {MSB / MSW:.6f}")
print(f"임계값 F(0.05, 2, 27) = {crit:.4f}  ->  "
      f"{'H0 기각' if MSB / MSW > crit else 'H0 기각 못 함'}")

# E[MSW] = sigma^2, E[MSB] = sigma^2 + sum n_i (mu_i - mubar)^2 / (k-1) 을
# 모의실험으로 확인한다. 평균이 서로 다른 세 집단(각 10개)을 쓴다.
sigma, mu = 0.6, np.array([5.0, 4.7, 5.5])
theory_MSW = sigma ** 2
theory_MSB = sigma ** 2 + (10 * (mu - mu.mean()) ** 2).sum() / df1
rng = np.random.default_rng(0)
Y = rng.normal(mu[None, :, None], sigma, size=(200_000, 3, 10))
gm = Y.mean(axis=2)
msw = ((Y - gm[:, :, None]) ** 2).sum(axis=(1, 2)) / df2
msb = (10 * (gm - gm.mean(axis=1, keepdims=True)) ** 2).sum(axis=1) / df1
print(f"\nE[MSW]  이론 {theory_MSW:.5f}   모의 {msw.mean():.5f}")
print(f"E[MSB]  이론 {theory_MSB:.5f}   모의 {msb.mean():.5f}")

# 귀무가설(세 평균이 모두 같음) 아래에서는 E[MSB] = E[MSW] = sigma^2 이지만
# E[F] 는 1 이 아니라 d2/(d2-2) 다. 비의 기댓값은 기댓값의 비가 아니다.
Y0 = rng.normal(5.0, sigma, size=(200_000, 3, 10))
gm0 = Y0.mean(axis=2)
msw0 = ((Y0 - gm0[:, :, None]) ** 2).sum(axis=(1, 2)) / df2
msb0 = (10 * (gm0 - gm0.mean(axis=1, keepdims=True)) ** 2).sum(axis=1) / df1
print(f"\nH0 아래  E[MSB] 모의 {msb0.mean():.5f}   E[MSW] 모의 {msw0.mean():.5f}")
print(f"         E[MSB]/E[MSW] = {msb0.mean() / msw0.mean():.4f}")
print(f"         E[F] 모의 {np.mean(msb0 / msw0):.4f}   이론 d2/(d2-2) = {df2 / (df2 - 2):.4f}")

SST = 14.25843 = SSW 10.49209 + SSB 3.76634   (차이 4.0e-15)
MSB = 1.883170,  MSW = 0.388596,  F = 4.846088
임계값 F(0.05, 2, 27) = 3.3541  ->  H0 기각

E[MSW]  이론 0.36000   모의 0.35980
E[MSB]  이론 1.99333   모의 1.99497

H0 아래  E[MSB] 모의 0.36119   E[MSW] 모의 0.35958
         E[MSB]/E[MSW] = 1.0045
         E[F] 모의 1.0847   이론 d2/(d2-2) = 1.0800

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
def load_data():
    """
    Load and preprocess plant growth data for ANOVA.
    """
    url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/1dcc2bf5f955cc1224a3e1307256e1fe86b68dae/csv/datasets/PlantGrowth.csv'
    df = pd.read_csv(url, usecols=[1, 2])
    group_data = df.groupby('group')
    data_ctrl = group_data.get_group('ctrl').weight
    data_trt1 = group_data.get_group('trt1').weight
    data_trt2 = group_data.get_group('trt2').weight
    data = (data_ctrl, data_trt1, data_trt2)
    total_samples = data_ctrl.shape[0] + data_trt1.shape[0] + data_trt2.shape[0]
    num_groups = len(data)
    df1 = num_groups - 1
    df2 = total_samples - num_groups
    return df, data, df1, df2
def perform_anova(df):
    """
    Perform one-way ANOVA using statsmodels.
    """
    model = ols('weight ~ C(group)', data=df).fit()
    anova_results = anova_lm(model)
    statistic = anova_results['F'].iloc[0]
    p_value = anova_results['PR(>F)'].iloc[0]
    print("\nANOVA Results:\n", anova_results)
    return statistic, p_value
def plot_data(data_ctrl, data_trt1, data_trt2, df1, df2, statistic, p_value):
    """
    Plot boxplot of weights for each group and F-distribution with critical region.
    """
    fig, (ax_box, ax_pdf) = plt.subplots(1, 2, figsize=(12, 4))
    ax_box.boxplot([data_ctrl, data_trt1, data_trt2], labels=['ctrl', 'trt1', 'trt2'])
    ax_box.set_ylim(3, 7)
    ax_box.set_xlabel('Group')
    ax_box.set_ylabel('Weight')
    x_vals = np.linspace(0, 6, 100)
    pdf_vals = stats.f(df1, df2).pdf(x_vals)
    ax_pdf.plot(x_vals, pdf_vals, label='F-distribution PDF')
    ax_pdf.fill_between(x_vals[x_vals >= statistic], pdf_vals[x_vals >= statistic], color='red', alpha=0.3)
    ax_pdf.spines['top'].set_visible(False)
    ax_pdf.spines['right'].set_visible(False)
    ax_pdf.set_title("F-distribution and Critical Region")
    ax_pdf.legend()
    ax_pdf.annotate(f'P-Value = {p_value:.2%}', xy=(5.0, 0.1), xytext=(5.0, 0.8),
                    arrowprops=dict(color='k', width=0.2, headwidth=8), fontsize=12)
    plt.tight_layout()
    plt.show()
# 자료 읽기 → 분산분석 → 그림 순으로 돌린다
df, (data_ctrl, data_trt1, data_trt2), df1, df2 = load_data()
statistic, p_value = perform_anova(df)
plot_data(data_ctrl, data_trt1, data_trt2, df1, df2, statistic, p_value)

ANOVA Results:
             df    sum_sq   mean_sq         F   PR(>F)
C(group)   2.0   3.76634  1.883170  4.846088  0.01591
Residual  27.0  10.49209  0.388596       NaN      NaN

import numpy as np
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

# 분산분석표의 네 칸을 집단평균과 집단크기만으로 다시 만든다.
model = ols('weight ~ C(group)', data=df).fit()
table = anova_lm(model)

grp = df.groupby('group').weight
n_i = grp.count().values
ybar_i = grp.mean().values
ybar = df.weight.mean()
k, N = len(n_i), n_i.sum()

SSB = (n_i * (ybar_i - ybar) ** 2).sum()
SSW = ((n_i - 1) * grp.var(ddof=1).values).sum()
SST = ((df.weight - ybar) ** 2).sum()
MSB, MSW = SSB / (k - 1), SSW / (N - k)
F_hand = MSB / MSW
p_hand = stats.f.sf(F_hand, k - 1, N - k)

print(f"집단크기 {n_i},  집단평균 {ybar_i.round(5)},  전체평균 {ybar:.5f}")
print(f"\nSSB  = {SSB:.6f}   (표: {table['sum_sq'].iloc[0]:.6f})")
print(f"SSW  = {SSW:.6f}   (표: {table['sum_sq'].iloc[1]:.6f})")
print(f"SST  = {SST:.6f} = SSB + SSW = {SSB + SSW:.6f}   (차이 {abs(SST - SSB - SSW):.1e})")
print(f"MSB  = {MSB:.6f}   (표: {table['mean_sq'].iloc[0]:.6f})")
print(f"MSW  = {MSW:.6f}   (표: {table['mean_sq'].iloc[1]:.6f})")
print(f"F    = {F_hand:.6f}   (표: {table['F'].iloc[0]:.6f})")
print(f"p    = {p_hand:.6f}   (표: {table['PR(>F)'].iloc[0]:.6f})")

# 제1·2·3 유형이 정말 같은 수를 주는가
print()
for t in (1, 2, 3):
    row = anova_lm(model, typ=t).loc['C(group)']
    print(f"typ={t}:  sum_sq = {row['sum_sq']:.6f},  F = {row['F']:.6f}")

# f_oneway 와도 맞는가
res = stats.f_oneway(*[g.values for _, g in grp])
print(f"\nf_oneway:  F = {res.statistic:.6f},  p = {res.pvalue:.6f}")
print(f"anova_lm:  F = {statistic:.6f},  p = {p_value:.6f}")

집단크기 [10 10 10],  집단평균 [5.032 4.661 5.526],  전체평균 5.07300

SSB  = 3.766340   (표: 3.766340)
SSW  = 10.492090   (표: 10.492090)
SST  = 14.258430 = SSB + SSW = 14.258430   (차이 1.8e-15)
MSB  = 1.883170   (표: 1.883170)
MSW  = 0.388596   (표: 0.388596)
F    = 4.846088   (표: 4.846088)
p    = 0.015910   (표: 0.015910)

typ=1:  sum_sq = 3.766340,  F = 4.846088
typ=2:  sum_sq = 3.766340,  F = 4.846088
typ=3:  sum_sq = 3.766340,  F = 4.846088

f_oneway:  F = 4.846088,  p = 0.015910
anova_lm:  F = 4.846088,  p = 0.015910

import scipy.stats as stats
# 세 집단의 자료
water_group = [19, 18, 17, 18, 20]
energy_drink_group = [20, 22, 19, 21, 20]
coffee_group = [18, 17, 16, 19, 20]
# 일원배치 분산분석
f_statistic, p_value = stats.f_oneway(water_group, energy_drink_group, coffee_group)
# 결과 출력
print(f"F-statistic: {f_statistic:.2f}")
print(f"P-value: {p_value:.4f}")

F-statistic: 4.86
P-value: 0.0284

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
# 주어진 F 값과 자유도
df_between = 3 - 1
df_within = 15 - 3
f_statistic = 4.88
# p-값은 F 분포의 오른쪽 꼬리 넓이다
p_value = stats.f.sf(f_statistic, df_between, df_within)
print(f"{f_statistic = :.04f}")
print(f"{p_value = :.04f}")
# 그림과 축을 만든다
fig, ax = plt.subplots(figsize=(12, 4))
# 통계량 왼쪽 구간
x = np.linspace(0, f_statistic, 500)
y = stats.f.pdf(x, df_between, df_within)
ax.plot(x, y, color='blue', linewidth=3)
# 왼쪽을 칠한다 — 기각하지 않는 쪽
x = np.concatenate([[0], x, [f_statistic], [0]])
y = np.concatenate([[0], y, [0], [0]])
ax.fill(x, y, color='blue', alpha=0.1)
# 통계량 오른쪽 꼬리
x = np.linspace(f_statistic, 20, 500)
y = stats.f.pdf(x, df_between, df_within)
ax.plot(x, y, color='red', linewidth=3)
# 오른쪽을 칠한다 — 이 넓이가 p-값이다
x = np.concatenate([[f_statistic], x, [20], [f_statistic]])
y = np.concatenate([[0], y, [0], [0]])
ax.fill(x, y, color='red', alpha=0.1)
# p-값을 화살표로 가리킨다
xy = ((f_statistic + 15.0) / 2, 0.01)
xytext = (f_statistic + 3, 0.5)
arrowprops = dict(color='black', width=0.2, headwidth=8)
ax.annotate(f'{p_value = :.02%}', xy, xytext=xytext, fontsize=15, arrowprops=arrowprops)
# 축과 테두리를 다듬는다
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('F-value')
ax.set_ylabel('Probability Density')
ax.set_title('F-distribution with Highlighted p-value Region')
plt.show()

f_statistic = 4.8800
p_value = 0.0281

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# 웹페이지 네 개의 체류시간 자료
four_sessions = pd.DataFrame({
    'Time': [164, 178, 175, 155, 172, 182, 180, 179, 165, 166,
             172, 161, 171, 173, 158, 161, 179, 159, 167, 162],
    'Page': ['Page 1']*5 + ['Page 2']*5 + ['Page 3']*5 + ['Page 4']*5
})
# 관측된 집단평균들의 분산. 이 값이 클수록 집단 차이가 크다는 뜻이다.
obs_variance = four_sessions.groupby('Page')['Time'].mean().var()
print(f"Observed variance of means: {obs_variance:.2f}")
# 순열검정: 집단 이름표를 뒤섞어 귀무분포를 만든다
def perm_test_anova(df, group_col='Page', value_col='Time', n_perms=3000):
    """
    Permutation test for ANOVA using variance of group means.
    Parameters:
    -----------
    df : DataFrame
        Data with groups and values
    group_col : str
        Name of column with group labels
    value_col : str
        Name of column with values
    n_perms : int
        Number of permutations
    Returns:
    --------
    p_value : float
        Permutation test p-value
    perm_vars : array
        Permuted variances
    """
    groups = df[group_col].unique()
    group_sizes = {g: (df[group_col] == g).sum() for g in groups}
    obs_var = df.groupby(group_col)[value_col].mean().var()
    perm_vars = np.zeros(n_perms)
    for i in range(n_perms):
        # 값을 뒤섞어 같은 크기의 집단으로 다시 나눈다
        shuffled_values = np.random.permutation(df[value_col].values)
        perm_df = df.copy()
        perm_df[value_col] = shuffled_values
        # 뒤섞은 자료에서 집단평균들의 분산을 구한다
        perm_vars[i] = perm_df.groupby(group_col)[value_col].mean().var()
    p_value = np.mean(perm_vars >= obs_var)
    return p_value, perm_vars, obs_var
# 순열검정 실행
np.random.seed(42)
p_val, perm_vars, obs_var = perm_test_anova(four_sessions)
print(f"Permutation test p-value: {p_val:.4f}")
print(f"Conclusion: {'Reject H0' if p_val < 0.05 else 'Fail to reject H0'}")
# 귀무분포에 관측값을 얹어 본다
fig, ax = plt.subplots(figsize=(10, 6))
ax.hist(perm_vars, bins=30, alpha=0.7, color='steelblue', edgecolor='black')
ax.axvline(obs_var, color='red', linewidth=2, label=f'Observed = {obs_var:.2f}')
ax.set_xlabel('Variance of Group Means')
ax.set_ylabel('Frequency')
ax.set_title(f'Permutation Distribution of Group Mean Variance (p={p_val:.3f})')
ax.legend()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
plt.show()

Observed variance of means: 14.92
Permutation test p-value: 0.3673
Conclusion: Fail to reject H0

y = four_sessions['Time'].values
N, k, n = len(y), four_sessions['Page'].nunique(), 5
S2 = y.var(ddof=1)

theory_mean = S2 / n
mc_se = perm_vars.std(ddof=1) / np.sqrt(len(perm_vars))
print(f"묶은 20개 값의 표본분산 S^2 = {S2:.4f}")
print(f"순열분포 평균   이론 S^2/n = {theory_mean:.4f}")
print(f"                모의       = {perm_vars.mean():.4f}  (MC 표준오차 {mc_se:.4f})")
print(f"                z = {(perm_vars.mean() - theory_mean) / mc_se:+.2f}")
print(f"관측 분산 {obs_var:.4f} 는 순열평균보다 {obs_var - theory_mean:+.4f} 크다")
print(f"순열 p-값 = {p_val:.4f}")

묶은 20개 값의 표본분산 S^2 = 68.0500
순열분포 평균   이론 S^2/n = 13.6100
                모의       = 13.6739  (MC 표준오차 0.1802)
                z = +0.35
관측 분산 14.9167 는 순열평균보다 +1.3067 크다
순열 p-값 = 0.3673

from scipy import stats
def perm_test_anova_f(df, group_col='Page', value_col='Time', n_perms=3000):
    """
    Permutation test for ANOVA using F-statistic.
    """
    groups = df[group_col].unique()
    group_sizes = {g: (df[group_col] == g).sum() for g in groups}
    # 관측된 F 값
    model = smf.ols(f'{value_col} ~ C({group_col})', data=df).fit()
    anova_table = sm.stats.anova_lm(model)
    f_obs = anova_table['F'].iloc[0]
    perm_f_stats = np.zeros(n_perms)
    for i in range(n_perms):
        # 값을 뒤섞는다
        shuffled_values = np.random.permutation(df[value_col].values)
        perm_df = df.copy()
        perm_df[value_col] = shuffled_values
        # 뒤섞은 자료의 F 값
        perm_model = smf.ols(f'{value_col} ~ C({group_col})', data=perm_df).fit()
        perm_anova = sm.stats.anova_lm(perm_model)
        perm_f_stats[i] = perm_anova['F'].iloc[0]
    p_value = np.mean(perm_f_stats >= f_obs)
    return p_value, perm_f_stats, f_obs
# F 를 통계량으로 쓴 순열검정
import statsmodels.formula.api as smf
import statsmodels.api as sm
p_val_f, perm_f, obs_f = perm_test_anova_f(four_sessions)
print(f"\nF-statistic based permutation test:")
print(f"Observed F: {obs_f:.4f}")
print(f"p-value: {p_val_f:.4f}")

F-statistic based permutation test:
Observed F: 1.1161
p-value: 0.3550

# 같은 순열을 두 통계량에 함께 먹인다. 그래야 '통계량의 차이'만 남는다.
y = four_sessions['Time'].values
N, k, n = len(y), 4, 5
SST = ((y - y.mean()) ** 2).sum()
V_obs = four_sessions.groupby('Page')['Time'].mean().var()
SSB_obs = n * (k - 1) * V_obs
F_obs = (SSB_obs / (k - 1)) / ((SST - SSB_obs) / (N - k))
print(f"SST = {SST:.4f} (순열에 불변)")
print(f"관측  V = {V_obs:.6f}   SSB = {SSB_obs:.6f}   F = {F_obs:.6f}")

rng = np.random.default_rng(11)
B = 20000
V_perm = np.empty(B)
F_perm = np.empty(B)
for b in range(B):
    gm = rng.permutation(y).reshape(k, n).mean(axis=1)
    V_perm[b] = gm.var(ddof=1)
    ssb = n * ((gm - y.mean()) ** 2).sum()
    F_perm[b] = (ssb / (k - 1)) / ((SST - ssb) / (N - k))

print(f"\n같은 {B} 개 순열로 잰 p-값")
print(f"  V 기준: {np.mean(V_perm >= V_obs):.6f}")
print(f"  F 기준: {np.mean(F_perm >= F_obs):.6f}")

# 두 통계량의 순서가 뒤바뀐 쌍이 하나라도 있는가.
order = np.argsort(V_perm, kind='mergesort')
Vs, Fs = V_perm[order], F_perm[order]
gap = Vs[1:] - Vs[:-1]
down = Fs[:-1] > Fs[1:]
print(f"  V 가 뚜렷이 커지는데 F 는 작아진 쌍: {np.sum((gap > 1e-10) & down)} 개")
print(f"  동점(차이 1e-10 이하)인데 뒤바뀐 쌍: {np.sum((gap <= 1e-10) & down)} 개"
      f"  (최대 차이 {gap[(gap <= 1e-10) & down].max():.1e})")

# F 를 V 로 직접 예측해 본다.
pred = (N - k) / (k - 1) * (n * (k - 1) * V_perm) / (SST - n * (k - 1) * V_perm)
print(f"  F 와 공식값의 최대 차이: {np.abs(pred - F_perm).max():.1e}")

print(f"\nB=3000, p≈0.36 에서 MC 표준오차 = {np.sqrt(0.36 * 0.64 / 3000):.4f}")
print(f"독립인 두 모의실험 차이의 표준오차 = {np.sqrt(2 * 0.36 * 0.64 / 3000):.4f}")
print(f"실제 차이 |0.3673 - 0.3550| = {abs(0.3673 - 0.3550):.4f}")

SST = 1292.9500 (순열에 불변)
관측  V = 14.916667   SSB = 223.750000   F = 1.116099

같은 20000 개 순열로 잰 p-값
  V 기준: 0.368000
  F 기준: 0.368000
  V 가 뚜렷이 커지는데 F 는 작아진 쌍: 0 개
  동점(차이 1e-10 이하)인데 뒤바뀐 쌍: 289 개  (최대 차이 3.6e-15)
  F 와 공식값의 최대 차이: 1.1e-14

B=3000, p≈0.36 에서 MC 표준오차 = 0.0088
독립인 두 모의실험 차이의 표준오차 = 0.0124
실제 차이 |0.3673 - 0.3550| = 0.0123

# 모수적 분산분석. 순열검정 결과와 견준다.
f_stat, p_param = stats.f_oneway(
    four_sessions[four_sessions.Page == 'Page 1']['Time'],
    four_sessions[four_sessions.Page == 'Page 2']['Time'],
    four_sessions[four_sessions.Page == 'Page 3']['Time'],
    four_sessions[four_sessions.Page == 'Page 4']['Time']
)
print(f"\nParametric ANOVA:")
print(f"F-statistic: {f_stat:.4f}")
print(f"p-value: {p_param:.4f}")
print(f"\nPermutation ANOVA (variance-based):")
print(f"p-value: {p_val:.4f}")

Parametric ANOVA:
F-statistic: 1.1161
p-value: 0.3718

Permutation ANOVA (variance-based):
p-value: 0.3673

import math

p_perm, B = p_val, 3000
se = math.sqrt(p_perm * (1 - p_perm) / B)
print(f"순열 p-값 = {p_perm:.4f},  MC 표준오차 = {se:.4f}")
print(f"95% 구간 = [{p_perm - 1.96 * se:.4f}, {p_perm + 1.96 * se:.4f}]")
print(f"모수적 p-값 {p_param:.4f} 와의 차이 = {abs(p_param - p_perm):.4f}"
      f" = {abs(p_param - p_perm) / se:.2f} SE")

n_assign = math.comb(20, 5) * math.comb(15, 5) * math.comb(10, 5)
print(f"\n가능한 배정의 수 = {n_assign:,}")

# 가정이 성립하는지도 살펴 둔다.
resid = four_sessions.Time - four_sessions.groupby('Page').Time.transform('mean')
print(f"\n잔차의 Shapiro-Wilk: W = {stats.shapiro(resid).statistic:.4f}, "
      f"p = {stats.shapiro(resid).pvalue:.4f}")
by_page = [four_sessions[four_sessions.Page == g].Time.values
           for g in four_sessions.Page.unique()]
lev = stats.levene(*by_page)
print(f"Levene(중앙값 기준): W = {lev.statistic:.4f}, p = {lev.pvalue:.4f}")

순열 p-값 = 0.3673,  MC 표준오차 = 0.0088
95% 구간 = [0.3501, 0.3846]
모수적 p-값 0.3718 와의 차이 = 0.0045 = 0.51 SE

가능한 배정의 수 = 11,732,745,024

잔차의 Shapiro-Wilk: W = 0.9468, p = 0.3211
Levene(중앙값 기준): W = 0.0736, p = 0.9733

import numpy as np
from scipy import stats

rng = np.random.default_rng(9090)

def welch_anova(*groups):
    k = len(groups)
    n = np.array([len(g) for g in groups], float)
    m = np.array([g.mean() for g in groups])
    v = np.array([g.var(ddof=1) for g in groups])
    w = n / v
    W = w.sum()
    m_tilde = (w * m).sum() / W
    tmp = np.sum((1 - w / W)**2 / (n - 1))
    F = ((w * (m - m_tilde)**2).sum() / (k - 1)) \
        / (1 + 2 * (k - 2) / (k * k - 1) * tmp)
    return F, stats.f.sf(F, k - 1, (k * k - 1) / (3 * tmp))

def perm_f(groups, B, rng):
    """F 통계량을 쓰는 순열검정."""
    obs = stats.f_oneway(*groups).statistic
    sizes = [len(g) for g in groups]
    pooled = np.concatenate(groups)
    cnt = 1
    for _ in range(B):
        rng.shuffle(pooled)
        parts = np.split(pooled, np.cumsum(sizes)[:-1])
        cnt += stats.f_oneway(*parts).statistic >= obs - 1e-12
    return cnt / (B + 1)

M = 2_000
print(f"{'상황':>20s} {'고전 F':>8s} {'순열 F':>8s} {'Welch':>8s} {'K-W':>8s}")
cases = [("정규 등분산 H0", (0, 0, 0), (1, 1, 1), (10, 10, 10)),
         ("정규 이분산 H0", (0, 0, 0), (1, 1, 4), (10, 10, 10)),
         ("이분산·불균형 H0", (0, 0, 0), (1, 1, 4), (20, 20, 5)),
         ("지수분포 H0", None, None, (10, 10, 10)),
         ("로그정규 H0", "ln", None, (15, 15, 15)),
         ("정규 등분산 H1", (0, 1, 2), (1, 1, 1), (10, 10, 10))]
for label, mu, sig, ns in cases:
    a = b = c = d = 0
    for _ in range(M):
        if mu is None:
            g = [rng.exponential(1, n) for n in ns]
        elif mu == "ln":
            g = [rng.lognormal(0, 1, n) for n in ns]
        else:
            g = [rng.normal(m, s, n) for m, s, n in zip(mu, sig, ns)]
        a += stats.f_oneway(*g).pvalue < 0.05
        b += perm_f(g, 299, rng) < 0.05
        c += welch_anova(*g)[1] < 0.05
        d += stats.kruskal(*g).pvalue < 0.05
    print(f"{label:>20s} {a / M:8.4f} {b / M:8.4f} {c / M:8.4f} {d / M:8.4f}")

                  상황     고전 F     순열 F    Welch      K-W
           정규 등분산 H0   0.0530   0.0480   0.0530   0.0410
           정규 이분산 H0   0.0895   0.0915   0.0455   0.0720
          이분산·불균형 H0   0.3405   0.3280   0.0625   0.1210
             지수분포 H0   0.0520   0.0520   0.0530   0.0475
             로그정규 H0   0.0315   0.0465   0.0375   0.0425
           정규 등분산 H1   0.9740   0.9675   0.9630   0.9590

import numpy as np
from scipy import stats

groups = {
    "A": np.array([4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14]),
    "B": np.array([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]),
    "C": np.array([6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26]),
}
names = list(groups)
k, n = len(names), 10
N = k * n
m = np.array([groups[x].mean() for x in names])
grand = np.concatenate(list(groups.values())).mean()
SSE = sum(((g - g.mean())**2).sum() for g in groups.values())
MSE = SSE / (N - k)
SST = n * ((m - grand)**2).sum()
print(f"SST = {SST:.4f},  MSE = {MSE:.4f},  "
      f"옴니버스 F = {SST / (k - 1) / MSE:.4f}\n")

contrasts = {"A 대 (B,C) 평균": np.array([2., -1, -1]),
             "B 대 C        ": np.array([0., 1, -1])}
total = 0.0
for label, c in contrasts.items():
    L = (c * m).sum()
    ss = L**2 / ((c**2).sum() / n)
    total += ss
    print(f"  {label}: L = {L:+.4f},  SS = {ss:.4f},  "
          f"F = {ss / MSE:.4f},  p = {stats.f.sf(ss / MSE, 1, N - k):.4f}")
print(f"\n  대비 SS 의 합 = {total:.4f}   (SST = {SST:.4f})")
c1, c2 = list(contrasts.values())
print(f"  직교성 확인: Σ c1·c2 / n = {(c1 * c2).sum() / n:.1f}")

SST = 3.7663,  MSE = 0.3886,  옴니버스 F = 4.8461

  A 대 (B,C) 평균: L = -0.1230,  SS = 0.0252,  F = 0.0649,  p = 0.8009
  B 대 C        : L = -0.8650,  SS = 3.7411,  F = 9.6273,  p = 0.0045

  대비 SS 의 합 = 3.7663   (SST = 3.7663)
  직교성 확인: Σ c1·c2 / n = 0.0

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

rng = np.random.default_rng(2323)
M, n, sigma, k = 5_000, 15, 1.0, 4
c_linear = np.array([-3., -1, 1, 3])      # 선형 추세
c_last = np.array([-1., -1, -1, 3])       # 마지막 대 나머지

for label, mu in [("선형 추세 (0,0.4,0.8,1.2)", np.arange(4) * 0.4),
                  ("하나만 다름 (0,0,0,1.2)", np.array([0., 0, 0, 1.2]))]:
    a = b = c = 0
    for _ in range(M):
        g = [rng.normal(mm, sigma, n) for mm in mu]
        m_hat = np.array([x.mean() for x in g])
        MSE = sum(((x - x.mean())**2).sum() for x in g) / (k * n - k)
        a += stats.f_oneway(*g).pvalue < 0.05
        for cc, which in [(c_linear, 0), (c_last, 1)]:
            ss = (cc * m_hat).sum()**2 / ((cc**2).sum() / n)
            hit = stats.f.sf(ss / MSE, 1, k * n - k) < 0.05
            if which == 0:
                b += hit
            else:
                c += hit
    print(f"{label:>26s}:  옴니버스 {a / M:.4f}   "
          f"선형대비 {b / M:.4f}   마지막대비 {c / M:.4f}")

     선형 추세 (0,0.4,0.8,1.2):  옴니버스 0.8108   선형대비 0.9274   마지막대비 0.7538
        하나만 다름 (0,0,0,1.2):  옴니버스 0.9198   선형대비 0.8632   마지막대비 0.9756

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

import numpy as np
from scipy import stats

def n_per_group(f, k, power=0.80, alpha=0.05):
    for n in range(3, 100_000):
        N = k * n
        lam = f**2 * N
        crit = stats.f.ppf(1 - alpha, k - 1, N - k)
        if stats.ncf.sf(crit, k - 1, N - k, lam) >= power:
            return n

print(f"{'k':>4s} " + " ".join(f"{'f=' + str(f):>9s}"
                               for f in [0.10, 0.25, 0.40]))
for k in [2, 3, 4, 5, 6, 8, 10]:
    row = " ".join(f"{n_per_group(f, k):9d}" for f in [0.10, 0.25, 0.40])
    print(f"{k:4d} {row}")

print(f"\n{'k':>4s} {'총 N (f=0.25)':>14s}")
for k in [2, 3, 4, 5, 6, 8, 10]:
    print(f"{k:4d} {k * n_per_group(0.25, k):14d}")

   k     f=0.1    f=0.25     f=0.4
   2       394        64        26
   3       323        53        22
   4       274        45        19
   5       240        40        16
   6       215        36        15
   8       181        30        13
  10       158        26        11

   k   총 N (f=0.25)
   2            128
   3            159
   4            180
   5            200
   6            216
   8            240
  10            260

print(f"{'k':>4s} {'μ 배치':>28s} {'f':>8s}")
for k, mu in [(3, [0, 0.5, 1.0]), (3, [0, 0, 1.0]),
              (5, [0, 0.25, 0.5, 0.75, 1.0]), (5, [0, 0, 0, 0, 1.0])]:
    mu = np.array(mu, float)
    f = np.sqrt(((mu - mu.mean())**2).sum() / k)
    print(f"{k:4d} {str(mu.tolist()):>28s} {f:8.4f}")

   k                         μ 배치        f
   3              [0.0, 0.5, 1.0]   0.4082
   3              [0.0, 0.0, 1.0]   0.4714
   5  [0.0, 0.25, 0.5, 0.75, 1.0]   0.3536
   5    [0.0, 0.0, 0.0, 0.0, 1.0]   0.4000

세 집단 이상의 평균을 비교한다
    │
    ├─ ① 자료 구조 확인
    │     ├─ 반복측정인가 → 반복측정 분산분석 / 혼합모형
    │     ├─ 군집 구조가 있는가 → 혼합모형
    │     └─ 독립 표본 → 계속
    │
    ├─ ② 그림을 그린다
    │     상자그림, 집단별 정규분위수그림, 잔차 그림
    │
    ├─ ③ 검정을 고른다 (사전에)
    │     ├─ 등분산이 합당 ────→ 고전 F
    │     ├─ 이분산이 의심 ────→ Welch 분산분석  ← 기본 권장
    │     ├─ 비정규 + 등분산 ──→ 순열 F
    │     └─ 순위가 관심 ──────→ 크러스컬·월리스
    │
    ├─ ④ 사전에 정한 대비가 있으면 그것을 먼저 본다
    │
    ├─ ⑤ 옴니버스 검정 + 효과크기 (η², ω², f)
    │
    └─ ⑥ 유의하면 사후비교
          ├─ 모든 쌍 ──────→ 투키 (등분산) / 게임스·하월 (이분산)
          ├─ 대조군 대비 ──→ 더넷
          └─ 자료를 보고 만든 대비 → 셰페