import numpy as np
from scipy import stats

# 2群のデータ
np.random.seed(42)
group1 = np.random.normal(50, 5, 20)   # σ = 5
group2 = np.random.normal(50, 8, 25)   # σ = 8

# F検定
s1_sq = np.var(group1, ddof=1)
s2_sq = np.var(group2, ddof=1)

# 大きい分散を分子にする
if s1_sq >= s2_sq:
    f_stat = s1_sq / s2_sq
    df1, df2 = len(group1) - 1, len(group2) - 1
else:
    f_stat = s2_sq / s1_sq
    df1, df2 = len(group2) - 1, len(group1) - 1

# 両側検定のp値
p_value = 2 * min(stats.f.cdf(f_stat, df1, df2),
                  1 - stats.f.cdf(f_stat, df1, df2))

print("=== 等分散のF検定 ===")
print(f"群1: n={len(group1)}, s² = {s1_sq:.4f}")
print(f"群2: n={len(group2)}, s² = {s2_sq:.4f}")
print(f"F = {f_stat:.4f}")
print(f"自由度: ({df1}, {df2})")
print(f"p値 (両側): {p_value:.4f}")
print(f"α = 0.05 での判定: {'棄却（分散に差あり）' if p_value < 0.05 else '棄却しない（等分散）'}")

import numpy as np
from scipy import stats

# 3種類の肥料のデータ
fertilizer_a = np.array([20, 22, 19, 21])
fertilizer_b = np.array([23, 25, 22, 24])
fertilizer_c = np.array([18, 20, 17, 19])

# scipy での一元配置ANOVA
f_stat, p_value = stats.f_oneway(fertilizer_a, fertilizer_b, fertilizer_c)

print("=== 一元配置分散分析 (scipy) ===")
print(f"群A: 平均 = {np.mean(fertilizer_a):.2f}")
print(f"群B: 平均 = {np.mean(fertilizer_b):.2f}")
print(f"群C: 平均 = {np.mean(fertilizer_c):.2f}")
print(f"F = {f_stat:.4f}")
print(f"p値: {p_value:.6f}")
print(f"α = 0.05 での判定: {'棄却（群間に差あり）' if p_value < 0.05 else '棄却しない'}")

import numpy as np
from scipy import stats

def one_way_anova(*groups):
    """
    一元配置分散分析をスクラッチで実装

    Parameters
    ----------
    *groups : array-like
        各群のデータ

    Returns
    -------
    anova_table : dict
        ANOVA表の各成分
    """
    k = len(groups)
    all_data = np.concatenate(groups)
    n_total = len(all_data)
    grand_mean = np.mean(all_data)

    # 群間平方和 (SSB)
    ssb = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)

    # 群内平方和 (SSW)
    ssw = sum(np.sum((g - np.mean(g))**2) for g in groups)

    # 全体平方和 (SST)
    sst = np.sum((all_data - grand_mean)**2)

    # 自由度
    df_between = k - 1
    df_within = n_total - k
    df_total = n_total - 1

    # 平均平方
    msb = ssb / df_between
    msw = ssw / df_within

    # F統計量
    f_stat = msb / msw

    # p値
    p_value = 1 - stats.f.cdf(f_stat, df_between, df_within)

    return {
        'SSB': ssb, 'SSW': ssw, 'SST': sst,
        'df_between': df_between, 'df_within': df_within, 'df_total': df_total,
        'MSB': msb, 'MSW': msw,
        'F': f_stat, 'p_value': p_value
    }


# 肥料データで実行
result = one_way_anova(
    np.array([20, 22, 19, 21]),
    np.array([23, 25, 22, 24]),
    np.array([18, 20, 17, 19])
)

print("=== ANOVA表 ===")
print(f"{'要因':<8} {'平方和':>10} {'自由度':>6} {'平均平方':>10} {'F値':>8} {'p値':>10}")
print("-" * 60)
print(f"{'群間':<8} {result['SSB']:>10.3f} {result['df_between']:>6} {result['MSB']:>10.3f} {result['F']:>8.3f} {result['p_value']:>10.6f}")
print(f"{'群内':<8} {result['SSW']:>10.3f} {result['df_within']:>6} {result['MSW']:>10.3f}")
print(f"{'全体':<8} {result['SST']:>10.3f} {result['df_total']:>6}")

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

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: F分布の形状
x = np.linspace(0, 6, 500)
for (d1, d2), color in zip([(2, 9), (3, 20), (5, 30), (10, 50)],
                            ['red', 'blue', 'green', 'purple']):
    pdf = stats.f.pdf(x, d1, d2)
    axes[0].plot(x, pdf, color=color, linewidth=2,
                 label=f'$F_{{{d1},{d2}}}$')

axes[0].set_xlabel('$x$', fontsize=12)
axes[0].set_ylabel('Probability Density', fontsize=12)
axes[0].set_title('F-distribution', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, 1.0)

# 右: ANOVAの検定結果
f_obs = 15.2
df1, df2 = 2, 9
x = np.linspace(0, 20, 500)
pdf = stats.f.pdf(x, df1, df2)

axes[1].plot(x, pdf, 'k-', linewidth=2, label=f'$F_{{{df1},{df2}}}$')

# 棄却域
f_crit = stats.f.ppf(0.95, df1, df2)
x_reject = x[x >= f_crit]
axes[1].fill_between(x_reject, stats.f.pdf(x_reject, df1, df2),
                     color='red', alpha=0.3, label=f'Rejection ($\\alpha$ = 0.05)')
axes[1].axvline(f_crit, color='red', linestyle='--', linewidth=1.5,
                label=f'$F_{{crit}}$ = {f_crit:.3f}')

# 観測されたF値
p_value = 1 - stats.f.cdf(f_obs, df1, df2)
axes[1].axvline(f_obs, color='blue', linewidth=2.5,
                label=f'$F_{{obs}}$ = {f_obs:.1f} (p = {p_value:.5f})')

axes[1].set_xlabel('$F$', fontsize=12)
axes[1].set_ylabel('Probability Density', fontsize=12)
axes[1].set_title('One-way ANOVA: Fertilizer Example', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(bottom=0)

plt.tight_layout()
plt.show()

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

# より多いデータで ANOVA を実施
np.random.seed(42)
group_a = np.random.normal(50, 5, 30)
group_b = np.random.normal(55, 5, 30)
group_c = np.random.normal(52, 5, 30)
group_d = np.random.normal(48, 5, 30)

f_stat, p_value = stats.f_oneway(group_a, group_b, group_c, group_d)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: 箱ひげ図
data = [group_a, group_b, group_c, group_d]
labels = ['Group A', 'Group B', 'Group C', 'Group D']
bp = axes[0].boxplot(data, labels=labels, patch_artist=True,
                     boxprops=dict(facecolor='lightblue', edgecolor='black'))
axes[0].set_ylabel('Value', fontsize=12)
axes[0].set_title(f'Box Plot (F = {f_stat:.2f}, p = {p_value:.4f})', fontsize=13)
axes[0].grid(True, alpha=0.3, axis='y')

# 各群の平均をプロット
means = [np.mean(g) for g in data]
axes[0].scatter(range(1, 5), means, color='red', zorder=5, s=80,
                label='Mean', marker='D')
axes[0].legend(fontsize=11)

# 右: 各群の分布
x_range = np.linspace(30, 70, 200)
colors = ['blue', 'red', 'green', 'orange']
for g, label, color in zip(data, labels, colors):
    mu, sigma = np.mean(g), np.std(g, ddof=1)
    axes[1].plot(x_range, stats.norm.pdf(x_range, mu, sigma),
                 color=color, linewidth=2, label=f'{label} ($\\bar{{x}}$={mu:.1f})')

axes[1].set_xlabel('Value', fontsize=12)
axes[1].set_ylabel('Density', fontsize=12)
axes[1].set_title('Estimated Distributions', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"F = {f_stat:.4f}, p = {p_value:.6f}")
print(f"α = 0.05 での判定: {'棄却（群間に差あり）' if p_value < 0.05 else '棄却しない'}")