"""Bayesian Friedman test."""
import numpy as np
from itertools import permutations
from typing import Any
from scipy.stats import f
from scipy.stats import rankdata
from compml.evaluation_scores import EvaluationScores
from compml.results.bayesian import (
BayesianFriedmanTestDecision,
BayesianFriedmanTestResult,
)
[docs]
class BayesianFriedmanTest:
"""Bayesian Friedman test for comparing multiple algorithms across multiple datasets.
This implementation follows the Bayesian Friedman test in Benavoli et al. (2015).
It uses a Dirichlet prior over the ranks of the algorithms and computes the posterior
distribution over the ranks. It can be used to test the null hypothesis of no difference
between the algorithms.
"""
def __init__(
self,
prior_strength: float = 1.0,
alpha: float = 0.05,
imprecise: bool = False,
):
"""Initialize the Bayesian Friedman test.
Parameters
----------
prior_strength : float, optional
Strength of the prior. Defaults to ``1.0``.
alpha : float, optional
Probability of wrongly rejecting equality if it holds. Defaults to ``0.05``.
imprecise : bool, optional
If ``True``, use the imprecise method of calculation. Defaults to ``False``.
"""
self.__validate_arguments(
prior_strength=prior_strength,
alpha=alpha,
imprecise=imprecise,
)
self.prior_strength = prior_strength
self.alpha = alpha
self.imprecise = imprecise
[docs]
def run(self, scores: EvaluationScores) -> BayesianFriedmanTestResult:
"""Run the Bayesian Friedman test on the given evaluation scores.
Parameters
----------
scores : EvaluationScores
The evaluation scores object containing algorithm performance over datasets.
Returns
-------
BayesianFriedmanTestResult
The result of the Bayesian Friedman test.
"""
self._validate_scores(scores)
n_datasets, n_algorithms = scores.n_datasets, scores.n_algorithms
posterior_mean_ranks, posterior_covariance, imprecise_decisions = (
None,
None,
None,
)
rank_matrix = self._compute_ranks(scores)
threshold = self._compute_decision_threshold(n_datasets, n_algorithms)
if self.imprecise:
decision, imprecise_decisions = self._evaluate_imprecise_hypothesis(
rank_matrix, threshold
)
else:
prior_mean_ranks = np.full(n_algorithms, (n_algorithms + 1.0) / 2.0)
posterior_mean_ranks = self._compute_posterior_mean_ranks(
rank_matrix, prior_mean_ranks
)
posterior_covariance = self._compute_posterior_covariance(
rank_matrix, prior_mean_ranks, posterior_mean_ranks
)
decision = self._evaluate_hypothesis(
prior_mean_ranks, posterior_mean_ranks, posterior_covariance, threshold
)
metadata = self._get_metadata()
return BayesianFriedmanTestResult(
metric=scores.metric,
procedure_name=self.__class__.__name__,
dataset_names=scores.dataset_names,
algorithm_names=scores.algorithm_names,
metadata=metadata,
decision=decision,
posterior_mean_ranks=posterior_mean_ranks,
posterior_covariance=posterior_covariance,
imprecise_decisions=imprecise_decisions,
)
def _compute_ranks(self, scores: EvaluationScores) -> np.ndarray:
"""
Compute the rank matrix for the given evaluation scores.
Parameters
----------
scores : EvaluationScores
The evaluation scores object containing algorithm performance over datasets.
Returns
-------
np.ndarray
Rank matrix of shape (n_algorithms, n_datasets), where each entry is
the rank (average for ties) of the corresponding algorithm on each dataset.
"""
mean_scores = scores.compute_means() # (n_datasets, n_algorithms)
n_datasets, _ = mean_scores.shape
ranks = np.empty_like(mean_scores, dtype=float)
for i in range(n_datasets):
ranks[i, :] = rankdata(mean_scores[i, :], method="average")
return ranks.T # (n_algorithms, n_datasets)
def _compute_decision_threshold(self, n_datasets: int, n_algorithms: int) -> float:
"""Compute the F based decision threshold rho.
Parameters
----------
n_datasets : int
Number of datasets.
n_algorithms : int
Number of algorithms.
Returns
-------
float
The F based decision threshold rho.
"""
n, m = n_datasets, n_algorithms
f_crit = float(
f.ppf(
1.0 - self.alpha,
dfn=(m - 1),
dfd=(n - m + 1),
)
)
rho = f_crit * (n - 1.0) * (m - 1.0) / (n - m + 1.0)
return rho
def _compute_posterior_mean_ranks(
self, rank_matrix: np.ndarray, prior_mean_ranks: np.ndarray
) -> np.ndarray:
"""Compute posterior mean ranks given a prior mean rank vector.
Parameters
----------
rank_matrix : np.ndarray
Rank matrix of shape (n_algorithms, n_datasets).
prior_mean_ranks : np.ndarray
Prior mean rank vector of shape (n_algorithms,).
Returns
-------
np.ndarray
Posterior mean ranks of shape (n_algorithms,).
"""
_, n_datasets = rank_matrix.shape
row_sums = np.sum(rank_matrix, axis=1)
posterior_mean_ranks = (row_sums + self.prior_strength * prior_mean_ranks) / (
self.prior_strength + n_datasets
)
return posterior_mean_ranks
def _compute_posterior_covariance(
self,
rank_matrix: np.ndarray,
prior_mean_ranks: np.ndarray,
posterior_mean_ranks: np.ndarray,
) -> np.ndarray:
"""Compute posterior covariance matrix for mean ranks.
Parameters
----------
rank_matrix : np.ndarray
Rank matrix of shape (n_algorithms, n_datasets).
prior_mean_ranks : np.ndarray
Prior mean rank vector of shape (n_algorithms,).
posterior_mean_ranks : np.ndarray
Posterior mean ranks of shape (n_algorithms,).
Returns
-------
np.ndarray
Posterior covariance matrix of shape (n_algorithms, n_algorithms).
"""
n_algorithms, n_datasets = rank_matrix.shape
s = self.prior_strength
n = n_datasets
# Covariance of Dirichlet weights over datasets
cov_w = (np.ones((n, n)) + np.eye(n)) / ((s + n) * (s + n + 1.0))
r0 = prior_mean_ranks.reshape(n_algorithms, 1)
r_mean = posterior_mean_ranks.reshape(n_algorithms, 1)
row_sums = np.sum(rank_matrix, axis=1, keepdims=True)
term1 = r0 @ r0.T * (s + 1.0)
term2 = r0 @ np.ones((1, n_datasets)) @ rank_matrix.T
term3 = row_sums @ r0.T
sigma = (
s * (term1 + term2 + term3) / ((s + n) * (s + n + 1.0))
+ rank_matrix @ cov_w @ rank_matrix.T
- r_mean @ r_mean.T
)
return sigma
def _evaluate_hypothesis(
self,
prior_mean_ranks: np.ndarray,
posterior_mean_ranks: np.ndarray,
posterior_covariance: np.ndarray,
threshold: float,
) -> BayesianFriedmanTestDecision:
"""Evaluate the hypothesis for a single prior mean rank vector.
Parameters
----------
prior_mean_ranks : np.ndarray
Prior mean ranks of shape (n_algorithms,).
posterior_mean_ranks : np.ndarray
Posterior mean ranks of shape (n_algorithms,).
posterior_covariance : np.ndarray
Posterior covariance matrix of shape (n_algorithms, n_algorithms).
threshold : float
Threshold for the F distribution.
Returns
-------
BayesianFriedmanTestDecision
DO_NOT_REJECT_NULL if the null is not rejected, REJECT_NULL otherwise.
"""
eigenvalues, _ = np.linalg.eigh(posterior_covariance)
idx = np.where(np.abs(eigenvalues) > 1e-15)[0]
m = posterior_mean_ranks.shape[0]
if idx.size == m and m > 0:
# Switch to an (m-1) dimensional subspace to avoid singularity
idx = idx[:-1]
if idx.size == 0:
# Degenerate covariance, accept if posterior mean equals prior mean
if np.allclose(posterior_mean_ranks, prior_mean_ranks):
return BayesianFriedmanTestDecision.DO_NOT_REJECT_NULL
return BayesianFriedmanTestDecision.REJECT_NULL
# Build reduced covariance and mean ranks
prior_mean_ranks_sub = prior_mean_ranks[idx]
posterior_mean_ranks_sub = posterior_mean_ranks[idx]
posterior_cov_sub = posterior_covariance[np.ix_(idx, idx)]
# Compute quadratic form
diff = posterior_mean_ranks_sub - prior_mean_ranks_sub
q_value = float(diff.T @ np.linalg.solve(posterior_cov_sub, diff))
return (
BayesianFriedmanTestDecision.DO_NOT_REJECT_NULL
if q_value < threshold
else BayesianFriedmanTestDecision.REJECT_NULL
)
def _evaluate_imprecise_hypothesis(
self,
rank_matrix: np.ndarray,
threshold: float,
) -> tuple[BayesianFriedmanTestDecision, list[BayesianFriedmanTestDecision]]:
"""Evaluate hypothesis over all permutations of prior mean ranks.
Parameters
----------
rank_matrix : np.ndarray
Rank matrix of shape (n_algorithms, n_datasets).
threshold : float
Threshold for the F distribution.
Returns
-------
tuple[BayesianFriedmanTestDecision, list[BayesianFriedmanTestDecision]]
Global decision and decisions for each permutation.
"""
n_algorithms, _ = rank_matrix.shape
decisions: list[BayesianFriedmanTestDecision] = []
for perm in permutations(range(1, n_algorithms + 1)):
prior_mean_ranks = np.array(perm, dtype=float)
posterior_mean_ranks = self._compute_posterior_mean_ranks(
rank_matrix, prior_mean_ranks
)
posterior_covariance = self._compute_posterior_covariance(
rank_matrix, prior_mean_ranks, posterior_mean_ranks
)
h = self._evaluate_hypothesis(
prior_mean_ranks, posterior_mean_ranks, posterior_covariance, threshold
)
decisions.append(h)
if len(set(decisions)) > 1:
global_decision = BayesianFriedmanTestDecision.INDETERMINATE
else:
global_decision = decisions[0]
return global_decision, decisions
def _validate_scores(self, scores: EvaluationScores):
"""Validate that the input scores meet the test assumptions.
Parameters
----------
scores : EvaluationScores
Evaluation scores container compatible with this test.
Raises
------
ValueError
If score structure (e.g. number of datasets/algorithms) is not
compatible with the test.
"""
if scores.n_datasets < 1:
raise ValueError("Bayesian Friedman test requires at least one dataset")
if scores.n_algorithms < 2:
raise ValueError("Bayesian Friedman test requires at least two algorithms")
if scores.n_algorithms >= 10:
raise ValueError(
"Bayesian Friedman test requires less than 10 algorithms to be run in "
"imprecise mode due to computational complexity. Consider using precise "
"mode or fewer algorithms."
)
if scores.n_datasets <= scores.n_algorithms:
raise ValueError(
"Bayesian Friedman test requires n_datasets > n_algorithms "
"to obtain valid F distribution degrees of freedom"
)
def _get_metadata(self) -> dict[str, Any]:
"""Return procedure specific metadata to include in the result.
Returns
-------
dict
Metadata mapping.
"""
return {
"prior_strength": self.prior_strength,
"alpha": self.alpha,
"imprecise": self.imprecise,
}
def __validate_arguments(
self,
prior_strength: float,
alpha: float,
imprecise: bool,
):
"""Validate initialization arguments.
Raises
------
ValueError
If any argument is invalid.
"""
if not isinstance(prior_strength, (float, int)) or prior_strength <= 0:
raise ValueError(
f"Prior strength must be a positive number, got {prior_strength}"
)
if not isinstance(alpha, float) or alpha <= 0 or alpha >= 1:
raise ValueError(f"Credibility level must be between 0 and 1, got {alpha}")
if not isinstance(imprecise, bool):
raise ValueError(f"Imprecise must be a boolean, got {imprecise}")