import numpy as np
from scipy import stats
from scipy.special import gamma

def cohens_d(group1, group2):
    """Cohen's d を計算する"""
    n1, n2 = len(group1), len(group2)
    s1, s2 = np.std(group1, ddof=1), np.std(group2, ddof=1)

    # プールされた標準偏差
    sp = np.sqrt(((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / (n1 + n2 - 2))

    d = (np.mean(group1) - np.mean(group2)) / sp
    return d

def glass_delta(group1, group2):
    """Glass' Δ を計算する（group2 を統制群とする）"""
    return (np.mean(group1) - np.mean(group2)) / np.std(group2, ddof=1)

def hedges_g(group1, group2):
    """Hedges' g を計算する（小標本バイアス補正）"""
    d = cohens_d(group1, group2)
    n1, n2 = len(group1), len(group2)
    df = n1 + n2 - 2

    # 補正係数の近似式
    j = 1 - 3 / (4 * df - 1)

    return d * j

def eta_squared(groups):
    """η²（イータ二乗）を計算する"""
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)

    ss_between = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)
    ss_total = np.sum((all_data - grand_mean)**2)

    return ss_between / ss_total

def omega_squared(groups):
    """ω²（オメガ二乗）を計算する"""
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)
    k = len(groups)
    N = len(all_data)

    ss_between = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)
    ss_within = sum(np.sum((g - np.mean(g))**2) for g in groups)
    ss_total = np.sum((all_data - grand_mean)**2)

    ms_within = ss_within / (N - k)

    return (ss_between - (k - 1) * ms_within) / (ss_total + ms_within)

def odds_ratio(a, b, c, d):
    """オッズ比を計算する"""
    OR = (a * d) / (b * c)
    log_OR = np.log(OR)
    se_log_OR = np.sqrt(1/a + 1/b + 1/c + 1/d)

    # 95%信頼区間
    ci_lower = np.exp(log_OR - 1.96 * se_log_OR)
    ci_upper = np.exp(log_OR + 1.96 * se_log_OR)

    return OR, ci_lower, ci_upper


# --- 具体例 ---
np.random.seed(42)

# 2群のデータ生成
group_a = np.random.normal(78.5, 10.2, 30)
group_b = np.random.normal(72.3, 11.5, 35)

print("=== 差の効果量 ===")
print(f"Cohen's d  = {cohens_d(group_a, group_b):.4f}")
print(f"Glass' Δ   = {glass_delta(group_a, group_b):.4f}")
print(f"Hedges' g  = {hedges_g(group_a, group_b):.4f}")

# t検定
t_stat, p_value = stats.ttest_ind(group_a, group_b)
print(f"\nt検定: t = {t_stat:.4f}, p = {p_value:.4f}")

# t統計量からCohen's dを逆算
d_from_t = t_stat * np.sqrt(1/len(group_a) + 1/len(group_b))
print(f"t統計量から算出したd = {d_from_t:.4f}")

# 分散分析の効果量
print("\n=== 関連の効果量（ANOVA） ===")
group_c = np.random.normal(75.0, 10.0, 28)
groups = [group_a, group_b, group_c]

print(f"η² = {eta_squared(groups):.4f}")
print(f"ω² = {omega_squared(groups):.4f}")

# ANOVA
f_stat, p_anova = stats.f_oneway(*groups)
print(f"\nANOVA: F = {f_stat:.4f}, p = {p_anova:.4f}")

# オッズ比
print("\n=== オッズ比 ===")
a, b, c, d_val = 40, 60, 20, 80  # 2×2分割表
OR, ci_low, ci_high = odds_ratio(a, b, c, d_val)
print(f"2×2分割表: a={a}, b={b}, c={c}, d={d_val}")
print(f"オッズ比 = {OR:.4f}")
print(f"95%信頼区間 = [{ci_low:.4f}, {ci_high:.4f}]")

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

np.random.seed(42)

# 異なる効果量でのCohens' dの可視化
effect_sizes = [0.2, 0.5, 0.8, 1.2]
labels = ['Small (d=0.2)', 'Medium (d=0.5)', 'Large (d=0.8)', 'Very Large (d=1.2)']
colors = ['#2196F3', '#FF9800', '#F44336', '#9C27B0']

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

x = np.linspace(-4, 6, 500)

for idx, (d, label, color) in enumerate(zip(effect_sizes, labels, colors)):
    # 群1: N(0, 1)、群2: N(d, 1)
    pdf1 = stats.norm.pdf(x, 0, 1)
    pdf2 = stats.norm.pdf(x, d, 1)

    axes[idx].fill_between(x, pdf1, alpha=0.3, color='blue', label='Group 1 ($\\mu=0$)')
    axes[idx].fill_between(x, pdf2, alpha=0.3, color=color, label=f'Group 2 ($\\mu={d}$)')
    axes[idx].plot(x, pdf1, 'b-', linewidth=2)
    axes[idx].plot(x, pdf2, '-', color=color, linewidth=2)

    # 重なりの割合（近似計算）
    overlap = np.trapz(np.minimum(pdf1, pdf2), x) * 100

    axes[idx].set_title(f"{label}\nOverlap: {overlap:.1f}%", fontsize=12)
    axes[idx].set_xlabel('Value', fontsize=11)
    axes[idx].set_ylabel('Density', fontsize=11)
    axes[idx].legend(fontsize=9)
    axes[idx].grid(True, alpha=0.3)

plt.suptitle("Cohen's d: Distribution Overlap for Different Effect Sizes",
             fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

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

np.random.seed(42)

# 固定された効果量で標本サイズを変えたときのp値
d_true = 0.3  # 小さな効果量
sample_sizes = np.arange(10, 510, 10)
n_simulations = 2000

mean_p_values = []
power_values = []
alpha = 0.05

for n in sample_sizes:
    p_vals = []
    rejects = 0

    for _ in range(n_simulations):
        group1 = np.random.normal(0, 1, n)
        group2 = np.random.normal(d_true, 1, n)
        _, p = stats.ttest_ind(group1, group2)
        p_vals.append(p)
        if p < alpha:
            rejects += 1

    mean_p_values.append(np.median(p_vals))
    power_values.append(rejects / n_simulations)

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

# p値の変化
ax1.plot(sample_sizes, mean_p_values, 'b-o', markersize=3, linewidth=1.5)
ax1.axhline(y=0.05, color='red', linestyle='--', linewidth=1.5, label='$\\alpha = 0.05$')
ax1.set_xlabel('Sample Size per Group', fontsize=12)
ax1.set_ylabel('Median p-value', fontsize=12)
ax1.set_title(f"Median p-value vs Sample Size (True $d$ = {d_true})", fontsize=13)
ax1.set_yscale('log')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 検出力の変化
ax2.plot(sample_sizes, power_values, 'r-s', markersize=3, linewidth=1.5)
ax2.axhline(y=0.8, color='gray', linestyle='--', linewidth=1.5, label='Power = 0.80')
ax2.set_xlabel('Sample Size per Group', fontsize=12)
ax2.set_ylabel('Power', fontsize=12)
ax2.set_title(f"Statistical Power vs Sample Size (True $d$ = {d_true})", fontsize=13)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(0, 1.05)

plt.tight_layout()
plt.show()

print(f"Cohen's d = {d_true} を80%の検出力で検出するには:")
# 理論値: n = (z_α/2 + z_β)² × 2 / d²
z_alpha = stats.norm.ppf(1 - 0.05/2)
z_beta = stats.norm.ppf(0.8)
n_required = 2 * ((z_alpha + z_beta) / d_true)**2
print(f"  理論上必要な標本サイズ（各群）: n ≈ {n_required:.0f}")

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

np.random.seed(42)

n_simulations = 5000
group_sizes = [10, 20, 50, 100]

# 真の効果量（母集団のη²）
# 3群: μ1=0, μ2=0, μ3=0（帰無仮説が真、η²_pop = 0）
sigma = 1.0

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.flatten()

for idx, n_per_group in enumerate(group_sizes):
    eta2_vals = []
    omega2_vals = []

    for _ in range(n_simulations):
        # H0のもとでデータ生成（3群、すべて同じ分布）
        g1 = np.random.normal(0, sigma, n_per_group)
        g2 = np.random.normal(0, sigma, n_per_group)
        g3 = np.random.normal(0, sigma, n_per_group)

        groups = [g1, g2, g3]
        all_data = np.concatenate(groups)
        grand_mean = np.mean(all_data)
        k = 3
        N = 3 * n_per_group

        ss_between = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)
        ss_within = sum(np.sum((g - np.mean(g))**2) for g in groups)
        ss_total = np.sum((all_data - grand_mean)**2)
        ms_within = ss_within / (N - k)

        eta2 = ss_between / ss_total
        omega2 = (ss_between - (k - 1) * ms_within) / (ss_total + ms_within)

        eta2_vals.append(eta2)
        omega2_vals.append(omega2)

    eta2_vals = np.array(eta2_vals)
    omega2_vals = np.array(omega2_vals)

    axes[idx].hist(eta2_vals, bins=50, density=True, alpha=0.6, color='blue',
                   label=f'$\\eta^2$ (mean={np.mean(eta2_vals):.4f})')
    axes[idx].hist(omega2_vals, bins=50, density=True, alpha=0.6, color='red',
                   label=f'$\\omega^2$ (mean={np.mean(omega2_vals):.4f})')
    axes[idx].axvline(x=0, color='gray', linestyle='--', linewidth=1.5,
                      label='True effect = 0')

    axes[idx].set_title(f'n per group = {n_per_group}', fontsize=13)
    axes[idx].set_xlabel('Effect Size', fontsize=11)
    axes[idx].set_ylabel('Density', fontsize=11)
    axes[idx].legend(fontsize=9)
    axes[idx].grid(True, alpha=0.3)

plt.suptitle('$\\eta^2$ vs $\\omega^2$ under $H_0$ (True Effect = 0)',
             fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

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

# さまざまなオッズ比の信頼区間を可視化
studies = [
    {'name': 'Study A', 'a': 40, 'b': 60, 'c': 20, 'd': 80},
    {'name': 'Study B', 'a': 15, 'b': 85, 'c': 10, 'd': 90},
    {'name': 'Study C', 'a': 55, 'b': 45, 'c': 25, 'd': 75},
    {'name': 'Study D', 'a': 30, 'b': 70, 'c': 28, 'd': 72},
    {'name': 'Study E', 'a': 8, 'b': 92, 'c': 3, 'd': 97},
]

fig, ax = plt.subplots(figsize=(10, 6))

for i, study in enumerate(studies):
    a, b, c, d = study['a'], study['b'], study['c'], study['d']
    OR = (a * d) / (b * c)
    log_OR = np.log(OR)
    se = np.sqrt(1/a + 1/b + 1/c + 1/d)
    ci_lower = np.exp(log_OR - 1.96 * se)
    ci_upper = np.exp(log_OR + 1.96 * se)

    color = 'red' if ci_lower > 1 or ci_upper < 1 else 'gray'
    ax.errorbar(OR, i, xerr=[[OR - ci_lower], [ci_upper - OR]],
                fmt='o', color=color, markersize=8, capsize=5,
                linewidth=2, elinewidth=2)
    ax.text(ci_upper + 0.3, i, f'OR = {OR:.2f} [{ci_lower:.2f}, {ci_upper:.2f}]',
            va='center', fontsize=10)

ax.axvline(x=1, color='gray', linestyle='--', linewidth=1.5, label='OR = 1 (No effect)')
ax.set_yticks(range(len(studies)))
ax.set_yticklabels([s['name'] for s in studies], fontsize=11)
ax.set_xlabel('Odds Ratio (95% CI)', fontsize=12)
ax.set_title('Forest Plot: Odds Ratios with 95% Confidence Intervals', fontsize=14)
ax.set_xscale('log')
ax.set_xlim(0.3, 15)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()