import numpy as np
from scipy import stats

def bootstrap_ci(
    data: np.ndarray,
    statistic=np.mean,
    n_bootstrap: int = 10_000,
    confidence: float = 0.95,
    method: str = "percentile",
    seed: int = 42,
) -> tuple[float, float, float]:
    """Compute bootstrap confidence interval.

    Args:
        data: 1D array of observations (e.g., accuracy per seed).
        statistic: Function to compute on each resample.
        method: 'percentile' or 'bca'.

    Returns:
        (point_estimate, lower, upper)
    """
    rng = np.random.default_rng(seed)
    observed = statistic(data)
    n = len(data)

    boot_stats = np.array([
        statistic(rng.choice(data, size=n, replace=True))
        for _ in range(n_bootstrap)
    ])

    alpha = 1 - confidence
    if method == "percentile":
        lower = np.percentile(boot_stats, 100 * alpha / 2)
        upper = np.percentile(boot_stats, 100 * (1 - alpha / 2))
    elif method == "bca":
        # Bias-correction
        z0 = stats.norm.ppf(np.mean(boot_stats < observed))
        # Acceleration (jackknife)
        jackknife = np.array([
            statistic(np.delete(data, i)) for i in range(n)
        ])
        jack_mean = jackknife.mean()
        acc = np.sum((jack_mean - jackknife) ** 3) / (
            6 * np.sum((jack_mean - jackknife) ** 2) ** 1.5
        )
        # Adjusted percentiles
        z_alpha = stats.norm.ppf(alpha / 2)
        z_1alpha = stats.norm.ppf(1 - alpha / 2)
        p_lower = stats.norm.cdf(z0 + (z0 + z_alpha) / (1 - acc * (z0 + z_alpha)))
        p_upper = stats.norm.cdf(z0 + (z0 + z_1alpha) / (1 - acc * (z0 + z_1alpha)))
        lower = np.percentile(boot_stats, 100 * p_lower)
        upper = np.percentile(boot_stats, 100 * p_upper)
    else:
        raise ValueError(f"Unknown method: {method}")

    return observed, lower, upper

# Usage
accuracies = np.array([92.1, 92.5, 91.8, 92.3, 92.0, 91.9, 92.4, 92.2, 92.6, 91.7])
point, lo, hi = bootstrap_ci(accuracies, method="bca")
print(f"Accuracy: {point:.2f} [{lo:.2f}, {hi:.2f}] (95% BCa CI)")

def paired_permutation_test(
    scores_a: np.ndarray,
    scores_b: np.ndarray,
    n_permutations: int = 100_000,
    seed: int = 42,
) -> tuple[float, float]:
    """Two-sided paired permutation test.

    Returns:
        (observed_diff, p_value) where diff = mean(a) - mean(b).
    """
    rng = np.random.default_rng(seed)
    diffs = scores_a - scores_b
    observed = np.mean(diffs)

    count = 0
    for _ in range(n_permutations):
        signs = rng.choice([-1, 1], size=len(diffs))
        perm_diff = np.mean(diffs * signs)
        if abs(perm_diff) >= abs(observed):
            count += 1

    p_value = count / n_permutations
    return observed, p_value

# Usage
model_a = np.array([92.1, 92.5, 91.8, 92.3, 92.0])
model_b = np.array([91.5, 91.8, 91.2, 91.9, 91.4])
diff, p = paired_permutation_test(model_a, model_b)
print(f"Mean diff: {diff:.2f}, p-value: {p:.4f}")

from statsmodels.stats.multitest import multipletests

def correct_pvalues(
    p_values: list[float],
    method: str = "holm",  # 'bonferroni', 'holm', 'fdr_bh'
    alpha: float = 0.05,
) -> dict:
    """Apply multiple comparison correction.

    Methods:
        'bonferroni': Conservative. alpha/n per test.
        'holm': Step-down Bonferroni. Less conservative, controls FWER.
        'fdr_bh': Benjamini-Hochberg. Controls false discovery rate.
    """
    reject, corrected, _, _ = multipletests(p_values, alpha=alpha, method=method)
    return {
        "reject": reject.tolist(),
        "corrected_p": corrected.tolist(),
        "method": method,
    }

# Usage: comparing model against 4 baselines
raw_p = [0.01, 0.04, 0.03, 0.08]
result = correct_pvalues(raw_p, method="holm")
# result['reject'] -> [True, True, True, False]
# result['corrected_p'] -> [0.04, 0.08, 0.09, 0.08]

from statsmodels.stats.power import TTestPower

def compute_required_seeds(
    effect_size: float,
    alpha: float = 0.05,
    power: float = 0.8,
) -> int:
    """How many seeds needed to detect a given effect size.

    effect_size: Cohen's d = (mean_a - mean_b) / pooled_std
        Small: 0.2, Medium: 0.5, Large: 0.8

    Typical ML scenario:
        If accuracy diff ~0.5% and std ~0.3%, d = 0.5/0.3 = 1.67 (large).
        If accuracy diff ~0.2% and std ~0.3%, d = 0.2/0.3 = 0.67 (medium).
    """
    analysis = TTestPower()
    n = analysis.solve_power(
        effect_size=effect_size,
        alpha=alpha,
        power=power,
        alternative="two-sided",
    )
    return int(np.ceil(n))

# How many seeds to detect 0.5% accuracy diff with 0.3% std?
d = 0.5 / 0.3  # ~1.67
n_seeds = compute_required_seeds(d)
print(f"Need {n_seeds} seeds per model")  # ~5

# Smaller effect: 0.2% diff, 0.3% std
d_small = 0.2 / 0.3  # ~0.67
n_seeds_small = compute_required_seeds(d_small)
print(f"Need {n_seeds_small} seeds per model")  # ~20