From ab79384a2c8c6daf73eb8483ddda60798fdb5f22 Mon Sep 17 00:00:00 2001 From: Kkkakania <200867803+Kkkakania@users.noreply.github.com> Date: Fri, 24 Jul 2026 11:03:34 +0800 Subject: [PATCH] Add Cramer's V contingency metric --- docs/source/api.rst | 1 + .../api/xskillscore.Contingency.cramers_v.rst | 6 ++ xskillscore/core/contingency.py | 58 ++++++++++++++++ xskillscore/tests/test_contingency.py | 66 +++++++++++++++++++ 4 files changed, 131 insertions(+) create mode 100644 docs/source/api/xskillscore.Contingency.cramers_v.rst diff --git a/docs/source/api.rst b/docs/source/api.rst index 9920673e..a0e46944 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -111,6 +111,7 @@ Multi-Category Metrics :toctree: api/ Contingency.accuracy + Contingency.cramers_v Contingency.gerrity_score Contingency.heidke_score Contingency.peirce_score diff --git a/docs/source/api/xskillscore.Contingency.cramers_v.rst b/docs/source/api/xskillscore.Contingency.cramers_v.rst new file mode 100644 index 00000000..d4e0a385 --- /dev/null +++ b/docs/source/api/xskillscore.Contingency.cramers_v.rst @@ -0,0 +1,6 @@ +xskillscore.Contingency.cramers\_v +================================== + +.. currentmodule:: xskillscore + +.. automethod:: Contingency.cramers_v diff --git a/xskillscore/core/contingency.py b/xskillscore/core/contingency.py index 21777c0f..f02d3f20 100644 --- a/xskillscore/core/contingency.py +++ b/xskillscore/core/contingency.py @@ -641,6 +641,64 @@ def accuracy(self) -> XArray: return corr / N + def cramers_v(self) -> XArray: + """Return Cramér's V for a contingency table. + + Cramér's V measures the association between categorical observations + and forecasts. It ranges from 0 (no association) to 1 (perfect + association). + + .. math:: + V = \\sqrt{\\frac{\\chi^2}{N \\min(r - 1, c - 1)}} + + where :math:`N` is the number of samples and :math:`r` and :math:`c` + are the numbers of non-empty observation and forecast categories, + respectively. + + Returns + ------- + xarray.Dataset or xarray.DataArray + An array containing Cramér's V. + + See Also + -------- + scipy.stats.contingency.association + + References + ---------- + https://en.wikipedia.org/wiki/Cram%C3%A9r%27s_V + """ + observation_dim = OBSERVATIONS_NAME + "_category" + forecast_dim = FORECASTS_NAME + "_category" + declared_degrees = min( + self.table.sizes[observation_dim] - 1, + self.table.sizes[forecast_dim] - 1, + ) + if declared_degrees < 1: + raise ValueError("Cramér's V requires at least two observation and forecast categories") + + total = self._sum_categories("total") + observation_totals = self.table.sum(dim=forecast_dim, skipna=True) + forecast_totals = self.table.sum(dim=observation_dim, skipna=True) + expected = observation_totals * forecast_totals / total + valid_expected = expected.where(expected > 0) + contributions = ((self.table - valid_expected) ** 2 / valid_expected).fillna(0) + chi_squared = contributions.sum( + dim=(observation_dim, forecast_dim), + skipna=True, + ) + degrees = ( + xr.apply_ufunc( + np.minimum, + (observation_totals > 0).sum(dim=observation_dim), + (forecast_totals > 0).sum(dim=forecast_dim), + dask="allowed", + ) + - 1 + ) + + return np.sqrt(chi_squared / (total * degrees)).where(degrees > 0) + def heidke_score(self) -> XArray: """Returns the Heidke skill score(s) for a contingency table with K categories diff --git a/xskillscore/tests/test_contingency.py b/xskillscore/tests/test_contingency.py index 5e1cf1fd..d5138070 100644 --- a/xskillscore/tests/test_contingency.py +++ b/xskillscore/tests/test_contingency.py @@ -1,6 +1,8 @@ import numpy as np import numpy.testing as npt import pytest +import xarray as xr +from scipy.stats.contingency import association from sklearn.metrics import confusion_matrix from xskillscore import Contingency @@ -91,3 +93,67 @@ def test_dichotomous_scores(dichotomous_Contingency_1d, method, expected): """ xs_score = getattr(dichotomous_Contingency_1d, method)().item() npt.assert_almost_equal(xs_score, expected) + + +def test_cramers_v_matches_scipy(dichotomous_Contingency_1d): + """Test Cramer's V against scipy.stats.contingency.association.""" + table = dichotomous_Contingency_1d.table.transpose( + "observations_category", "forecasts_category" + ) + + actual = dichotomous_Contingency_1d.cramers_v().item() + + npt.assert_allclose(actual, association(table.values, method="cramer")) + + +def test_cramers_v_dataset_and_dask(observation_3d_int, forecast_3d_int): + """Test that Cramer's V preserves non-reduced dimensions and laziness.""" + category_edges = np.array([-np.inf, 2.5, 5.5, np.inf]) + contingency = Contingency( + observation_3d_int.to_dataset(name="var").chunk(), + forecast_3d_int.to_dataset(name="var").chunk(), + category_edges, + category_edges, + dim="time", + ) + + result = contingency.cramers_v() + + assert result["var"].dims == ("lat", "lon") + assert result["var"].chunks is not None + assert ((result >= 0) & (result <= 1)).to_array().all().compute().item() + + +def test_cramers_v_ignores_empty_categories(): + """Test that declared but empty categories do not change Cramér's V.""" + observations = xr.DataArray([0, 0, 0, 1, 1, 1], dims="x") + forecasts = xr.DataArray([0, 0, 1, 0, 1, 1], dims="x") + category_edges = np.array([-0.5, 0.5, 1.5, 2.5]) + contingency = Contingency( + observations, + forecasts, + category_edges, + category_edges, + dim="x", + ) + + actual = contingency.cramers_v().item() + expected = association(np.array([[2, 1], [1, 2]]), method="cramer") + + npt.assert_allclose(actual, expected) + + +def test_cramers_v_requires_two_categories(): + """Test that Cramer's V rejects a table with only one category.""" + values = np.arange(4) + category_edges = np.array([-0.5, 3.5]) + contingency = Contingency( + xr.DataArray(values, dims="x"), + xr.DataArray(values, dims="x"), + category_edges, + category_edges, + dim="x", + ) + + with pytest.raises(ValueError, match="at least two observation and forecast categories"): + contingency.cramers_v()