fairBS and fairRPS

xskillscore/__init__.py

@@ -27,6 +27,7 @@
     crps_gaussian,
     crps_quadrature,
     discrimination,
+    fair_brier_score,
     rank_histogram,
     reliability,
     rps,

xskillscore/core/probabilistic.py

@@ -267,6 +267,63 @@ def brier_score(observations, forecasts, dim=None, weights=None, keep_attrs=Fals
         return res.mean(dim, keep_attrs=keep_attrs)
 
 
+def fair_brier_score(
+    observations,
+    forecasts,
+    dim=None,
+    weights=None,
+    keep_attrs=False,
+    member_dim="member",
+):
+    """Calculate the fair Brier score (FairBS).
+
+    .. math:
+        FairBS(p, k) = change(p_1 - k)^{2}
+
+    Parameters
+    ----------
+    observations : xarray.Dataset or xarray.DataArray
+        The observations or set of observations of the event.
+        Data should be boolean or logical \
+        (True or 1 for event occurance, False or 0 for non-occurance).
+    forecasts : xarray.Dataset or xarray.DataArray
+        The forecast likelihoods of the event. Data should be between 0 and 1.
+    dim : str or list of str, optional
+        Dimension over which to compute mean after computing ``brier_score``.
+        Defaults to None implying averaging over all dimensions.
+    weights : xr.DataArray with dimensions from dim, optional
+        Weights for `weighted.mean(dim)`.
+        Defaults to None, such that no weighting is applied.
+    keep_attrs : bool
+        If True, the attributes (attrs) will be copied
+        from the first input to the new one.
+        If False (default), the new object will
+        be returned without attributes.
+
+    Returns
+    -------
+    xarray.Dataset or xarray.DataArray
+
+    References
+    ----------
+    * https://www-miklip.dkrz.de/about/problems/
+
+    """
+    if member_dim in forecasts.dims:
+        M = forecasts[member_dim].size
+        e = (forecasts == 1).sum(member_dim, keep_attrs=keep_attrs)
+    else:
+        raise ValueError("need forecast with member dim")
+    o = observations
+    with xr.set_options(keep_attrs=keep_attrs):
+        res = e / M - o ** 2 - e * (M - e) / (M ** 2 * (M - 1))
+    res.attrs = observations.attrs  # dirty fix
+    if weights is not None:
+        return res.weighted(weights).mean(dim, keep_attrs=keep_attrs)
+    else:
+        return res.mean(dim, keep_attrs=keep_attrs)
+
+
 def threshold_brier_score(
     observations,
     forecasts,

xskillscore/tests/test_probabilistic.py

@@ -3,6 +3,7 @@
 import properscoring
 import pytest
 import xarray as xr
+from dask import is_dask_collection
 from scipy.stats import norm
 from sklearn.calibration import calibration_curve
 from xarray.tests import assert_allclose, assert_identical
@@ -13,6 +14,7 @@
     crps_gaussian,
     crps_quadrature,
     discrimination,
+    fair_brier_score,
     rank_histogram,
     reliability,
     rps,
@@ -260,6 +262,44 @@ def test_brier_score(o, f_prob, keep_attrs):
         assert actual.attrs == {}
 
 
+@pytest.mark.parametrize("keep_attrs", [True, False])
+def test_fair_brier_score(o, f_prob, keep_attrs):
+    actual = fair_brier_score(
+        (o > 0.5).assign_attrs(**o.attrs),
+        (f_prob > 0.5),
+        keep_attrs=keep_attrs,
+    )
+    assert actual.chunks is None or actual.chunks == ()
+    if keep_attrs:
+        assert actual.attrs == o.attrs
+    else:
+        assert actual.attrs == {}
+
+
+@pytest.mark.parametrize("dim", DIMS)
+def test_fair_brier_score_dim(o, f_prob, dim):
+    actual = fair_brier_score((o > 0.5), (f_prob > 0.5), dim=dim)
+    assert_only_dim_reduced(dim, actual, o)
+
+
+def test_brier_score_vs_fair_brier_score(o, f_prob):
+    dim = ["lon", "lat"]
+    fbs = fair_brier_score((o > 0.5), (f_prob > 0.5), dim=dim)
+    bs = brier_score((o > 0.5), (f_prob > 0.5).mean("member"), dim=dim)
+    print("fairBS", fbs, "\nBS", bs)
+    assert (fbs <= bs).all()
+    assert fbs <= 1
+    assert fbs >= 0
+
+
+def test_fair_brier_score_dask(o_dask, f_prob_dask):
+    actual = brier_score(
+        (o_dask > 0.5).assign_attrs(**o_dask.attrs),
+        (f_prob_dask > 0.5),
+    )
+    assert is_dask_collection(actual)
+
+
 @pytest.mark.parametrize("dim", DIMS)
 def test_threshold_brier_score_dim(o, f_prob, dim):
     actual = threshold_brier_score(o, f_prob, threshold=0.5, dim=dim)