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

# 正規分布での片側・両側検定のp値の可視化
z_obs = 2.0  # 観測された検定統計量

z = np.linspace(-4, 4, 1000)
pdf = stats.norm.pdf(z)

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

# 片側検定（右側）
ax = axes[0]
ax.plot(z, pdf, 'k-', linewidth=2)
z_fill = z[z >= z_obs]
ax.fill_between(z_fill, stats.norm.pdf(z_fill), alpha=0.4, color='red',
                label=f'p = {1 - stats.norm.cdf(z_obs):.4f}')
ax.axvline(x=z_obs, color='red', linestyle='--', linewidth=1.5)
ax.set_title(f'One-sided test ($z_{{obs}}$ = {z_obs})', fontsize=12)
ax.set_xlabel('z')
ax.set_ylabel('Density')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

# 両側検定
ax = axes[1]
ax.plot(z, pdf, 'k-', linewidth=2)
z_fill_right = z[z >= z_obs]
z_fill_left = z[z <= -z_obs]
ax.fill_between(z_fill_right, stats.norm.pdf(z_fill_right), alpha=0.4, color='blue')
ax.fill_between(z_fill_left, stats.norm.pdf(z_fill_left), alpha=0.4, color='blue',
                label=f'p = {2 * (1 - stats.norm.cdf(z_obs)):.4f}')
ax.axvline(x=z_obs, color='blue', linestyle='--', linewidth=1.5)
ax.axvline(x=-z_obs, color='blue', linestyle='--', linewidth=1.5)
ax.set_title(f'Two-sided test ($|z_{{obs}}|$ = {z_obs})', fontsize=12)
ax.set_xlabel('z')
ax.set_ylabel('Density')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("pvalue_one_two_sided.png", dpi=150, bbox_inches="tight")
plt.show()

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

np.random.seed(42)
n = 30
sigma = 1
mu0 = 0
n_simulations = 10000

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

# H0が真の場合のp値の分布
p_values_h0 = []
for _ in range(n_simulations):
    sample = np.random.normal(mu0, sigma, n)
    z = sample.mean() / (sigma / np.sqrt(n))
    p = 2 * (1 - stats.norm.cdf(abs(z)))  # 両側検定
    p_values_h0.append(p)

axes[0].hist(p_values_h0, bins=50, density=True, alpha=0.7, color='skyblue', edgecolor='white')
axes[0].axhline(y=1, color='red', linewidth=2, linestyle='--', label='Uniform(0,1)')
axes[0].set_xlabel('p-value', fontsize=12)
axes[0].set_ylabel('Density', fontsize=12)
axes[0].set_title('p-value Distribution under $H_0$ (should be Uniform)', fontsize=12)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# H1が真の場合のp値の分布
mu1 = 0.5
p_values_h1 = []
for _ in range(n_simulations):
    sample = np.random.normal(mu1, sigma, n)
    z = sample.mean() / (sigma / np.sqrt(n))
    p = 2 * (1 - stats.norm.cdf(abs(z)))
    p_values_h1.append(p)

axes[1].hist(p_values_h1, bins=50, density=True, alpha=0.7, color='lightcoral', edgecolor='white')
axes[1].axvline(x=0.05, color='green', linewidth=2, linestyle='--', label='$\\alpha$ = 0.05')
axes[1].set_xlabel('p-value', fontsize=12)
axes[1].set_ylabel('Density', fontsize=12)
axes[1].set_title(f'p-value Distribution under $H_1$ ($\\mu$ = {mu1})', fontsize=12)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("pvalue_distribution.png", dpi=150, bbox_inches="tight")
plt.show()

# 第1種・第2種の過誤率
alpha = 0.05
type1 = np.mean(np.array(p_values_h0) < alpha)
power = np.mean(np.array(p_values_h1) < alpha)
print(f"Type I Error Rate: {type1:.4f} (expected: {alpha})")
print(f"Power: {power:.4f}")

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

np.random.seed(42)

# 多重検定のシミュレーション
n_tests_list = [1, 5, 10, 20, 50, 100]
alpha = 0.05
n_simulations = 10000

# FWERの理論値と実験値
fwer_theory = [1 - (1 - alpha) ** m for m in n_tests_list]
fwer_empirical = []

for m in n_tests_list:
    false_positives = 0
    for _ in range(n_simulations):
        # m個の独立なz検定（すべてH0が真）
        z_values = np.random.normal(0, 1, m)
        p_values = 2 * (1 - stats.norm.cdf(np.abs(z_values)))
        if np.any(p_values < alpha):
            false_positives += 1
    fwer_empirical.append(false_positives / n_simulations)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(n_tests_list, fwer_theory, 'r-o', linewidth=2, markersize=8, label='Theoretical FWER')
ax.plot(n_tests_list, fwer_empirical, 'b--s', linewidth=2, markersize=8, label='Empirical FWER')
ax.axhline(y=alpha, color='green', linestyle=':', linewidth=2, label=f'$\\alpha$ = {alpha}')
ax.set_xlabel('Number of tests (m)', fontsize=12)
ax.set_ylabel('Family-Wise Error Rate', fontsize=12)
ax.set_title('Multiple Testing Problem: FWER increases with m', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.savefig("multiple_testing_fwer.png", dpi=150, bbox_inches="tight")
plt.show()

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

np.random.seed(42)

# シナリオ: 100個の検定のうち10個が真に有意
m = 100
m_true = 10  # 本当に効果があるもの
alpha = 0.05

# p値の生成
# 真にH0のもの: 一様分布
p_null = np.random.uniform(0, 1, m - m_true)
# 真にH1のもの: 小さなp値（効果あり）
p_alt = np.random.beta(1, 20, m_true)  # 0に偏った分布

p_values = np.concatenate([p_null, p_alt])
is_true_alt = np.concatenate([np.zeros(m - m_true), np.ones(m_true)]).astype(bool)

# ソートのインデックス
sort_idx = np.argsort(p_values)
p_sorted = p_values[sort_idx]
true_alt_sorted = is_true_alt[sort_idx]

# 補正なし
reject_none = p_values < alpha

# ボンフェローニ補正
p_bonf = np.minimum(p_values * m, 1.0)
reject_bonf = p_bonf < alpha

# BH法（Benjamini-Hochberg）
rank = np.arange(1, m + 1)
bh_threshold = rank / m * alpha
reject_bh_sorted = p_sorted <= bh_threshold
# 最大のkを見つける
if np.any(reject_bh_sorted):
    max_k = np.max(np.where(reject_bh_sorted)[0])
    reject_bh_sorted_final = np.zeros(m, dtype=bool)
    reject_bh_sorted_final[:max_k + 1] = True
else:
    reject_bh_sorted_final = np.zeros(m, dtype=bool)
# 元の順序に戻す
reject_bh = np.zeros(m, dtype=bool)
reject_bh[sort_idx] = reject_bh_sorted_final

# 結果の表示
def calc_metrics(reject, is_true_alt):
    TP = np.sum(reject & is_true_alt)
    FP = np.sum(reject & ~is_true_alt)
    FN = np.sum(~reject & is_true_alt)
    TN = np.sum(~reject & ~is_true_alt)
    return TP, FP, FN, TN

print(f"{'Method':<20s} {'TP':>4s} {'FP':>4s} {'FN':>4s} {'TN':>4s} {'FDR':>8s} {'Power':>8s}")
print("-" * 60)
for name, reject in [('No correction', reject_none),
                      ('Bonferroni', reject_bonf),
                      ('BH (FDR)', reject_bh)]:
    TP, FP, FN, TN = calc_metrics(reject, is_true_alt)
    fdr = FP / max(FP + TP, 1)
    power = TP / max(TP + FN, 1)
    print(f"{name:<20s} {TP:>4d} {FP:>4d} {FN:>4d} {TN:>4d} {fdr:>8.3f} {power:>8.3f}")

# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for ax, (name, reject) in zip(axes, [('No Correction', reject_none),
                                       ('Bonferroni', reject_bonf),
                                       ('BH (FDR)', reject_bh)]):
    colors = []
    for r, t in zip(reject, is_true_alt):
        if r and t:
            colors.append('green')    # 真陽性
        elif r and not t:
            colors.append('red')      # 偽陽性
        elif not r and t:
            colors.append('orange')   # 偽陰性
        else:
            colors.append('gray')     # 真陰性

    ax.scatter(range(m), -np.log10(p_values), c=colors, s=30, alpha=0.7)
    ax.axhline(y=-np.log10(alpha), color='blue', linestyle='--', label=f'$\\alpha$ = {alpha}')
    ax.set_xlabel('Test index', fontsize=10)
    ax.set_ylabel('$-\\log_{10}(p)$', fontsize=10)
    ax.set_title(name, fontsize=12)
    ax.grid(True, alpha=0.3)

    # カスタム凡例
    from matplotlib.patches import Patch
    legend_elements = [Patch(facecolor='green', label='True Positive'),
                       Patch(facecolor='red', label='False Positive'),
                       Patch(facecolor='orange', label='False Negative'),
                       Patch(facecolor='gray', label='True Negative')]
    ax.legend(handles=legend_elements, fontsize=8, loc='upper right')

plt.tight_layout()
plt.savefig("multiple_testing_correction.png", dpi=150, bbox_inches="tight")
plt.show()

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

# 異なる検定統計量の値に対するp値の変化
z_values = np.arange(0, 4.1, 0.5)

fig, axes = plt.subplots(2, 4, figsize=(18, 8))
axes = axes.ravel()

z_grid = np.linspace(-4, 4, 1000)
pdf = stats.norm.pdf(z_grid)

for i, z_obs in enumerate(z_values):
    if i >= 8:
        break
    ax = axes[i]
    ax.plot(z_grid, pdf, 'k-', linewidth=1.5)

    # p値の領域を塗る（両側検定）
    if z_obs > 0:
        mask_right = z_grid >= z_obs
        mask_left = z_grid <= -z_obs
        ax.fill_between(z_grid[mask_right], pdf[mask_right], alpha=0.4, color='red')
        ax.fill_between(z_grid[mask_left], pdf[mask_left], alpha=0.4, color='red')

    p_val = 2 * (1 - stats.norm.cdf(z_obs))
    ax.set_title(f'z = {z_obs:.1f}, p = {p_val:.4f}', fontsize=10)
    ax.set_xlim(-4, 4)
    ax.set_ylim(0, 0.45)
    ax.grid(True, alpha=0.3)

    # 有意かどうか
    if p_val < 0.05:
        ax.set_facecolor('#fff0f0')
    else:
        ax.set_facecolor('#f0fff0')

plt.suptitle('p-value as area under the curve (Two-sided test)', fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig("pvalue_intuition.png", dpi=150, bbox_inches="tight")
plt.show()