Module `neuralops.analysis.comparison.statistical_tests`

Complete statistical tests for model comparison – publication‑grade.

References -----= - Demšar (2006): Statistical Comparisons of Classifiers over Multiple Data Sets (JMLR) - Benavoli et al. (2017): Time for a Change – Bayesian Comparison of Classifiers (JMLR) - García & Herrera (2008): An Extension on Statistical Comparisons of Classifiers - García et al. (2010): Advanced nonparametric tests for multiple comparisons - Derrac et al. (2011): A practical tutorial on the use of nonparametric tests - Iman & Davenport (1980): Approximations of the critical region for Friedman - Quade (1979): Using weighted rankings in analysis of complete blocks - Cliff (1993): Dominance statistics – ordinal analyses to answer ordinal questions - Vargha & Delaney (2000): Critique and improvement of CL common language effect size - Cohen (1988): Statistical Power Analysis for the Behavioral Sciences - Hedges & Olkin (1985): Statistical Methods for Meta-Analysis - Nemenyi (1963): Distribution-free Multiple Comparisons - Dunn (1961): Multiple Comparisons among Means - Holm (1979): A simple sequentially rejective multiple test procedure - Hochberg (1988): A sharper Bonferroni procedure - Hommel (1988): A stagewise rejective multiple test procedure - Holland & Copenhaver (1987): An improved sequentially rejective Bonferroni test - Finner (1993): On a monotonicity problem in step-down multiple test procedures - Bergmann & Hommel (1988): Improvements of general multiple test procedures

Functions

def bayesian_correlated_ttest(x: numpy.ndarray, y: numpy.ndarray, rope: float = 0.01, runs: int = 1, nsamples: int = 50000, seed: int = 42) ‑> dict[str, float]

Expand source code

def bayesian_correlated_ttest(
    x: np.ndarray,
    y: np.ndarray,
    rope: float = 0.01,
    runs: int = 1,
    nsamples: int = 50000,
    seed: int = 42,
) -> dict[str, float]:
    """
    Bayesian correlated t‑test for k‑fold cross‑validation results.

    Based on: Benavoli et al. (2017) "Time for a Change"
    ROPE = Region of Practical Equivalence

    Returns probabilities: p_left (x < y), p_rope (x ≈ y), p_right (x > y).

    Note: This is a simplified implementation using a Monte Carlo approximation.
    For exact results, use the `baycomp` Python package.
    """
    np.random.seed(seed)
    diffs = x - y
    n = len(diffs)
    mean_diff = np.mean(diffs)

    # Standard error corrected for correlation in CV (Nadeau & Bengio, 2003)
    var_diff = np.var(diffs, ddof=1)
    # For k-fold CV, the correction factor depends on the number of folds.
    # A common approximation for 10‑fold CV is n_test / n_train ≈ 1/9.
    correction = 1.0 / runs  # simplified for runs
    se = np.sqrt(var_diff / n * (1 / correction))

    # Posterior sampling (non‑informative prior)
    posterior = np.random.normal(mean_diff, se, nsamples)

    p_left = np.mean(posterior < -rope)
    p_rope = np.mean((posterior >= -rope) & (posterior <= rope))
    p_right = np.mean(posterior > rope)

    winner = "x" if p_left > p_right else "y" if p_right > p_left else "tie"

    return {
        "method": "Bayesian correlated t-test",
        "p_left": round(float(p_left), 4),
        "p_rope": round(float(p_rope), 4),
        "p_right": round(float(p_right), 4),
        "winner": winner,
        "rope": rope,
        "mean_diff": round(float(mean_diff), 6),
    }

Bayesian correlated t‑test for k‑fold cross‑validation results.

Based on: Benavoli et al. (2017) "Time for a Change" ROPE = Region of Practical Equivalence

Returns probabilities: p_left (x < y), p_rope (x ≈ y), p_right (x > y).

Note: This is a simplified implementation using a Monte Carlo approximation. For exact results, use the baycomp Python package.

def bayesian_signed_rank_test(x: numpy.ndarray, y: numpy.ndarray, rope: float = 0.01, nsamples: int = 50000, seed: int = 42) ‑> dict[str, float]

Expand source code

def bayesian_signed_rank_test(
    x: np.ndarray,
    y: np.ndarray,
    rope: float = 0.01,
    nsamples: int = 50000,
    seed: int = 42,
) -> dict[str, float]:
    """
    Bayesian signed‑rank test – Bayesian equivalent of Wilcoxon.

    Based on: Benavoli et al. (2014, 2017)
    Simplified implementation using Dirichlet process prior approximation.
    """
    np.random.seed(seed)
    diffs = x - y
    n = len(diffs)

    # Exclude zeros for rank computation
    nonzero = diffs[diffs != 0]
    n_nonzero = len(nonzero)
    if n_nonzero == 0:
        return {"p_left": 0.0, "p_rope": 1.0, "p_right": 0.0, "winner": "tie"}

    # Signed rank sum (Wilcoxon style)
    abs_diffs = np.abs(nonzero)
    ranks = stats.rankdata(abs_diffs)
    np.sum(ranks[nonzero > 0])

    # Dirichlet posterior for the probability that a positive difference occurs
    # Using the method from Benavoli et al.
    # Prior: symmetric Dirichlet
    alpha = np.ones(2)  # uniform prior
    # Counts
    n_pos = np.sum(diffs > 0)
    n_neg = np.sum(diffs < 0)

    # Posterior samples
    posterior = np.random.dirichlet(alpha + np.array([n_pos, n_neg]), nsamples)
    theta = posterior[:, 0]  # P(X > Y)

    # Probabilities
    p_left = np.mean(theta < 0.5 - rope)
    p_rope = np.mean((theta >= 0.5 - rope) & (theta <= 0.5 + rope))
    p_right = np.mean(theta > 0.5 + rope)

    winner = "y" if p_left > p_right else "x" if p_right > p_left else "tie"

    return {
        "method": "Bayesian signed-rank",
        "p_left": round(float(p_left), 4),
        "p_rope": round(float(p_rope), 4),
        "p_right": round(float(p_right), 4),
        "winner": winner,
        "rope": rope,
        "n": int(n),
        "n_nonzero": int(n_nonzero),
    }

Bayesian signed‑rank test – Bayesian equivalent of Wilcoxon.

Based on: Benavoli et al. (2014, 2017) Simplified implementation using Dirichlet process prior approximation.

def bergmann_hommel_posthoc(avg_ranks: numpy.ndarray, n_datasets: int, k: int = None, alpha: float = 0.05) ‑> pandas.DataFrame

Expand source code

def bergmann_hommel_posthoc(
    avg_ranks: np.ndarray, n_datasets: int, k: int = None, alpha: float = 0.05
) -> pd.DataFrame:
    """
    Bergmann-Hommel post-hoc – the most powerful post-hoc for Friedman test.
    Uses the logical structure of hypotheses for dynamic correction.

    References: Bergmann & Hommel (1988), García & Herrera (2008)
    """
    k = len(avg_ranks) if k is None else k
    N = n_datasets

    # Standard error of rank difference
    se = np.sqrt(k * (k + 1) / (6 * N))

    # Compute all pairwise z-statistics
    z_matrix = np.zeros((k, k))
    for i in range(k):
        for j in range(k):
            if i != j:
                z_matrix[i, j] = abs(avg_ranks[i] - avg_ranks[j]) / se

    # Standard normal two-tailed p-values
    p_matrix = 2 * (1 - norm.cdf(z_matrix))

    # Bergmann-Hommel dynamic correction
    # Number of simultaneous hypotheses
    M = k * (k - 1) // 2
    # Sorted p-values
    all_p = []
    for i in range(k):
        for j in range(i + 1, k):
            all_p.append(p_matrix[i, j])
    all_p = np.sort(all_p)

    # Adjusted p-values (step‑up)
    adjusted_list = np.zeros(len(all_p))
    for idx in range(len(all_p)):
        adjusted_list[idx] = min(all_p[idx] * M / (idx + 1), 1.0)

    # Reconstruct matrix
    adj_matrix = np.ones((k, k))
    idx = 0
    for i in range(k):
        for j in range(i + 1, k):
            adj_matrix[i, j] = adj_matrix[j, i] = adjusted_list[idx]
            idx += 1

    return pd.DataFrame(adj_matrix)

Bergmann-Hommel post-hoc – the most powerful post-hoc for Friedman test. Uses the logical structure of hypotheses for dynamic correction.

References: Bergmann & Hommel (1988), García & Herrera (2008)

def bonferroni_dunn_posthoc(p_values: numpy.ndarray | None = None, k: int = 0, n: int = 0, alpha: float = 0.05) ‑> dict[str, typing.Any]

Expand source code

def bonferroni_dunn_posthoc(
    p_values: np.ndarray | None = None, k: int = 0, n: int = 0, alpha: float = 0.05
) -> dict[str, Any]:
    """
    Bonferroni-Dunn post-hoc: control method vs. all others.

    Adjusted alpha = alpha / (k-1)

    Reference: Dunn (1961)
    """
    return {"adjusted_alpha": alpha / (k - 1), "method": "Bonferroni-Dunn"}

Bonferroni-Dunn post-hoc: control method vs. all others.

Adjusted alpha = alpha / (k-1)

Reference: Dunn (1961)

def cliffs_delta(a: numpy.ndarray, b: numpy.ndarray) ‑> dict[str, float]

Expand source code

def cliffs_delta(a: np.ndarray, b: np.ndarray) -> dict[str, float]:
    """
    Cliff's delta – non‑parametric effect size for ordinal data.
    Measures the probability that a randomly selected value from group A
    exceeds a randomly selected value from group B.

    delta = (#(a > b) - #(a < b)) / (n_a * n_b)

    Reference: Cliff (1993)
    """
    n_a, n_b = len(a), len(b)
    gt, lt = 0, 0
    for ai in a:
        gt += np.sum(ai > b)
        lt += np.sum(ai < b)
    delta = (gt - lt) / (n_a * n_b) if n_a * n_b > 0 else 0.0
    # Interpretation: |delta| < 0.147 = negligible, < 0.33 = small, < 0.474 = medium, ≥ 0.474 = large
    d_abs = abs(delta)
    magnitude = (
        "negligible"
        if d_abs < 0.147
        else "small"
        if d_abs < 0.330
        else "medium"
        if d_abs < 0.474
        else "large"
    )
    return {"cliffs_delta": round(delta, 4), "magnitude": magnitude}

Cliff's delta – non‑parametric effect size for ordinal data. Measures the probability that a randomly selected value from group A exceeds a randomly selected value from group B.

delta = (#(a > b) - #(a < b)) / (n_a * n_b)

Reference: Cliff (1993)

def cohens_d(a: numpy.ndarray, b: numpy.ndarray, paired: bool = True) ‑> dict[str, float]

Expand source code

def cohens_d(a: np.ndarray, b: np.ndarray, paired: bool = True) -> dict[str, float]:
    """
    Cohen's d – standardised mean difference.

    Reference: Cohen (1988)
    """
    diff = np.mean(a) - np.mean(b)
    if paired:
        # Use standard deviation of differences
        sd = np.std(a - b, ddof=1)
    else:
        # Pooled standard deviation
        n1, n2 = len(a), len(b)
        sd = np.sqrt(((n1 - 1) * np.var(a, ddof=1) + (n2 - 1) * np.var(b, ddof=1)) / (n1 + n2 - 2))
    d = diff / sd if sd > 0 else 0.0
    # Interpretation: small=0.2, medium=0.5, large=0.8
    magnitude = (
        "negligible"
        if abs(d) < 0.2
        else ("small" if abs(d) < 0.5 else ("medium" if abs(d) < 0.8 else "large"))
    )
    return {"cohens_d": round(d, 4), "magnitude": magnitude}

Cohen's d – standardised mean difference.

Reference: Cohen (1988)

def compute_aligned_ranks(scores_matrix: numpy.ndarray) ‑> numpy.ndarray

Expand source code

def compute_aligned_ranks(scores_matrix: np.ndarray) -> np.ndarray:
    """
    Compute aligned ranks (Hodges-Lehmann / Friedman Aligned Ranks).

    The aligned rank for each observation is its rank among ALL observations
    after subtracting the row mean.  This allows comparison across rows.

    Reference: Hodges & Lehmann (1962)
    """
    n_rows, n_cols = scores_matrix.shape
    # Compute mean performance per dataset
    row_means = np.mean(scores_matrix, axis=1, keepdims=True)
    # Subtract row mean (align)
    aligned = scores_matrix - row_means
    # Rank all aligned values together
    flat = aligned.flatten()
    flat_ranks = stats.rankdata(flat, method="average")
    # Reshape back
    aligned_ranks = flat_ranks.reshape(n_rows, n_cols)
    return aligned_ranks

Compute aligned ranks (Hodges-Lehmann / Friedman Aligned Ranks).

The aligned rank for each observation is its rank among ALL observations after subtracting the row mean. This allows comparison across rows.

Reference: Hodges & Lehmann (1962)

def compute_ranks(scores_matrix: numpy.ndarray) ‑> numpy.ndarray

Expand source code

def compute_ranks(scores_matrix: np.ndarray) -> np.ndarray:
    """
    Compute average ranks per column (algorithm) across rows (datasets).
    Within each row, the best algorithm gets rank 1 (highest score = best).

    Args:
        scores_matrix: shape (n_datasets, n_algorithms) – higher is better

    Returns:
        ranks: shape (n_algorithms,) – average rank
    """
    n_rows, n_cols = scores_matrix.shape
    ranks = np.zeros(n_cols)
    for i in range(n_rows):
        row = scores_matrix[i]
        # Use average ranking for ties
        row_ranks = stats.rankdata(-row, method="average")  # negate so highest=rank 1
        ranks += row_ranks
    return ranks / n_rows

Compute average ranks per column (algorithm) across rows (datasets). Within each row, the best algorithm gets rank 1 (highest score = best).

Args
-----=
scores_matrix: shape (n_datasets, n_algorithms) – higher is better
Returns
-----=
ranks: shape (n_algorithms,) – average rank

def critical_difference(k: int, N: int, alpha: float = 0.05) ‑> float

Expand source code

def critical_difference(k: int, N: int, alpha: float = 0.05) -> float:
    """
    Nemenyi Critical Difference (CD).

    CD = q_alpha * sqrt(k * (k + 1) / (6 * N))

    Two algorithms are significantly different if their average rank
    difference exceeds CD.

    Reference: Nemenyi (1963), Demšar (2006)
    """
    q_alpha = studentized_range.ppf(1 - alpha, k, np.inf)
    cd = q_alpha * np.sqrt(k * (k + 1) / (6 * N))
    return cd

Nemenyi Critical Difference (CD).

CD = q_alpha * sqrt(k * (k + 1) / (6 * N))

Two algorithms are significantly different if their average rank difference exceeds CD.

Reference: Nemenyi (1963), Demšar (2006)

def finner_correction(p_values: list[float] | numpy.ndarray) ‑> numpy.ndarray

Expand source code

def finner_correction(p_values: list[float] | np.ndarray) -> np.ndarray:
    """
    Finner step‑down correction – most powerful but less common.

    Reference: Finner (1993)
    """
    p = np.array(p_values, dtype=float)
    m = len(p)
    order = np.argsort(p)
    adjusted = np.ones(m)
    for i, idx in enumerate(order):
        adjusted[idx] = 1 - (1 - p[idx]) ** (m - i)
        if i > 0:
            adjusted[idx] = max(adjusted[idx], adjusted[order[i - 1]])
    return adjusted

Finner step‑down correction – most powerful but less common.

Reference: Finner (1993)

def friedman_aligned_ranks_test(scores_matrix: numpy.ndarray) ‑> dict[str, float]

Expand source code

def friedman_aligned_ranks_test(scores_matrix: np.ndarray) -> dict[str, float]:
    """
    Friedman Aligned Ranks test – more powerful than the standard Friedman test.

    Reference: García et al. (2010), Hodges & Lehmann (1962)
    """
    n, k = scores_matrix.shape
    aligned_ranks = compute_aligned_ranks(scores_matrix)

    # Sum of aligned ranks per algorithm
    R_j = np.sum(aligned_ranks, axis=0)
    total = n * k
    # Expected rank sum under H0
    expected = k * n * (total + 1) / (2 * k)
    # Test statistic
    T = ((k - 1) * np.sum((R_j - expected) ** 2)) / (
        np.sum((aligned_ranks.flatten() - (total + 1) / 2) ** 2)
    )

    p_value = 1 - chi2_dist.cdf(T, k - 1)
    significant = bool(p_value < 0.05)

    return {
        "method": "Friedman Aligned Ranks",
        "statistic": float(T),
        "p_value": float(p_value),
        "significant": significant,
        "aligned_rank_sums": R_j.tolist(),
    }

Friedman Aligned Ranks test – more powerful than the standard Friedman test.

Reference: García et al. (2010), Hodges & Lehmann (1962)

def friedman_test(scores_matrix: numpy.ndarray) ‑> dict[str, float]

Expand source code

def friedman_test(scores_matrix: np.ndarray) -> dict[str, float]:
    """
    Friedman test – the classic omnibus test for comparing k algorithms
    over N datasets.  H0: all algorithms perform equivalently.
    """
    n, k = scores_matrix.shape
    chi2_stat, p_value = friedmanchisquare(*[scores_matrix[:, j] for j in range(k)])

    # Iman-Davenport correction (recommended over plain Friedman)
    iman_davenport = (
        ((n - 1) * chi2_stat) / (n * (k - 1) - chi2_stat) if chi2_stat < n * (k - 1) else np.inf
    )
    p_iman = 1 - f_dist.cdf(iman_davenport, k - 1, (k - 1) * (n - 1))

    # Average ranks
    avg_ranks = compute_ranks(scores_matrix)

    return {
        "method": "Friedman",
        "chi2": float(chi2_stat),
        "p_value": float(p_value),
        "iman_davenport_F": float(iman_davenport),
        "p_value_iman_davenport": float(p_iman),
        "avg_ranks": avg_ranks.tolist(),
        "significant": bool(p_iman < 0.05),
    }

Friedman test – the classic omnibus test for comparing k algorithms over N datasets. H0: all algorithms perform equivalently.

def glass_rank_biserial(a: numpy.ndarray, b: numpy.ndarray) ‑> dict[str, float]

Expand source code

def glass_rank_biserial(a: np.ndarray, b: np.ndarray) -> dict[str, float]:
    """
    Glass rank-biserial correlation – effect size for Mann‑Whitney / Wilcoxon.
    r_rb = 2 * (A12 - 0.5)

    Reference: Glass (1966)
    """
    a12 = vargha_delaney_a12(a, b)["vargha_delaney_A12"]
    r_rb = 2 * (a12 - 0.5)
    return {"glass_rank_biserial": round(r_rb, 4)}

Glass rank-biserial correlation – effect size for Mann‑Whitney / Wilcoxon. r_rb = 2 * (A12 - 0.5)

Reference: Glass (1966)

def hedges_g(a: numpy.ndarray, b: numpy.ndarray, paired: bool = True) ‑> dict[str, float]

Expand source code

def hedges_g(a: np.ndarray, b: np.ndarray, paired: bool = True) -> dict[str, float]:
    """
    Hedges' g – bias‑corrected Cohen's d for small samples.

    Reference: Hedges & Olkin (1985)
    """
    d = cohens_d(a, b, paired)["cohens_d"]
    n = len(a) + len(b) if not paired else len(a)
    # Correction factor
    J = 1 - 3 / (4 * (n - 1) - 1) if n > 2 else 1.0
    g = d * J
    magnitude = (
        "negligible"
        if abs(g) < 0.2
        else ("small" if abs(g) < 0.5 else ("medium" if abs(g) < 0.8 else "large"))
    )
    return {"hedges_g": round(g, 4), "magnitude": magnitude}

Hedges' g – bias‑corrected Cohen's d for small samples.

Reference: Hedges & Olkin (1985)

def hochberg_correction(p_values: list[float] | numpy.ndarray) ‑> numpy.ndarray

Expand source code

def hochberg_correction(p_values: list[float] | np.ndarray) -> np.ndarray:
    """
    Hochberg step‑up correction – generally more powerful than Holm.

    Reference: Hochberg (1988)
    """
    p = np.array(p_values, dtype=float)
    m = len(p)
    order = np.argsort(p)[::-1]  # descending
    adjusted = np.ones(m)
    for i, idx in enumerate(order):
        adjusted[idx] = min(p[idx] * (i + 1), 1.0)
        if i > 0:
            adjusted[idx] = min(adjusted[idx], adjusted[order[i - 1]])
    return adjusted

Hochberg step‑up correction – generally more powerful than Holm.

Reference: Hochberg (1988)

def holland_copenhaver_correction(p_values: list[float] | numpy.ndarray) ‑> numpy.ndarray

Expand source code

def holland_copenhaver_correction(
    p_values: list[float] | np.ndarray,
) -> np.ndarray:
    """
    Holland & Copenhaver correction – improved sequentially rejective Bonferroni.

    Reference: Holland & Copenhaver (1987)
    """
    p = np.array(p_values, dtype=float)
    m = len(p)
    order = np.argsort(p)
    adjusted = np.ones(m)
    for i, idx in enumerate(order):
        adjusted[idx] = 1 - (1 - p[idx]) ** (m - i)
        if i > 0:
            adjusted[idx] = max(adjusted[idx], adjusted[order[i - 1]])
    return adjusted

Holland & Copenhaver correction – improved sequentially rejective Bonferroni.

Reference: Holland & Copenhaver (1987)

def holm_correction(p_values: list[float] | numpy.ndarray) ‑> numpy.ndarray

Expand source code

def holm_correction(p_values: list[float] | np.ndarray) -> np.ndarray:
    """
    Holm step-down correction – controls family-wise error rate.
    More powerful than Bonferroni.

    Reference: Holm (1979)
    """
    p = np.array(p_values, dtype=float)
    m = len(p)
    order = np.argsort(p)
    adjusted = np.ones(m)
    for i, idx in enumerate(order):
        adjusted[idx] = min(p[idx] * (m - i), 1.0)
        if i > 0:
            adjusted[idx] = max(adjusted[idx], adjusted[order[i - 1]])
    return adjusted

Holm step-down correction – controls family-wise error rate. More powerful than Bonferroni.

Reference: Holm (1979)

def hommel_correction(p_values: list[float] | numpy.ndarray) ‑> numpy.ndarray

Expand source code

def hommel_correction(p_values: list[float] | np.ndarray) -> np.ndarray:
    """
    Hommel correction – stagewise rejective procedure.

    Reference: Hommel (1988)
    """
    p = np.array(p_values, dtype=float)
    m = len(p)
    order = np.argsort(p)
    adjusted = np.ones(m)
    for i in range(m):
        j = order[i]
        candidates = [p[j]]
        for r in range(1, m - i + 1):
            candidates.append(r * p[j] / (r - 0.5) if r > 1 else p[j])
        adjusted[j] = min(1.0, min(candidates))
    return adjusted

Hommel correction – stagewise rejective procedure.

Reference: Hommel (1988)

def mcnemar_test(y_true: numpy.ndarray, pred_a: numpy.ndarray, pred_b: numpy.ndarray) ‑> dict[str, typing.Any]

Expand source code

def mcnemar_test(y_true: np.ndarray, pred_a: np.ndarray, pred_b: np.ndarray) -> dict[str, Any]:
    """
    McNemar's test for comparing two classifiers on a single test set.
    Uses a 2x2 contingency table of correct/incorrect predictions.

    Reference: McNemar (1947)
    """
    correct_a = pred_a == y_true
    correct_b = pred_b == y_true

    # Count discordant pairs
    n01 = np.sum((correct_a == 0) & (correct_b == 1))  # B wins, A loses
    n10 = np.sum((correct_a == 1) & (correct_b == 0))  # A wins, B loses

    if n01 + n10 == 0:
        return {"statistic": 0.0, "p_value": 1.0, "significant": False}

    # Continuity correction
    chi2 = (abs(n01 - n10) - 1) ** 2 / (n01 + n10) if n01 + n10 > 0 else 0
    p_value = 1 - chi2_dist.cdf(chi2, 1)

    return {
        "method": "McNemar",
        "chi2": float(chi2),
        "p_value": float(p_value),
        "significant": bool(p_value < 0.05),
        "n01_A_wrong_B_right": int(n01),
        "n10_A_right_B_wrong": int(n10),
    }

McNemar's test for comparing two classifiers on a single test set. Uses a 2x2 contingency table of correct/incorrect predictions.

Reference: McNemar (1947)

def nemenyi_posthoc(avg_ranks: numpy.ndarray, n_datasets: int, alpha: float = 0.05) ‑> pandas.DataFrame

Expand source code

def nemenyi_posthoc(avg_ranks: np.ndarray, n_datasets: int, alpha: float = 0.05) -> pd.DataFrame:
    """
    Nemenyi post-hoc test: all vs. all pairwise comparisons.

    Returns a DataFrame of p-values (approximated from CD).
    """
    k = len(avg_ranks)
    cd = critical_difference(k, n_datasets, alpha)
    # Normalised difference matrix: |rank_i - rank_j| / CD
    diff_matrix = np.zeros((k, k))
    for i in range(k):
        for j in range(k):
            diff_matrix[i, j] = abs(avg_ranks[i] - avg_ranks[j]) / cd
    names = [f"Alg {i + 1}" for i in range(k)]
    return pd.DataFrame(diff_matrix, index=names, columns=names)

Nemenyi post-hoc test: all vs. all pairwise comparisons.

Returns a DataFrame of p-values (approximated from CD).

def quade_test(scores_matrix: numpy.ndarray) ‑> dict[str, float]

Expand source code

def quade_test(scores_matrix: np.ndarray) -> dict[str, float]:
    """
    Quade test – weighted ranking test that gives more weight to datasets
    where larger differences between algorithms exist.

    References: Quade (1979), Conover (1999)
    """
    n, k = scores_matrix.shape

    # Step 1: Compute range per dataset (row)
    row_ranges = np.max(scores_matrix, axis=1) - np.min(scores_matrix, axis=1)
    # Handle zero-range rows
    row_ranges[row_ranges == 0] = 1e-10

    # Step 2: Rank the ranges
    range_ranks = stats.rankdata(row_ranges, method="average")

    # Step 3: Compute ranks within each row, then weight by range rank
    Q = np.zeros(k)
    S = np.zeros(k)
    row_ranks_matrix = np.zeros_like(scores_matrix)
    for i in range(n):
        row_ranks = stats.rankdata(-scores_matrix[i], method="average")
        row_ranks_matrix[i] = row_ranks
        for j in range(k):
            weighted = range_ranks[i] * (row_ranks[j] - (k + 1) / 2)
            Q[j] += weighted
    S = np.sum(Q**2)

    # Step 4: Total sum of squares
    A = np.sum((row_ranks_matrix - (k + 1) / 2) ** 2 * range_ranks[:, np.newaxis] ** 2)

    # Step 5: Test statistic
    F_quade = ((n - 1) * S) / (A - S) if A > S else np.inf
    p_value = 1 - f_dist.cdf(F_quade, k - 1, (k - 1) * (n - 1))

    return {
        "method": "Quade",
        "statistic": float(F_quade),
        "p_value": float(p_value),
        "significant": bool(p_value < 0.05),
    }

Quade test – weighted ranking test that gives more weight to datasets where larger differences between algorithms exist.

References: Quade (1979), Conover (1999)

def run_full_statistical_analysis(scores_matrix: numpy.ndarray, algorithm_names: list[str] | None = None, alpha: float = 0.05, generate_posthoc: bool = True) ‑> dict[str, typing.Any]

Expand source code

def run_full_statistical_analysis(
    scores_matrix: np.ndarray,
    algorithm_names: list[str] | None = None,
    alpha: float = 0.05,
    generate_posthoc: bool = True,
) -> dict[str, Any]:
    """
    Run the complete statistical pipeline recommended by Demšar (2006)
    and extended by García & Herrera (2008, 2010):

    1. Friedman test (with Iman-Davenport correction)
    2. Friedman Aligned Ranks test
    3. Quade test
    4. If significant: Nemenyi post‑hoc + Bonferroni-Dunn + Holm + Hochberg
    5. Average ranks + ranking

    Args:
        scores_matrix: shape (n_datasets, n_algorithms), higher = better
        algorithm_names: list of algorithm name strings
        alpha: significance level
        generate_posthoc: whether to compute post‑hoc tests

    Returns:
        Dictionary with all results
    """
    n, k = scores_matrix.shape
    if algorithm_names is None:
        algorithm_names = [f"Algorithm {i + 1}" for i in range(k)]

    # Omnibus tests
    friedman = friedman_test(scores_matrix)
    friedman_aligned = friedman_aligned_ranks_test(scores_matrix)
    quade = quade_test(scores_matrix)

    result = {
        "n_datasets": n,
        "n_algorithms": k,
        "algorithm_names": algorithm_names,
        "alpha": alpha,
        "friedman": friedman,
        "friedman_aligned_ranks": friedman_aligned,
        "quade": quade,
        "ranking": [
            (name, round(rank, 3)) for name, rank in zip(algorithm_names, friedman["avg_ranks"])
        ],
        "best_algorithm": algorithm_names[np.argmin(friedman["avg_ranks"])],
    }

    if generate_posthoc and friedman["significant"]:
        avg_ranks = np.array(friedman["avg_ranks"])
        result["nemenyi"] = {
            "critical_distance": round(critical_difference(k, n, alpha), 4),
            "matrix": nemenyi_posthoc(avg_ranks, n, alpha),
        }
        result["bonferroni_dunn"] = bonferroni_dunn_posthoc(k=k, n=n, alpha=alpha)
        result["bergmann_hommel"] = bergmann_hommel_posthoc(avg_ranks, n, k, alpha)

    return result

Run the complete statistical pipeline recommended by Demšar (2006) and extended by García & Herrera (2008, 2010):

Friedman test (with Iman-Davenport correction)
Friedman Aligned Ranks test
Quade test
If significant: Nemenyi post‑hoc + Bonferroni-Dunn + Holm + Hochberg
Average ranks + ranking

Args
-----=
scores_matrix: shape (n_datasets, n_algorithms), higher = better
algorithm_names: list of algorithm name strings
alpha: significance level
generate_posthoc: whether to compute post‑hoc tests

Returns -----= Dictionary with all results

def sign_test(scores_a: numpy.ndarray, scores_b: numpy.ndarray) ‑> dict[str, typing.Any]

Expand source code

def sign_test(scores_a: np.ndarray, scores_b: np.ndarray) -> dict[str, Any]:
    """
    Simple sign test – counts wins/losses with binomial distribution.
    Less powerful than Wilcoxon but assumption‑free.
    """
    diffs = scores_a - scores_b
    wins = np.sum(diffs > 0)
    losses = np.sum(diffs < 0)
    ties = np.sum(diffs == 0)
    n = wins + losses
    p_value = 2 * stats.binom.sf(max(wins, losses) - 1, n, 0.5) if n > 1 else 1.0
    return {
        "method": "Sign test",
        "wins_a": int(wins),
        "wins_b": int(losses),
        "ties": int(ties),
        "p_value": float(p_value),
        "significant": bool(p_value < 0.05),
    }

Simple sign test – counts wins/losses with binomial distribution. Less powerful than Wilcoxon but assumption‑free.

def vargha_delaney_a12(a: numpy.ndarray, b: numpy.ndarray) ‑> dict[str, float]

Expand source code

def vargha_delaney_a12(a: np.ndarray, b: np.ndarray) -> dict[str, float]:
    """
    Vargha-Delaney's A12 – stochastic superiority measure.
    A12 = P(X > Y) + 0.5 * P(X = Y)

    Interpretation: 0.56=small, 0.64=medium, 0.71=large

    Reference: Vargha & Delaney (2000)
    """
    n_a, n_b = len(a), len(b)
    gt, eq = 0, 0
    for ai in a:
        gt += np.sum(ai > b)
        eq += np.sum(ai == b)
    a12 = (gt + 0.5 * eq) / (n_a * n_b) if n_a * n_b > 0 else 0.5
    # Normalise so that >= 0.5 (A is better if A12 > 0.5)
    if a12 < 0.5:
        a12 = 1 - a12
    magnitude = (
        "negligible"
        if a12 < 0.56
        else "small"
        if a12 < 0.64
        else "medium"
        if a12 < 0.71
        else "large"
    )
    return {"vargha_delaney_A12": round(a12, 4), "magnitude": magnitude}

Vargha-Delaney's A12 – stochastic superiority measure. A12 = P(X > Y) + 0.5 * P(X = Y)

Interpretation: 0.56=small, 0.64=medium, 0.71=large

Reference: Vargha & Delaney (2000)

def wilcoxon_signed_rank_test(scores_a: numpy.ndarray, scores_b: numpy.ndarray) ‑> dict[str, typing.Any]

Expand source code

def wilcoxon_signed_rank_test(scores_a: np.ndarray, scores_b: np.ndarray) -> dict[str, Any]:
    """
    Wilcoxon signed-rank test – the recommended test for comparing
    two classifiers over multiple datasets (Demšar 2006).
    """
    statistic, p_value = wilcoxon(
        scores_a, scores_b, zero_method="wilcox", alternative="two-sided", method="auto"
    )
    n = len(scores_a)
    # Classification of effect (Demšar 2006 heuristic)
    z = (statistic - n * (n + 1) / 4) / np.sqrt(n * (n + 1) * (2 * n + 1) / 24)
    return {
        "method": "Wilcoxon signed-rank",
        "statistic": float(statistic),
        "p_value": float(p_value),
        "z_approx": float(z),
        "significant": bool(p_value < 0.05),
        "n": n,
    }

Wilcoxon signed-rank test – the recommended test for comparing two classifiers over multiple datasets (Demšar 2006).