from typing import Optional, Tuple, Union, cast
import numpy as np
import scipy
from sklearn.utils import check_random_state
from sklearn.utils.validation import check_array, column_or_1d
from ._machine_precision import EPSILON
from ._typing import ArrayLike, NDArray
from .utils import (calc_bins, check_alpha, check_array_inf, check_array_nan,
check_array_shape_classification,
check_array_shape_regression, check_arrays_length,
check_binary_zero_one, check_lower_upper_bounds,
check_nb_intervals_sizes, check_nb_sets_sizes,
check_number_bins, check_split_strategy)
[docs]def regression_coverage_score(
y_true: ArrayLike,
y_pred_low: ArrayLike,
y_pred_up: ArrayLike,
) -> float:
"""
Effective coverage score obtained by the prediction intervals.
The effective coverage is obtained by estimating the fraction
of true labels that lie within the prediction intervals.
Parameters
----------
y_true: ArrayLike of shape (n_samples,)
True labels.
y_pred_low: ArrayLike of shape (n_samples,)
Lower bound of prediction intervals.
y_pred_up: ArrayLike of shape (n_samples,)
Upper bound of prediction intervals.
Returns
-------
float
Effective coverage obtained by the prediction intervals.
Examples
--------
>>> from mapie.metrics import regression_coverage_score
>>> import numpy as np
>>> y_true = np.array([5, 7.5, 9.5, 10.5, 12.5])
>>> y_pred_low = np.array([4, 6, 9, 8.5, 10.5])
>>> y_pred_up = np.array([6, 9, 10, 12.5, 12])
>>> print(regression_coverage_score(y_true, y_pred_low, y_pred_up))
0.8
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_pred_low = cast(NDArray, column_or_1d(y_pred_low))
y_pred_up = cast(NDArray, column_or_1d(y_pred_up))
check_arrays_length(y_true, y_pred_low, y_pred_up)
check_lower_upper_bounds(y_true, y_pred_low, y_pred_up)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_pred_low)
check_array_inf(y_pred_low)
check_array_nan(y_pred_up)
check_array_inf(y_pred_up)
coverage = np.mean(
((y_pred_low <= y_true) & (y_pred_up >= y_true))
)
return float(coverage)
[docs]def classification_coverage_score(
y_true: ArrayLike,
y_pred_set: ArrayLike
) -> float:
"""
Effective coverage score obtained by the prediction sets.
The effective coverage is obtained by estimating the fraction
of true labels that lie within the prediction sets.
Parameters
----------
y_true: ArrayLike of shape (n_samples,)
True labels.
y_pred_set: ArrayLike of shape (n_samples, n_class)
Prediction sets given by booleans of labels.
Returns
-------
float
Effective coverage obtained by the prediction sets.
Examples
--------
>>> from mapie.metrics import classification_coverage_score
>>> import numpy as np
>>> y_true = np.array([3, 3, 1, 2, 2])
>>> y_pred_set = np.array([
... [False, False, True, True],
... [False, True, False, True],
... [False, True, True, False],
... [False, False, True, True],
... [False, True, False, True]
... ])
>>> print(classification_coverage_score(y_true, y_pred_set))
0.8
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_pred_set = cast(
NDArray,
check_array(
y_pred_set, force_all_finite=True, dtype=["bool"]
)
)
check_arrays_length(y_true, y_pred_set)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_pred_set)
check_array_inf(y_pred_set)
coverage = np.take_along_axis(
y_pred_set, y_true.reshape(-1, 1), axis=1
).mean()
return float(coverage)
[docs]def regression_mean_width_score(
y_pred_low: ArrayLike,
y_pred_up: ArrayLike
) -> float:
"""
Effective mean width score obtained by the prediction intervals.
Parameters
----------
y_pred_low: ArrayLike of shape (n_samples,)
Lower bound of prediction intervals.
y_pred_up: ArrayLike of shape (n_samples,)
Upper bound of prediction intervals.
Returns
-------
float
Effective mean width of the prediction intervals.
Examples
--------
>>> from mapie.metrics import regression_mean_width_score
>>> import numpy as np
>>> y_pred_low = np.array([4, 6, 9, 8.5, 10.5])
>>> y_pred_up = np.array([6, 9, 10, 12.5, 12])
>>> print(regression_mean_width_score(y_pred_low, y_pred_up))
2.3
"""
y_pred_low = cast(NDArray, column_or_1d(y_pred_low))
y_pred_up = cast(NDArray, column_or_1d(y_pred_up))
check_arrays_length(y_pred_low, y_pred_up)
check_array_nan(y_pred_low)
check_array_inf(y_pred_low)
check_array_nan(y_pred_up)
check_array_inf(y_pred_up)
mean_width = np.abs(y_pred_up - y_pred_low).mean()
return float(mean_width)
[docs]def classification_mean_width_score(y_pred_set: ArrayLike) -> float:
"""
Mean width of prediction set output by
:class:`~mapie.classification.MapieClassifier`.
Parameters
----------
y_pred_set: ArrayLike of shape (n_samples, n_class)
Prediction sets given by booleans of labels.
Returns
-------
float
Mean width of the prediction set.
Examples
--------
>>> from mapie.metrics import classification_mean_width_score
>>> import numpy as np
>>> y_pred_set = np.array([
... [False, False, True, True],
... [False, True, False, True],
... [False, True, True, False],
... [False, False, True, True],
... [False, True, False, True]
... ])
>>> print(classification_mean_width_score(y_pred_set))
2.0
"""
y_pred_set = cast(
NDArray,
check_array(
y_pred_set, force_all_finite=True, dtype=["bool"]
)
)
check_array_nan(y_pred_set)
check_array_inf(y_pred_set)
mean_width = y_pred_set.sum(axis=1).mean()
return float(mean_width)
[docs]def expected_calibration_error(
y_true: ArrayLike,
y_scores: ArrayLike,
num_bins: int = 50,
split_strategy: Optional[str] = None,
) -> float:
"""
The expected calibration error, which is the difference between
the confidence scores and accuracy per bin [1].
[1] Naeini, Mahdi Pakdaman, Gregory Cooper, and Milos Hauskrecht.
"Obtaining well calibrated probabilities using bayesian binning."
Twenty-Ninth AAAI Conference on Artificial Intelligence. 2015.
Parameters
----------
y_true: ArrayLike of shape (n_samples,)
The target values for the calibrator.
y_score: ArrayLike of shape (n_samples,) or (n_samples, n_classes)
The predictions scores.
num_bins: int
Number of bins to make the split in the y_score. The allowed
values are num_bins above 0.
split_strategy: str
The way of splitting the predictions into different bins.
The allowed split strategies are "uniform", "quantile" and
"array split".
Returns
-------
float
The score of ECE (Expected Calibration Error).
"""
split_strategy = check_split_strategy(split_strategy)
num_bins = check_number_bins(num_bins)
y_true_ = check_binary_zero_one(y_true)
y_scores = cast(NDArray, y_scores)
check_arrays_length(y_true_, y_scores)
check_array_nan(y_true_)
check_array_inf(y_true_)
check_array_nan(y_scores)
check_array_inf(y_scores)
if np.size(y_scores.shape) == 2:
y_score = cast(
NDArray, column_or_1d(np.nanmax(y_scores, axis=1))
)
else:
y_score = cast(NDArray, column_or_1d(y_scores))
_, bin_accs, bin_confs, bin_sizes = calc_bins(
y_true_, y_score, num_bins, split_strategy
)
return np.divide(
np.sum(bin_sizes * np.abs(bin_accs - bin_confs)),
np.sum(bin_sizes)
)
[docs]def top_label_ece(
y_true: ArrayLike,
y_scores: ArrayLike,
y_score_arg: Optional[ArrayLike] = None,
num_bins: int = 50,
split_strategy: Optional[str] = None,
classes: Optional[ArrayLike] = None,
) -> float:
"""
The Top-Label ECE which is a method adapted to fit the
ECE to a Top-Label setting [2].
[2] Gupta, Chirag, and Aaditya K. Ramdas.
"Top-label calibration and multiclass-to-binary reductions."
arXiv preprint arXiv:2107.08353 (2021).
Parameters
----------
y_true: ArrayLike of shape (n_samples,)
The target values for the calibrator.
y_scores: ArrayLike of shape (n_samples, n_classes)
or (n_samples,)
The predictions scores, either the maximum score and the
argmax needs to be inputted or in the form of the prediction
probabilities.
y_score_arg: Optional[ArrayLike] of shape (n_samples,)
If only the maximum is provided in the y_scores, the argmax must
be provided here. This is optional and could be directly infered
from the y_scores.
num_bins: int
Number of bins to make the split in the y_score. The allowed
values are num_bins above 0.
split_strategy: str
The way of splitting the predictions into different bins.
The allowed split strategies are "uniform", "quantile" and
"array split".
classes: ArrayLike of shape (n_samples,)
The different classes, in order of the indices that would be
present in a pred_proba.
Returns
-------
float
The ECE score adapted in the top label setting.
"""
y_scores = cast(NDArray, y_scores)
y_true = cast(NDArray, y_true)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_scores)
check_array_inf(y_scores)
if y_score_arg is None:
check_arrays_length(y_true, y_scores)
else:
y_score_arg = cast(NDArray, y_score_arg)
check_array_nan(y_score_arg)
check_array_inf(y_score_arg)
check_arrays_length(y_true, y_scores, y_score_arg)
ece = float(0.)
split_strategy = check_split_strategy(split_strategy)
num_bins = check_number_bins(num_bins)
y_true = cast(NDArray, column_or_1d(y_true))
if y_score_arg is None:
y_score = cast(
NDArray, column_or_1d(np.nanmax(y_scores, axis=1))
)
if classes is None:
y_score_arg = cast(
NDArray, column_or_1d(np.nanargmax(y_scores, axis=1))
)
else:
classes = cast(NDArray, classes)
y_score_arg = cast(
NDArray, column_or_1d(classes[np.nanargmax(y_scores, axis=1)])
)
else:
y_score = cast(NDArray, column_or_1d(y_scores))
y_score_arg = cast(NDArray, column_or_1d(y_score_arg))
labels = np.unique(y_score_arg)
for label in labels:
label_ind = np.where(label == y_score_arg)[0]
y_true_ = np.array(y_true[label_ind] == label, dtype=int)
ece += expected_calibration_error(
y_true_,
y_scores=y_score[label_ind],
num_bins=num_bins,
split_strategy=split_strategy
)
ece /= len(labels)
return ece
[docs]def regression_coverage_score_v2(
y_true: NDArray,
y_intervals: NDArray,
) -> NDArray:
"""
Effective coverage score obtained by the prediction intervals.
The effective coverage is obtained by estimating the fraction
of true labels that lie within the prediction intervals.
It is different from ``regression_coverage_score`` because it uses
directly the output of ``predict`` method and can compute the
coverage for each alpha.
Parameters
----------
y_true: NDArray of shape (n_samples, n_alpha) or (n_samples,)
True labels.
y_intervals: NDArray of shape (n_samples, 2, n_alpha)
Lower and upper bound of prediction intervals
with different alpha risks.
Returns
-------
NDArray of shape (n_alpha,)
Effective coverage obtained by the prediction intervals.
"""
check_arrays_length(y_true, y_intervals)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_intervals)
check_array_inf(y_intervals)
y_intervals = check_array_shape_regression(y_true, y_intervals)
if len(y_true.shape) != 2:
y_true = cast(NDArray, column_or_1d(y_true))
y_true = np.expand_dims(y_true, axis=1)
coverages = np.mean(
np.logical_and(
np.less_equal(y_intervals[:, 0, :], y_true),
np.greater_equal(y_intervals[:, 1, :], y_true)
),
axis=0
)
return coverages
[docs]def classification_coverage_score_v2(
y_true: NDArray,
y_pred_set: NDArray
) -> NDArray:
"""
Effective coverage score obtained by the prediction sets.
The effective coverage is obtained by estimating the fraction
of true labels that lie within the prediction sets.
It is different from ``classification_coverage_score`` because it uses
directly the output of ``predict`` method and can compute the
coverage for each alpha.
Parameters
----------
y_true: NDArray of shape (n_samples, n_alpha) or (n_samples,)
True labels.
y_pred_set: NDArray of shape (n_samples, n_class, n_alpha)
Prediction sets given by booleans of labels.
Returns
-------
NDArray of shape (n_alpha,)
Effective coverage obtained by the prediction sets.
"""
check_arrays_length(y_true, y_pred_set)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_pred_set)
check_array_inf(y_pred_set)
y_pred_set = check_array_shape_classification(y_true, y_pred_set)
if len(y_true.shape) != 2:
y_true = cast(NDArray, column_or_1d(y_true))
y_true = np.expand_dims(y_true, axis=1)
y_true = np.expand_dims(y_true, axis=1)
coverage = np.nanmean(
np.take_along_axis(y_pred_set, y_true, axis=1),
axis=0
)
return coverage[0]
[docs]def regression_ssc(
y_true: NDArray,
y_intervals: NDArray,
num_bins: int = 3
) -> NDArray:
"""
Compute Size-Stratified Coverage metrics proposed in [3] that is
the conditional coverage conditioned by the size of the intervals.
The intervals are ranked by their size (ascending) and then divided into
num_bins groups: one value of coverage by groups is computed.
Warning: This metric should be used only with non constant intervals
(intervals of different sizes), with constant intervals the result
may be misinterpreted.
[3] Angelopoulos, A. N., & Bates, S. (2021).
A gentle introduction to conformal prediction and
distribution-free uncertainty quantification.
arXiv preprint arXiv:2107.07511.
Parameters
----------
y_true: NDArray of shape (n_samples,)
True labels.
y_intervals: NDArray of shape (n_samples, 2, n_alpha) or (n_samples, 2)
Prediction intervals given by booleans of labels.
num_bins: int n
Number of groups. Should be less than the number of different
interval widths.
Returns
-------
NDArray of shape (n_alpha, num_bins)
Examples
--------
>>> from mapie.metrics import regression_ssc
>>> import numpy as np
>>> y_true = np.array([5, 7.5, 9.5])
>>> y_intervals = np.array([
... [4, 6],
... [6.0, 9.0],
... [9, 10.0]
... ])
>>> print(regression_ssc(y_true, y_intervals, num_bins=2))
[[1. 1.]]
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_intervals = check_array_shape_regression(y_true, y_intervals)
check_number_bins(num_bins)
widths = np.abs(y_intervals[:, 1, :] - y_intervals[:, 0, :])
check_nb_intervals_sizes(widths, num_bins)
check_arrays_length(y_true, y_intervals)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_intervals)
check_array_inf(y_intervals)
indexes_sorted = np.argsort(widths, axis=0)
indexes_bybins = np.array_split(indexes_sorted, num_bins, axis=0)
coverages = np.zeros((y_intervals.shape[2], num_bins))
for i, indexes in enumerate(indexes_bybins):
intervals_binned = np.stack([
np.take_along_axis(y_intervals[:, 0, :], indexes, axis=0),
np.take_along_axis(y_intervals[:, 1, :], indexes, axis=0)
], axis=1)
coverages[:, i] = regression_coverage_score_v2(y_true[indexes],
intervals_binned)
return coverages
[docs]def regression_ssc_score(
y_true: NDArray,
y_intervals: NDArray,
num_bins: int = 3
) -> NDArray:
"""
Aggregate by the minimum for each alpha the Size-Stratified Coverage [3]:
returns the maximum violation of the conditional coverage
(with the groups defined).
Warning: This metric should be used only with non constant intervals
(intervals of different sizes), with constant intervals the result
may be misinterpreted.
Parameters
----------
y_true: NDArray of shape (n_samples,)
True labels.
y_intervals: NDArray of shape (n_samples, 2, n_alpha) or (n_samples, 2)
Prediction intervals given by booleans of labels.
num_bins: int n
Number of groups. Should be less than the number of different
interval widths.
Returns
-------
NDArray of shape (n_alpha,)
Examples
--------
>>> from mapie.metrics import regression_ssc
>>> import numpy as np
>>> y_true = np.array([5, 7.5, 9.5])
>>> y_intervals = np.array([
... [[4, 4], [6, 7.5]],
... [[6.0, 8], [9.0, 10]],
... [[9, 9], [10.0, 10.0]]
... ])
>>> print(regression_ssc_score(y_true, y_intervals, num_bins=2))
[1. 0.5]
"""
return np.min(regression_ssc(y_true, y_intervals, num_bins), axis=1)
[docs]def classification_ssc(
y_true: NDArray,
y_pred_set: NDArray,
num_bins: Union[int, None] = None
) -> NDArray:
"""
Compute Size-Stratified Coverage metrics proposed in [3] that is
the conditional coverage conditioned by the size of the predictions sets.
The sets are ranked by their size (ascending) and then divided into
num_bins groups: one value of coverage by groups is computed.
[3] Angelopoulos, A. N., & Bates, S. (2021).
A gentle introduction to conformal prediction and
distribution-free uncertainty quantification.
arXiv preprint arXiv:2107.07511.
Parameters
----------
y_true: NDArray of shape (n_samples,)
True labels.
y_pred_set: NDArray of shape (n_samples, n_class, n_alpha)
or (n_samples, n_class)
Prediction sets given by booleans of labels.
num_bins: int or None
Number of groups. If None, one value of coverage by possible
size of sets (n_classes +1) is computed. Should be less than the
number of different set sizes.
Returns
-------
NDArray of shape (n_alpha, num_bins)
Examples
--------
>>> from mapie.metrics import classification_ssc
>>> import numpy as np
>>> y_true = y_true_class = np.array([3, 3, 1, 2, 2])
>>> y_pred_set = np.array([
... [True, True, True, True],
... [False, True, False, True],
... [True, True, True, False],
... [False, False, True, True],
... [True, True, False, True]])
>>> print(classification_ssc(y_true, y_pred_set, num_bins=2))
[[1. 0.66666667]]
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_pred_set = check_array_shape_classification(y_true, y_pred_set)
check_arrays_length(y_true, y_pred_set)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_pred_set)
check_array_inf(y_pred_set)
sizes = np.sum(y_pred_set, axis=1)
n_classes = y_pred_set.shape[1]
if num_bins is None:
bins = list(range(n_classes + 1))
else:
check_nb_sets_sizes(sizes, num_bins)
check_number_bins(num_bins)
bins = [
b[0] for b in np.array_split(range(n_classes + 1), num_bins)
]
digitized_sizes = np.digitize(sizes, bins)
coverages = np.zeros((y_pred_set.shape[2], len(bins)))
for alpha in range(y_pred_set.shape[2]):
indexes_bybins = [
np.argwhere(digitized_sizes[:, alpha] == i)
for i in range(1, len(bins)+1)
]
for i, indexes in enumerate(indexes_bybins):
coverages[alpha, i] = classification_coverage_score_v2(
y_true[indexes],
np.take_along_axis(
y_pred_set[:, :, alpha],
indexes,
axis=0
)
)
return coverages
[docs]def classification_ssc_score(
y_true: NDArray,
y_pred_set: NDArray,
num_bins: Union[int, None] = None
) -> NDArray:
"""
Aggregate by the minimum for each alpha the Size-Stratified Coverage [3]:
returns the maximum violation of the conditional coverage
(with the groups defined).
Parameters
----------
y_true: NDArray of shape (n_samples,)
True labels.
y_pred_set: NDArray of shape (n_samples, n_class, n_alpha)
or (n_samples, n_class)
Prediction sets given by booleans of labels.
num_bins: int or None
Number of groups. If None, one value of coverage by possible
size of sets (n_classes +1) is computed. Should be less than
the number of different set sizes.
Returns
-------
NDArray of shape (n_alpha,)
Examples
--------
>>> from mapie.metrics import classification_ssc_score
>>> import numpy as np
>>> y_true = y_true_class = np.array([3, 3, 1, 2, 2])
>>> y_pred_set = np.array([
... [True, True, True, True],
... [False, True, False, True],
... [True, True, True, False],
... [False, False, True, True],
... [True, True, False, True]])
>>> print(classification_ssc_score(y_true, y_pred_set, num_bins=2))
[0.66666667]
"""
check_arrays_length(y_true, y_pred_set)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_pred_set)
check_array_inf(y_pred_set)
return np.nanmin(classification_ssc(y_true, y_pred_set, num_bins), axis=1)
def _gaussian_kernel(
x: NDArray,
kernel_size: int
) -> NDArray:
"""
Computes the gaussian kernel of x. (Used in hsic function)
Parameters
----------
x: NDArray
The values from which to compute the gaussian kernel.
kernel_size: int
The variance (sigma), this coefficient controls the width of the curve.
"""
norm_x = x ** 2
dist = -2 * np.matmul(x, x.transpose((0, 2, 1))) \
+ norm_x + norm_x.transpose((0, 2, 1))
return np.exp(-dist / kernel_size)
[docs]def hsic(
y_true: NDArray,
y_intervals: NDArray,
kernel_sizes: ArrayLike = (1, 1)
) -> NDArray:
"""
Compute the square root of the hsic coefficient. HSIC is Hilbert-Schmidt
independence criterion that is a correlation measure. Here we use it as
proposed in [4], to compute the correlation between the indicator of
coverage and the interval size.
If hsic is 0, the two variables (the indicator of coverage and the
interval size) are independant.
Warning: This metric should be used only with non constant intervals
(intervals of different sizes), with constant intervals the result
may be misinterpreted.
[4] Feldman, S., Bates, S., & Romano, Y. (2021).
Improving conditional coverage via orthogonal quantile regression.
Advances in Neural Information Processing Systems, 34, 2060-2071.
Parameters
----------
y_true: NDArray of shape (n_samples,)
True labels.
y_intervals: NDArray of shape (n_samples, 2, n_alpha) or (n_samples, 2)
Prediction sets given by booleans of labels.
kernel_sizes: ArrayLike of size (2,)
The variance (sigma) for each variable (the indicator of coverage and
the interval size), this coefficient controls the width of the curve.
Returns
-------
NDArray of shape (n_alpha,)
One hsic correlation coefficient by alpha.
Raises
------
ValueError
If kernel_sizes has a length different from 2
and if it has negative or null values.
Examples
--------
>>> from mapie.metrics import hsic
>>> import numpy as np
>>> y_true = np.array([9.5, 10.5, 12.5])
>>> y_intervals = np.array([
... [[9, 9], [10.0, 10.0]],
... [[8.5, 9], [12.5, 12]],
... [[10.5, 10.5], [12.0, 12]]
... ])
>>> print(hsic(y_true, y_intervals))
[0.31787614 0.2962914 ]
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_intervals = check_array_shape_regression(y_true, y_intervals)
check_arrays_length(y_true, y_intervals)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_intervals)
check_array_inf(y_intervals)
kernel_sizes = cast(NDArray, column_or_1d(kernel_sizes))
if len(kernel_sizes) != 2:
raise ValueError(
"kernel_sizes should be an ArrayLike of length 2"
)
if (kernel_sizes <= 0).any():
raise ValueError(
"kernel_size should be positive"
)
n_samples, _, n_alpha = y_intervals.shape
y_true_per_alpha = np.tile(y_true, (n_alpha, 1)).transpose()
widths = np.expand_dims(
np.abs(y_intervals[:, 1, :] - y_intervals[:, 0, :]).transpose(),
axis=2
)
cov_ind = np.expand_dims(
np.int_(
((y_intervals[:, 0, :] <= y_true_per_alpha) &
(y_intervals[:, 1, :] >= y_true_per_alpha))
).transpose(),
axis=2
)
k_mat = _gaussian_kernel(widths, kernel_sizes[0])
l_mat = _gaussian_kernel(cov_ind, kernel_sizes[1])
h_mat = np.eye(n_samples) - 1 / n_samples * np.ones((n_samples, n_samples))
hsic_mat = np.matmul(l_mat, np.matmul(h_mat, np.matmul(k_mat, h_mat)))
hsic_mat /= ((n_samples - 1) ** 2)
coef_hsic = np.sqrt(np.matrix.trace(hsic_mat, axis1=1, axis2=2))
return coef_hsic
def coverage_width_based(
y_true: ArrayLike,
y_pred_low: ArrayLike,
y_pred_up: ArrayLike,
eta: float,
alpha: float
) -> float:
"""
Coverage Width-based Criterion (CWC) obtained by the prediction intervals.
The effective coverage score is a criterion used to evaluate the quality
of prediction intervals (PIs) based on their coverage and width.
Khosravi, Abbas, Saeid Nahavandi, and Doug Creighton.
"Construction of optimal prediction intervals for load forecasting
problems."
IEEE Transactions on Power Systems 25.3 (2010): 1496-1503.
Parameters
----------
Coverage score : float
Prediction interval coverage probability (Coverage score), which is
the estimated fraction of true labels that lie within the prediction
intervals.
Mean Width Score : float
Prediction interval normalized average width (Mean Width Score),
calculated as the average width of the prediction intervals.
eta : int
A user-defined parameter that balances the contributions of
Mean Width Score and Coverage score in the CWC calculation.
alpha : float
A user-defined parameter representing the designed confidence level of
the PI.
Returns
-------
float
Effective coverage score (CWC) obtained by the prediction intervals.
Notes
-----
The effective coverage score (CWC) is calculated using the following
formula:
CWC = (1 - Mean Width Score) * exp(-eta * (Coverage score - (1-alpha))**2)
The CWC penalizes under- and overcoverage in the same way and summarizes
the quality of the prediction intervals in a single value.
High Eta (Large Positive Value):
When eta is a high positive value, it will strongly
emphasize the contribution of (1-Mean Width Score). This means that the
algorithm will prioritize reducing the average width of the prediction
intervals (Mean Width Score) over achieving a high coverage probability
(Coverage score). The exponential term np.exp(-eta*(Coverage score -
(1-alpha))**2) will have a sharp decline as Coverage score deviates
from (1-alpha). So, achieving a high Coverage score becomes less important
compared to minimizing Mean Width Score.
The impact will be narrower prediction intervals on average, which may
result in more precise but less conservative predictions.
Low Eta (Small Positive Value):
When eta is a low positive value, it will still
prioritize reducing the average width of the prediction intervals
(Mean Width Score) but with less emphasis compared to higher
eta values.
The exponential term will be less steep, meaning that deviations of
Coverage score from (1-alpha) will have a moderate impact.
You'll get a balance between prediction precision and coverage, but the
exact balance will depend on the specific value of eta.
Negative Eta (Any Negative Value):
When eta is negative, it will have a different effect on the formula.
Negative values of eta will cause the exponential term
np.exp(-eta*(Coverage score - (1-alpha))**2) to become larger as
Coverage score deviates from (1-alpha). This means that
a negative eta prioritizes achieving a high coverage probability
(Coverage score) over minimizing Mean Width Score.
In this case, the algorithm will aim to produce wider prediction intervals
to ensure a higher likelihood of capturing the true values within those
intervals, even if it sacrifices precision.
Negative eta values might be used in scenarios where avoiding errors or
outliers is critical.
Null Eta (Eta = 0):
Specifically, when eta is zero, the CWC score becomes equal to
(1 - Mean Width Score), which is equivalent to
(1 - average width of the prediction intervals).
Therefore, in this case, the CWC score is primarily based on the size of
the prediction interval.
Examples
--------
>>> y_true = np.array([5, 7.5, 9.5, 10.5, 12.5])
>>> y_preds_low = np.array([4, 6, 9, 8.5, 10.5])
>>> y_preds_up = np.array([6, 9, 10, 12.5, 12])
>>> eta = 0.01
>>> alpha = 0.1
>>> cwb = coverage_width_based(y_true, y_preds_low, y_preds_up, eta, alpha)
>>> print(np.round(cwb ,2))
0.69
"""
y_true = cast(NDArray, column_or_1d(y_true))
y_pred_low = cast(NDArray, column_or_1d(y_pred_low))
y_pred_up = cast(NDArray, column_or_1d(y_pred_up))
check_alpha(1-alpha)
coverage_score = regression_coverage_score(
y_true,
y_pred_low,
y_pred_up
)
mean_width = regression_mean_width_score(
y_pred_low,
y_pred_up
)
ref_length = np.subtract(
float(y_true.max()),
float(y_true.min())
)
avg_length = mean_width / ref_length
cwc = (1-avg_length)*np.exp(-eta*(coverage_score-(1-alpha))**2)
return float(cwc)
def add_jitter(
x: NDArray,
noise_amplitude: float = 1e-8,
random_state: Optional[Union[int, np.random.RandomState]] = None
) -> NDArray:
"""
Add a tiny normal distributed perturbation to an array x.
Parameters
----------
x : NDArray
The array to jitter.
noise_amplitude : float, optional
The tiny relative noise amplitude to add, by default 1e-8.
random_state: Optional[Union[int, RandomState]]
Pseudo random number generator state used for random sampling.
Pass an int for reproducible output across multiple function calls.
Returns
-------
NDArray
The array x jittered.
Examples
--------
>>> import numpy as np
>>> from mapie.metrics import add_jitter
>>> x = np.array([0, 1, 2, 3, 4])
>>> res = add_jitter(x, random_state=1)
>>> res
array([0. , 0.99999999, 1.99999999, 2.99999997, 4.00000003])
"""
n = len(x)
random_state_np = check_random_state(random_state)
noise = noise_amplitude * random_state_np.normal(size=n)
x_jittered = x * (1 + noise)
return x_jittered
def sort_xy_by_y(x: NDArray, y: NDArray) -> Tuple[NDArray, NDArray]:
"""
Sort two arrays x and y according to y values.
Parameters
----------
x : NDArray of size (n_samples,)
The array to sort according to y.
y : NDArray of size (n_samples,)
The array to sort.
Returns
-------
Tuple[NDArray, NDArray]
Both arrays sorted.
Examples
--------
>>> import numpy as np
>>> from mapie.metrics import sort_xy_by_y
>>> x = np.array([1, 2, 3, 4, 5])
>>> y = np.array([5, 4, 3, 1, 2])
>>> x_sorted, y_sorted = sort_xy_by_y(x, y)
>>> print(x_sorted)
[4 5 3 2 1]
>>> print(y_sorted)
[1 2 3 4 5]
"""
x = column_or_1d(x)
y = column_or_1d(y)
sort_index = np.argsort(y)
x_sorted = x[sort_index]
y_sorted = y[sort_index]
return x_sorted, y_sorted
[docs]def cumulative_differences(
y_true: NDArray,
y_score: NDArray,
noise_amplitude: float = 1e-8,
random_state: Optional[Union[int, np.random.RandomState]] = 1
) -> NDArray:
"""
Compute the cumulative difference between y_true and y_score, both ordered
according to y_scores array.
Parameters
----------
y_true : NDArray of size (n_samples,)
An array of ground truths.
y_score : NDArray of size (n_samples,)
An array of scores.
noise_amplitude : float, optional
The tiny relative noise amplitude to add, by default 1e-8.
random_state: Optional[Union[int, RandomState]]
Pseudo random number generator state used for random sampling.
Pass an int for reproducible output across multiple function calls.
Returns
-------
NDArray
The mean cumulative difference between y_true and y_score.
References
----------
Arrieta-Ibarra I, Gujral P, Tannen J, Tygert M, Xu C.
Metrics of calibration for probabilistic predictions.
The Journal of Machine Learning Research.
2022 Jan 1;23(1):15886-940.
Examples
--------
>>> import numpy as np
>>> from mapie.metrics import cumulative_differences
>>> y_true = np.array([1, 0, 0])
>>> y_score = np.array([0.7, 0.3, 0.6])
>>> cum_diff = cumulative_differences(y_true, y_score)
>>> print(len(cum_diff))
3
>>> print(np.max(cum_diff) <= 1)
True
>>> print(np.min(cum_diff) >= -1)
True
>>> cum_diff
array([-0.1, -0.3, -0.2])
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
n = len(y_true)
y_score_jittered = add_jitter(
y_score,
noise_amplitude=noise_amplitude,
random_state=random_state
)
y_true_sorted, y_score_sorted = sort_xy_by_y(y_true, y_score_jittered)
cumulative_differences = np.cumsum(y_true_sorted - y_score_sorted)/n
return cumulative_differences
[docs]def length_scale(s: NDArray) -> float:
"""
Compute the mean square root of the sum of s * (1 - s).
This is basically the standard deviation of the
cumulative differences.
Parameters
----------
s : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
The length_scale array.
References
----------
Arrieta-Ibarra I, Gujral P, Tannen J, Tygert M, Xu C.
Metrics of calibration for probabilistic predictions.
The Journal of Machine Learning Research.
2022 Jan 1;23(1):15886-940.
Examples
--------
>>> import numpy as np
>>> from mapie.metrics import length_scale
>>> s = np.array([0, 0, 0.4, 0.3, 0.8])
>>> res = length_scale(s)
>>> print(np.round(res, 2))
0.16
"""
n = len(s)
length_scale = np.sqrt(np.sum(s * (1 - s)))/n
return length_scale
[docs]def kolmogorov_smirnov_statistic(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Kolmogorov-smirnov's statistic for calibration test.
Also called ECCE-MAD
(Estimated Cumulative Calibration Errors - Maximum Absolute Deviation).
The closer to zero, the better the scores are calibrated.
Indeed, if the scores are perfectly calibrated,
the cumulative differences between ``y_true`` and ``y_score``
should share the same properties of a standard Brownian motion
asymptotically.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores..
Returns
-------
float
Kolmogorov-smirnov's statistic.
References
----------
Arrieta-Ibarra I, Gujral P, Tannen J, Tygert M, Xu C.
Metrics of calibration for probabilistic predictions.
The Journal of Machine Learning Research.
2022 Jan 1;23(1):15886-940.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import kolmogorov_smirnov_statistic
>>> y_true = np.array([0, 1, 0, 1, 0])
>>> y_score = np.array([0.1, 0.9, 0.21, 0.9, 0.5])
>>> print(np.round(kolmogorov_smirnov_statistic(y_true, y_score), 3))
0.978
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
y_true = column_or_1d(y_true)
y_score = column_or_1d(y_score)
cum_diff = cumulative_differences(y_true, y_score)
sigma = length_scale(y_score)
ks_stat = np.max(np.abs(cum_diff)) / sigma
return ks_stat
[docs]def kolmogorov_smirnov_cdf(x: float) -> float:
"""
Compute the Kolmogorov-smirnov cumulative distribution
function (CDF) for the float x.
This is interpreted as the CDF of the maximum absolute value
of the standard Brownian motion over the unit interval [0, 1].
The function is approximated by its power series, truncated so as to hit
machine precision error.
Parameters
----------
x : float
The float x to compute the cumulative distribution function on.
Returns
-------
float
The Kolmogorov-smirnov cumulative distribution function.
References
----------
Tygert M.
Calibration of P-values for calibration and for deviation
of a subpopulation from the full population.
arXiv preprint arXiv:2202.00100.
2022 Jan 31.
D. A. Darling. A. J. F. Siegert.
The First Passage Problem for a Continuous Markov Process.
Ann. Math. Statist. 24 (4) 624 - 639, December,
1953.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import kolmogorov_smirnov_cdf
>>> print(np.round(kolmogorov_smirnov_cdf(1), 4))
0.3708
"""
kmax = np.ceil(
0.5 + x * np.sqrt(2) / np.pi * np.sqrt(np.log(4 / (np.pi*EPSILON)))
)
c = 0.0
for k in range(int(kmax)):
kplus = k + 1 / 2
c += (-1)**k / kplus * np.exp(-kplus**2 * np.pi**2 / (2 * x**2))
c *= 2 / np.pi
return c
[docs]def kolmogorov_smirnov_p_value(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Kolmogorov Smirnov p-value.
Deduced from the corresponding statistic and CDF.
It represents the probability of the observed statistic
under the null hypothesis of perfect calibration.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
The Kolmogorov Smirnov p-value.
References
----------
Tygert M.
Calibration of P-values for calibration and for deviation
of a subpopulation from the full population.
arXiv preprint arXiv:2202.00100.
2022 Jan 31.
D. A. Darling. A. J. F. Siegert.
The First Passage Problem for a Continuous Markov Process.
Ann. Math. Statist. 24 (4) 624 - 639, December,
1953.
Example
-------
>>> import pandas as pd
>>> from mapie.metrics import kolmogorov_smirnov_p_value
>>> y_true = np.array([1, 0, 1, 0, 1, 0])
>>> y_score = np.array([0.8, 0.3, 0.5, 0.5, 0.7, 0.1])
>>> ks_p_value = kolmogorov_smirnov_p_value(y_true, y_score)
>>> print(np.round(ks_p_value, 4))
0.7857
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
ks_stat = kolmogorov_smirnov_statistic(y_true, y_score)
ks_p_value = 1 - kolmogorov_smirnov_cdf(ks_stat)
return ks_p_value
[docs]def kuiper_statistic(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Kuiper's statistic for calibration test.
Also called ECCE-R (Estimated Cumulative Calibration Errors - Range).
The closer to zero, the better the scores are calibrated.
Indeed, if the scores are perfectly calibrated,
the cumulative differences between ``y_true`` and ``y_score``
should share the same properties of a standard Brownian motion
asymptotically.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
Kuiper's statistic.
References
----------
Arrieta-Ibarra I, Gujral P, Tannen J, Tygert M, Xu C.
Metrics of calibration for probabilistic predictions.
The Journal of Machine Learning Research.
2022 Jan 1;23(1):15886-940.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import kuiper_statistic
>>> y_true = np.array([0, 1, 0, 1, 0])
>>> y_score = np.array([0.1, 0.9, 0.21, 0.9, 0.5])
>>> print(np.round(kuiper_statistic(y_true, y_score), 3))
0.857
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
y_true = column_or_1d(y_true)
y_score = column_or_1d(y_score)
cum_diff = cumulative_differences(y_true, y_score)
sigma = length_scale(y_score)
ku_stat = (np.max(cum_diff) - np.min(cum_diff)) / sigma
return ku_stat
[docs]def kuiper_cdf(x: float) -> float:
"""
Compute the Kuiper cumulative distribution function (CDF) for the float x.
This is interpreted as the CDF of the range
of the standard Brownian motion over the unit interval [0, 1].
The function is approximated by its power series, truncated so as to hit
machine precision error.
Parameters
----------
x : float
The float x to compute the cumulative distribution function.
Returns
-------
float
The Kuiper cumulative distribution function.
References
----------
Tygert M.
Calibration of P-values for calibration and for deviation
of a subpopulation from the full population.
arXiv preprint arXiv:2202.00100.
2022 Jan 31.
William Feller.
The Asymptotic Distribution of the Range of Sums of
Independent Random Variables.
Ann. Math. Statist. 22 (3) 427 - 432
September, 1951.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import kuiper_cdf
>>> print(np.round(kuiper_cdf(1), 4))
0.0634
"""
kmax = np.ceil(
(
0.5 + x / (np.pi * np.sqrt(2)) *
np.sqrt(
np.log(
4 / (np.sqrt(2 * np.pi) * EPSILON) * (1 / x + x / np.pi**2)
)
)
)
)
c = 0.0
for k in range(int(kmax)):
kplus = k + 1 / 2
c += (
(8 / x**2 + 2 / kplus**2 / np.pi**2) *
np.exp(-2 * kplus**2 * np.pi**2 / x**2)
)
return c
[docs]def kuiper_p_value(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Kuiper statistic p-value.
Deduced from the corresponding statistic and CDF.
It represents the probability of the observed statistic
under the null hypothesis of perfect calibration.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
The Kuiper p-value.
References
----------
Tygert M.
Calibration of P-values for calibration and for deviation
of a subpopulation from the full population.
arXiv preprint arXiv:2202.00100.
2022 Jan 31.
William Feller.
The Asymptotic Distribution of the Range of Sums of
Independent Random Variables.
Ann. Math. Statist. 22 (3) 427 - 432
September, 1951.
Example
-------
>>> import pandas as pd
>>> from mapie.metrics import kuiper_p_value
>>> y_true = np.array([1, 0, 1, 0, 1, 0])
>>> y_score = np.array([0.8, 0.3, 0.5, 0.5, 0.7, 0.1])
>>> ku_p_value = kuiper_p_value(y_true, y_score)
>>> print(np.round(ku_p_value, 4))
0.9684
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
ku_stat = kuiper_statistic(y_true, y_score)
ku_p_value = 1 - kuiper_cdf(ku_stat)
return ku_p_value
[docs]def spiegelhalter_statistic(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Spiegelhalter's statistic for calibration test.
The closer to zero, the better the scores are calibrated.
Indeed, if the scores are perfectly calibrated,
the Brier score simplifies to an expression whose expectancy
and variance are easy to compute. The statistic is no more that
a z-score on this normalized expression.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
Spiegelhalter's statistic.
References
----------
Spiegelhalter DJ.
Probabilistic prediction in patient management and clinical trials.
Statistics in medicine.
1986 Sep;5(5):421-33.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import spiegelhalter_statistic
>>> y_true = np.array([0, 1, 0, 1, 0])
>>> y_score = np.array([0.1, 0.9, 0.21, 0.9, 0.5])
>>> print(np.round(spiegelhalter_statistic(y_true, y_score), 3))
-0.757
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
y_true = column_or_1d(y_true)
y_score = column_or_1d(y_score)
numerator = np.sum(
(y_true - y_score) * (1 - 2 * y_score)
)
denominator = np.sqrt(
np.sum(
(1 - 2 * y_score) ** 2 * y_score * (1 - y_score)
)
)
sp_stat = numerator/denominator
return sp_stat
[docs]def spiegelhalter_p_value(y_true: NDArray, y_score: NDArray) -> float:
"""
Compute Spiegelhalter statistic p-value.
Deduced from the corresponding statistic and CDF,
which is no more than the normal distribution.
It represents the probability of the observed statistic
under the null hypothesis of perfect calibration.
Parameters
----------
y_true : NDArray of shape (n_samples,)
An array of ground truth.
y_score : NDArray of shape (n_samples,)
An array of scores.
Returns
-------
float
The Spiegelhalter statistic p_value.
References
----------
Spiegelhalter DJ.
Probabilistic prediction in patient management and clinical trials.
Statistics in medicine.
1986 Sep;5(5):421-33.
Example
-------
>>> import numpy as np
>>> from mapie.metrics import spiegelhalter_p_value
>>> y_true = np.array([1, 0, 1, 0, 1, 0])
>>> y_score = np.array([0.8, 0.3, 0.5, 0.5, 0.7, 0.1])
>>> sp_p_value = spiegelhalter_p_value(y_true, y_score)
>>> print(np.round(sp_p_value, 4))
0.8486
"""
check_arrays_length(y_true, y_score)
check_array_nan(y_true)
check_array_inf(y_true)
check_array_nan(y_score)
check_array_inf(y_score)
sp_stat = spiegelhalter_statistic(y_true, y_score)
sp_p_value = 1 - scipy.stats.norm.cdf(sp_stat)
return sp_p_value
def regression_mwi_score(
y_true: NDArray,
y_pis: NDArray,
alpha: float
) -> float:
"""
The Winkler score, proposed by Winkler (1972), is a measure used to
evaluate prediction intervals, combining the length of the interval
with a penalty that increases proportionally to the distance of an
observation outside the interval.
Parameters
----------
y_true: ArrayLike of shape (n_samples,)
Ground truth values
y_pis: ArrayLike of shape (n_samples, 2, 1)
Lower and upper bounds of prediction intervals
output from a MAPIE regressor
alpha: float
The value of alpha
Returns
-------
float
The mean Winkler interval score
References
----------
[1] Robert L. Winkler
"A Decision-Theoretic Approach to Interval Estimation",
Journal of the American Statistical Association,
volume 67, pages 187-191 (1972)
(https://doi.org/10.1080/01621459.1972.10481224)
[2] Tilmann Gneiting and Adrian E Raftery
"Strictly Proper Scoring Rules, Prediction, and Estimation",
Journal of the American Statistical Association,
volume 102, pages 359-378 (2007)
(https://doi.org/10.1198/016214506000001437) (Section 6.2)
"""
# Undo any possible quantile crossing
y_pred_low = np.minimum(y_pis[:, 0, 0], y_pis[:, 1, 0])
y_pred_up = np.maximum(y_pis[:, 0, 0], y_pis[:, 1, 0])
check_arrays_length(y_true, y_pred_low, y_pred_up)
# Checking for NaN and inf values
for array in (y_true, y_pred_low, y_pred_up):
check_array_nan(array)
check_array_inf(array)
width = np.sum(y_pred_up) - np.sum(y_pred_low)
error_above = np.sum((y_true - y_pred_up)[y_true > y_pred_up])
error_below = np.sum((y_pred_low - y_true)[y_true < y_pred_low])
total_error = error_above + error_below
mwi = (width + total_error * 2 / alpha) / len(y_true)
return mwi