diff --git a/docs/paper_metrics.md b/docs/paper_metrics.md
new file mode 100644
index 0000000..a863625
--- /dev/null
+++ b/docs/paper_metrics.md
@@ -0,0 +1,51 @@
+# Posterior metrics
+
+These reusable metrics score caller-supplied posterior draws and finite loss
+tables. They do not establish calibration, online utility, or a scientific
+workflow.
+
+## Marginal CRPS
+
+`marginal_crps(..., estimator="empirical")` computes the exact CRPS of a
+finite empirical distribution, optionally with nonnegative draw weights.
+Weights are max-rescaled before summation, so finite very-large weights remain
+well-defined; zero-mass rows are rejected.
+
+Contract: `samples` has shape `(cases, draws, targets)` and `truth` has shape
+`(cases, targets)`, with no broadcasting. Optional `weights` has shape
+`(cases, draws)` and is finite, nonnegative, and strictly positive per row.
+The return value has shape `(cases, targets)`. The `"unbiased"` estimator
+requires unweighted iid draws and at least two draws.
+
+`estimator="unbiased"` is the iid finite-sample U-estimator with a single
+fixed truth value. The triangle inequality gives
+`sum_{i != j}|x_i - x_j| <= 2*(M-1)*sum_i|x_i - y|`, so it is nonnegative for
+every finite sample in exact arithmetic; floating-point roundoff can cause
+tiny negative values. Its unbiasedness does not require negative realizations.
+
+## Decisions
+
+`brier_score` accepts finite nonnegative probabilities of shape
+`(cases, classes)` whose rows sum to one, plus integer `labels` of shape
+`(cases,)`. It returns shape `(cases,)`, with the sum of squared errors over
+all classes. For two classes this is exactly factor 2 times the scalar binary
+convention `(p_positive - y)**2`.
+
+`expected_action_losses` accepts a finite loss table of shape
+`(draws, actions)` and optional finite nonnegative `weights` of shape
+`(draws,)` (strictly positive total). It returns shape `(actions,)`, the
+posterior-weighted mean loss per action. `bayes_action_set` has the same input
+contract and returns `(min_risk, minimizers)`, where `minimizers` is a sorted
+integer array containing every exact minimizing action. Their weights use the
+same stable normalization and validation rules.
+
+Tiny executable example:
+
+```python
+import numpy as np
+from skbel.metrics import bayes_action_set, marginal_crps
+
+draws = np.array([[[0.0], [2.0]]])
+print(marginal_crps(draws, np.array([[1.0]])))  # [[0.5]]
+print(bayes_action_set(np.array([[1.0, 3.0], [1.0, 3.0]])))  # (1.0, [0])
+```
diff --git a/skbel/metrics/__init__.py b/skbel/metrics/__init__.py
new file mode 100644
index 0000000..85a1edb
--- /dev/null
+++ b/skbel/metrics/__init__.py
@@ -0,0 +1,20 @@
+#  Copyright (c) 2026. Robin Thibaut, Ghent University
+
+"""General-purpose posterior-evaluation and finite-action decision-risk metrics.
+
+NumPy-only. No TensorFlow/scikit-learn dependency.
+"""
+
+from .posterior import (
+    bayes_action_set,
+    brier_score,
+    expected_action_losses,
+    marginal_crps,
+)
+
+__all__ = [
+    "bayes_action_set",
+    "brier_score",
+    "expected_action_losses",
+    "marginal_crps",
+]
diff --git a/skbel/metrics/posterior.py b/skbel/metrics/posterior.py
new file mode 100644
index 0000000..f1d36d2
--- /dev/null
+++ b/skbel/metrics/posterior.py
@@ -0,0 +1,337 @@
+#  Copyright (c) 2026. Robin Thibaut, Ghent University
+
+"""Posterior-sample scoring and finite-action decision-risk primitives.
+
+Everything here is general-purpose evaluation arithmetic: it consumes arrays
+the caller already has (posterior draws, class probabilities, a loss table)
+and returns per-case scores or risk summaries. Nothing here runs a simulator,
+a prior, or model training, and nothing here supplies a hydrological (or any
+other domain) loss function -- that is an application responsibility.
+
+Conventions enforced throughout, deliberately strict (no silent broadcasting,
+no default clipping of bad inputs, no NaN hiding):
+
+- Shapes are explicit and checked exactly; mismatches raise ``ValueError``.
+- Weights, where accepted, must be finite, non-negative, and sum to a
+  strictly positive value per case (or overall, for 1-D weights); they are
+  then explicitly renormalized by this module (callers do not need to
+  pre-normalize).
+- Outputs are checked for finiteness before being returned.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+__all__ = [
+    "bayes_action_set",
+    "brier_score",
+    "expected_action_losses",
+    "marginal_crps",
+]
+
+
+def _as_float_array(x, name: str) -> np.ndarray:
+    arr = np.asarray(x)
+    if arr.dtype.kind not in "fiu":
+        raise TypeError(f"{name} must be a numeric array, got dtype {arr.dtype}")
+    return arr.astype(np.float64, copy=False)
+
+
+def _check_finite(arr: np.ndarray, name: str) -> None:
+    if not np.all(np.isfinite(arr)):
+        raise ValueError(f"{name} contains non-finite values (NaN or inf)")
+
+
+def _validate_weights_2d(weights, cases: int, draws: int) -> np.ndarray:
+    """Validate and explicitly normalize per-case weights of shape (cases, draws)."""
+    w = _as_float_array(weights, "weights")
+    if w.shape != (cases, draws):
+        raise ValueError(f"weights must have shape {(cases, draws)}, got {w.shape}")
+    _check_finite(w, "weights")
+    if np.any(w < 0):
+        raise ValueError("weights must be non-negative")
+    row_max = np.max(w, axis=1)
+    if np.any(row_max <= 0) or not np.all(np.isfinite(row_max)):
+        raise ValueError("every case's weights must have a strictly positive sum")
+    scaled = w / row_max[:, None]
+    row_sums = scaled.sum(axis=1)
+    if np.any(~np.isfinite(row_sums)) or np.any(row_sums <= 0):
+        raise ValueError("every case's weights must have a finite, strictly positive sum")
+    normalized = scaled / row_sums[:, None]
+    if np.any(~np.isfinite(normalized)) or not np.all(np.sum(normalized, axis=1) > 0):
+        raise ValueError("normalized weights must have finite positive mass")
+    return normalized
+
+
+def _validate_weights_1d(weights, draws: int) -> np.ndarray:
+    """Validate and explicitly normalize weights of shape (draws,)."""
+    w = _as_float_array(weights, "weights")
+    if w.shape != (draws,):
+        raise ValueError(f"weights must have shape {(draws,)}, got {w.shape}")
+    _check_finite(w, "weights")
+    if np.any(w < 0):
+        raise ValueError("weights must be non-negative")
+    scale = np.max(w)
+    if scale <= 0 or not np.isfinite(scale):
+        raise ValueError("weights must have a strictly positive sum")
+    scaled = w / scale
+    total = scaled.sum()
+    if not np.isfinite(total) or total <= 0:
+        raise ValueError("weights must have a finite, strictly positive sum")
+    normalized = scaled / total
+    if np.any(~np.isfinite(normalized)) or normalized.sum() <= 0:
+        raise ValueError("normalized weights must have finite positive mass")
+    return normalized
+
+
+def _uniform_weights_2d(cases: int, draws: int) -> np.ndarray:
+    return np.full((cases, draws), 1.0 / draws, dtype=np.float64)
+
+
+def _weighted_pairwise_abs_diff_sum_sorted(x_1d: np.ndarray, w_1d: np.ndarray) -> float:
+    """Return sum_i sum_j w_i * w_j * |x_i - x_j| via an O(M log M) sort.
+
+    Derivation (sorted ascending x_(1) <= ... <= x_(M), weights normalized to
+    sum to 1, W_{i-1} = cumulative weight strictly before position i):
+
+        S = 2 * sum_i w_i * x_i * (2*W_{i-1} + w_i - 1)
+
+    which follows from expanding S = 2 * sum_{i<j} w_i w_j (x_j - x_i) and
+    using sum_{i<j} w_i = W_{i-1}, sum_{j>i} w_j = 1 - W_i, W_i = W_{i-1} + w_i.
+    """
+    order = np.argsort(x_1d, kind="mergesort")
+    x_sorted = x_1d[order]
+    w_sorted = w_1d[order]
+    cum_w = np.cumsum(w_sorted)
+    w_prev = cum_w - w_sorted  # exclusive cumulative weight before each position
+    # Center to reduce cancellation: CRPS pairwise-difference terms are
+    # translation invariant, so shifting x by a constant before the
+    # multiply-and-sum does not change S but keeps the summed magnitudes small.
+    shift = x_sorted.mean()
+    x_centered = x_sorted - shift
+    terms = w_sorted * x_centered * (2.0 * w_prev + w_sorted - 1.0)
+    return float(2.0 * terms.sum())
+
+
+def _unweighted_pairwise_abs_diff_sum_sorted(x_1d: np.ndarray) -> float:
+    """Return sum_i sum_j |x_i - x_j| (raw, unnormalized) via an O(M log M) sort.
+
+    Closed form for sorted ascending x_(1) <= ... <= x_(M), 1-indexed i:
+
+        sum_{i,j} |x_i - x_j| = 2 * sum_i (2*i - M - 1) * x_(i)
+    """
+    m = x_1d.shape[0]
+    order = np.argsort(x_1d, kind="mergesort")
+    x_sorted = x_1d[order]
+    shift = x_sorted.mean()
+    x_centered = x_sorted - shift
+    i = np.arange(1, m + 1, dtype=np.float64)
+    coeff = 2.0 * i - m - 1.0
+    return float(2.0 * np.sum(coeff * x_centered))
+
+
+def marginal_crps(
+    samples: np.ndarray,
+    truth: np.ndarray,
+    weights: np.ndarray | None = None,
+    estimator: str = "empirical",
+) -> np.ndarray:
+    """Per-case, per-target marginal continuous ranked probability score.
+
+    :param samples: posterior draws, shape ``(cases, draws, targets)``.
+    :param truth: realized/true values, shape ``(cases, targets)``.
+    :param weights: optional per-draw weights, shape ``(cases, draws)``.
+        Must be finite, non-negative, with a strictly positive sum per case;
+        this function explicitly renormalizes them (they need not already
+        sum to 1). If ``None``, draws are weighted uniformly.
+    :param estimator: ``"empirical"`` (default) computes the exact CRPS of
+        the (possibly weighted) empirical distribution against ``truth`` --
+        valid for any non-negative weights and any number of draws
+        ``M >= 1``. ``"unbiased"`` computes the finite-sample U-statistic
+        ("fair") estimator of CRPS against the fixed value ``truth``, under
+        the assumption that draws are unweighted i.i.d. samples from the
+        predictive distribution; it requires ``weights is None`` and
+        ``M >= 2``, and uses the pair denominator ``M*(M-1)`` (no diagonal
+        terms). For a fixed observation ``y``, the triangle inequality gives
+        ``sum_{i != j}|x_i - x_j| <= 2*(M-1)*sum_i|x_i - y|``. Therefore this
+        U-estimator is non-negative in exact arithmetic; floating-point
+        roundoff can produce tiny negative values.
+    :return: array of shape ``(cases, targets)``, the marginal score for
+        each case/target pair (each target scored independently -- this is
+        NOT a joint multivariate score and does not by itself establish
+        joint calibration).
+    """
+    if estimator not in ("empirical", "unbiased"):
+        raise ValueError(f"estimator must be 'empirical' or 'unbiased', got {estimator!r}")
+
+    samples_arr = _as_float_array(samples, "samples")
+    truth_arr = _as_float_array(truth, "truth")
+
+    if samples_arr.ndim != 3:
+        raise ValueError(
+            f"samples must have shape (cases, draws, targets), got ndim={samples_arr.ndim}"
+        )
+    if truth_arr.ndim != 2:
+        raise ValueError(f"truth must have shape (cases, targets), got ndim={truth_arr.ndim}")
+
+    cases, draws, targets = samples_arr.shape
+    if cases == 0 or draws == 0 or targets == 0:
+        raise ValueError(f"samples must have no empty axes, got shape {samples_arr.shape}")
+    if truth_arr.shape != (cases, targets):
+        raise ValueError(
+            f"truth shape {truth_arr.shape} is not compatible with samples shape "
+            f"{samples_arr.shape}: expected exactly {(cases, targets)} -- no broadcasting is performed"
+        )
+
+    _check_finite(samples_arr, "samples")
+    _check_finite(truth_arr, "truth")
+
+    if estimator == "unbiased":
+        if weights is not None:
+            raise ValueError(
+                "estimator='unbiased' requires unweighted i.i.d. draws; got explicit weights. "
+                "There is no supported weighted-U-statistic rule in this function."
+            )
+        if draws < 2:
+            raise ValueError("estimator='unbiased' requires at least 2 draws (M >= 2)")
+
+        result = np.empty((cases, targets), dtype=np.float64)
+        for c in range(cases):
+            for t in range(targets):
+                x = samples_arr[c, :, t]
+                y = truth_arr[c, t]
+                mean_abs_diff_to_y = np.mean(np.abs(x - y))
+                raw_pair_sum = _unweighted_pairwise_abs_diff_sum_sorted(x)
+                result[c, t] = mean_abs_diff_to_y - raw_pair_sum / (2.0 * draws * (draws - 1))
+        _check_finite(result, "marginal_crps result")
+        return result
+
+    # estimator == "empirical"
+    if weights is None:
+        w_norm = _uniform_weights_2d(cases, draws)
+    else:
+        w_norm = _validate_weights_2d(weights, cases, draws)
+
+    result = np.empty((cases, targets), dtype=np.float64)
+    for t in range(targets):
+        x_t = samples_arr[:, :, t]  # (cases, draws)
+        y_t = truth_arr[:, t]  # (cases,)
+        weighted_abs_diff_to_y = np.sum(w_norm * np.abs(x_t - y_t[:, None]), axis=1)  # (cases,)
+        for c in range(cases):
+            pair_sum = _weighted_pairwise_abs_diff_sum_sorted(x_t[c], w_norm[c])
+            result[c, t] = weighted_abs_diff_to_y[c] - 0.5 * pair_sum
+
+    _check_finite(result, "marginal_crps result")
+    return result
+
+
+def brier_score(probabilities: np.ndarray, labels: np.ndarray) -> np.ndarray:
+    """Multiclass Brier score, summed over all classes, returned per case.
+
+    :param probabilities: shape ``(cases, classes)``, finite, non-negative,
+        each row summing to 1 (within ``1e-8`` absolute tolerance) -- not
+        auto-normalized; malformed rows raise.
+    :param labels: integer class labels, shape ``(cases,)``, each in
+        ``[0, classes)``.
+    :return: array of shape ``(cases,)``: ``sum_c (p[c] - 1{c == label})**2``.
+
+        Convention, stated explicitly: this sums the squared error over
+        *every* class, including the label's own class and every other
+        class. For the binary case (``classes == 2``) this reads as exactly
+        twice the single-probability binary Brier-score convention
+        ``(p_positive - y)**2``, because both the positive and negative
+        class terms contribute equally by symmetry. Callers who want the
+        single-probability binary convention should divide this function's
+        binary-case output by 2, rather than assuming this function already
+        returns it.
+    """
+    probs = _as_float_array(probabilities, "probabilities")
+    labels_arr = np.asarray(labels)
+
+    if probs.ndim != 2:
+        raise ValueError(f"probabilities must have shape (cases, classes), got ndim={probs.ndim}")
+    cases, classes = probs.shape
+    if cases == 0 or classes == 0:
+        raise ValueError(f"probabilities must have no empty axes, got shape {probs.shape}")
+
+    if labels_arr.ndim != 1:
+        raise ValueError(f"labels must have shape (cases,), got ndim={labels_arr.ndim}")
+    if labels_arr.shape[0] != cases:
+        raise ValueError(
+            f"labels shape {labels_arr.shape} incompatible with probabilities cases={cases}"
+        )
+    if labels_arr.dtype.kind not in "iu":
+        raise TypeError(f"labels must be an integer array, got dtype {labels_arr.dtype}")
+
+    _check_finite(probs, "probabilities")
+    if np.any(probs < 0):
+        raise ValueError("probabilities must be non-negative")
+    row_sums = probs.sum(axis=1)
+    if not np.allclose(row_sums, 1.0, atol=1e-8, rtol=0.0):
+        raise ValueError("every probability row must sum to 1 (within 1e-8 absolute tolerance)")
+
+    if np.any(labels_arr < 0) or np.any(labels_arr >= classes):
+        raise ValueError(
+            f"labels must lie in [0, {classes}), got range [{labels_arr.min()}, {labels_arr.max()}]"
+        )
+
+    one_hot = np.zeros((cases, classes), dtype=np.float64)
+    one_hot[np.arange(cases), labels_arr] = 1.0
+
+    result = np.sum((probs - one_hot) ** 2, axis=1)
+    _check_finite(result, "brier_score result")
+    return result
+
+
+def expected_action_losses(losses: np.ndarray, weights: np.ndarray | None = None) -> np.ndarray:
+    """Posterior-expected loss per action, from an explicit per-draw loss table.
+
+    This is posterior decision-risk *reduction* over actions the caller has
+    already enumerated and already supplied losses for -- not an online
+    acquisition/design optimizer, and it carries no hydrological (or any
+    other domain) semantics; the loss table's meaning is entirely the
+    caller's responsibility.
+
+    :param losses: shape ``(draws, actions)``, finite.
+    :param weights: optional shape ``(draws,)``, finite, non-negative, with
+        a strictly positive sum; explicitly renormalized by this function.
+        If ``None``, draws are weighted uniformly.
+    :return: array of shape ``(actions,)``, the weighted-mean loss for each
+        action.
+    """
+    losses_arr = _as_float_array(losses, "losses")
+    if losses_arr.ndim != 2:
+        raise ValueError(f"losses must have shape (draws, actions), got ndim={losses_arr.ndim}")
+    draws, actions = losses_arr.shape
+    if draws == 0 or actions == 0:
+        raise ValueError(f"losses must have no empty axes, got shape {losses_arr.shape}")
+    _check_finite(losses_arr, "losses")
+
+    if weights is None:
+        w_norm = np.full(draws, 1.0 / draws, dtype=np.float64)
+    else:
+        w_norm = _validate_weights_1d(weights, draws)
+
+    result = w_norm @ losses_arr
+    _check_finite(result, "expected_action_losses result")
+    return result
+
+
+def bayes_action_set(
+    losses: np.ndarray, weights: np.ndarray | None = None
+) -> tuple[float, np.ndarray]:
+    """Minimum posterior-expected loss and the set of ALL exact minimizers.
+
+    :param losses: shape ``(draws, actions)``, finite.
+    :param weights: see :func:`expected_action_losses`.
+    :return: ``(min_risk, minimizing_action_indices)`` where
+        ``minimizing_action_indices`` is a sorted ``int64`` array containing
+        every action index whose expected loss exactly equals ``min_risk``
+        (a single best action returns a length-1 array; exact ties return
+        every tied index -- there is no arbitrary tie-break).
+    """
+    expected = expected_action_losses(losses, weights)
+    min_risk = float(np.min(expected))
+    minimizers = np.flatnonzero(expected == min_risk).astype(np.int64)
+    return min_risk, minimizers
diff --git a/skbel/testing/test_posterior_metrics.py b/skbel/testing/test_posterior_metrics.py
new file mode 100644
index 0000000..f2b66de
--- /dev/null
+++ b/skbel/testing/test_posterior_metrics.py
@@ -0,0 +1,475 @@
+"""Unit tests for skbel.metrics.posterior.
+
+Import strategy: normal tests require the complete package import. An import
+failure intentionally fails test collection rather than silently bypassing
+``skbel.__init__`` with a direct file import.
+"""
+
+from __future__ import annotations
+
+import unittest
+
+import numpy as np
+
+from skbel.metrics import bayes_action_set, brier_score, expected_action_losses, marginal_crps
+
+# --------------------------------------------------------------------------
+# Independent oracle implementations (naive O(M^2)/O(M) double loops, NOT the
+# production sorted formula) used to cross-check the module under arbitrary
+# random inputs, not just the exact hand-worked tiny cases below.
+# --------------------------------------------------------------------------
+
+
+def _oracle_empirical_crps_1d(x: np.ndarray, w: np.ndarray, y: float) -> float:
+    x = np.asarray(x, dtype=float)
+    w = np.asarray(w, dtype=float)
+    w = w / w.sum()
+    m = len(x)
+    term1 = sum(w[i] * abs(x[i] - y) for i in range(m))
+    term2 = sum(w[i] * w[j] * abs(x[i] - x[j]) for i in range(m) for j in range(m))
+    return term1 - 0.5 * term2
+
+
+def _oracle_unbiased_crps_1d(x: np.ndarray, y: float) -> float:
+    x = np.asarray(x, dtype=float)
+    m = len(x)
+    term1 = sum(abs(xi - y) for xi in x) / m
+    term2 = sum(abs(x[i] - x[j]) for i in range(m) for j in range(m) if i != j)
+    return term1 - term2 / (2.0 * m * (m - 1))
+
+
+def _oracle_brier(probs: np.ndarray, labels: np.ndarray) -> np.ndarray:
+    cases, classes = probs.shape
+    out = np.zeros(cases)
+    for c in range(cases):
+        s = 0.0
+        for k in range(classes):
+            target = 1.0 if k == labels[c] else 0.0
+            s += (probs[c, k] - target) ** 2
+        out[c] = s
+    return out
+
+
+def _oracle_expected_action_losses(losses: np.ndarray, w: np.ndarray) -> np.ndarray:
+    draws, actions = losses.shape
+    w = np.asarray(w, dtype=float)
+    w = w / w.sum()
+    out = np.zeros(actions)
+    for a in range(actions):
+        out[a] = sum(w[i] * losses[i, a] for i in range(draws))
+    return out
+
+
+# --------------------------------------------------------------------------
+# marginal_crps
+# --------------------------------------------------------------------------
+
+
+class TestMarginalCrpsExactTinyCases(unittest.TestCase):
+    """Hand-worked exact cases with analytically determined values."""
+
+    def test_pointmass_pair_empirical(self):
+        # samples = [0, 2], truth = 1 -> empirical CRPS = 0.5
+        samples = np.array([[[0.0], [2.0]]])  # (cases=1, draws=2, targets=1)
+        truth = np.array([[1.0]])
+        result = marginal_crps(samples, truth, estimator="empirical")
+        self.assertEqual(result.shape, (1, 1))
+        self.assertAlmostEqual(result[0, 0], 0.5, places=12)
+
+    def test_pointmass_pair_unbiased_is_zero(self):
+        samples = np.array([[[0.0], [2.0]]])
+        truth = np.array([[1.0]])
+        result = marginal_crps(samples, truth, estimator="unbiased")
+        self.assertAlmostEqual(result[0, 0], 0.0, places=12)
+
+    def test_weighted_pair_masses_quarter_three_quarter(self):
+        # samples = [0, 1], weights = [.25, .75], truth = 0 -> 0.5625
+        samples = np.array([[[0.0], [1.0]]])
+        truth = np.array([[0.0]])
+        weights = np.array([[0.25, 0.75]])
+        result = marginal_crps(samples, truth, weights=weights, estimator="empirical")
+        self.assertAlmostEqual(result[0, 0], 0.5625, places=12)
+
+
+class TestMarginalCrpsOracleCrossCheck(unittest.TestCase):
+    """Cross-check against the independent naive oracle on random data."""
+
+    def test_random_cases_targets_against_oracle_empirical_unweighted(self):
+        rng = np.random.default_rng(20260910)
+        cases, draws, targets = 4, 7, 3
+        samples = rng.normal(size=(cases, draws, targets))
+        truth = rng.normal(size=(cases, targets))
+        result = marginal_crps(samples, truth, estimator="empirical")
+        for c in range(cases):
+            for t in range(targets):
+                expected = _oracle_empirical_crps_1d(
+                    samples[c, :, t], np.full(draws, 1.0 / draws), truth[c, t]
+                )
+                self.assertAlmostEqual(result[c, t], expected, places=9)
+
+    def test_random_cases_targets_against_oracle_empirical_weighted(self):
+        rng = np.random.default_rng(7)
+        cases, draws, targets = 3, 9, 2
+        samples = rng.normal(size=(cases, draws, targets))
+        truth = rng.normal(size=(cases, targets))
+        raw_weights = rng.uniform(0.01, 1.0, size=(cases, draws))
+        result = marginal_crps(samples, truth, weights=raw_weights, estimator="empirical")
+        for c in range(cases):
+            for t in range(targets):
+                expected = _oracle_empirical_crps_1d(samples[c, :, t], raw_weights[c], truth[c, t])
+                self.assertAlmostEqual(result[c, t], expected, places=9)
+
+    def test_random_against_oracle_unbiased(self):
+        rng = np.random.default_rng(99)
+        cases, draws, targets = 3, 12, 2
+        samples = rng.normal(size=(cases, draws, targets))
+        truth = rng.normal(size=(cases, targets))
+        result = marginal_crps(samples, truth, estimator="unbiased")
+        for c in range(cases):
+            for t in range(targets):
+                expected = _oracle_unbiased_crps_1d(samples[c, :, t], truth[c, t])
+                self.assertAlmostEqual(result[c, t], expected, places=9)
+
+    def test_weights_need_not_be_prenormalized(self):
+        # Same relative masses (1:3) scaled up -> identical result to .25/.75
+        samples = np.array([[[0.0], [1.0]]])
+        truth = np.array([[0.0]])
+        weights = np.array([[10.0, 30.0]])  # sums to 40, ratio 1:3
+        result = marginal_crps(samples, truth, weights=weights, estimator="empirical")
+        self.assertAlmostEqual(result[0, 0], 0.5625, places=12)
+
+    def test_huge_weights_are_rescaled_before_sum(self):
+        # Direct summation of these finite weights overflows to inf and used
+        # to produce an all-zero normalized vector and an incorrect CRPS.
+        samples = np.array([[[0.0], [2.0]]])
+        truth = np.array([[1.0]])
+        result = marginal_crps(samples, truth, weights=np.array([[1e308, 1e308]]))
+        self.assertAlmostEqual(result[0, 0], 0.5, places=12)
+
+    def test_huge_and_tiny_weight_rescalings_preserve_relative_mass(self):
+        samples = np.array([[[0.0], [1.0]]])
+        truth = np.array([[0.0]])
+        expected = marginal_crps(samples, truth, weights=np.array([[1.0, 3.0]]))
+        for weights in (np.array([[1e307, 3e307]]), np.array([[1e-320, 3e-320]])):
+            np.testing.assert_allclose(
+                marginal_crps(samples, truth, weights=weights), expected, atol=1e-12
+            )
+
+
+class TestMarginalCrpsShapesAndPermutations(unittest.TestCase):
+    def test_no_silent_broadcasting_truth_shape_mismatch(self):
+        samples = np.zeros((2, 5, 3))
+        bad_truth = np.zeros((2, 4))  # wrong targets count
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, bad_truth)
+
+    def test_samples_must_be_3d(self):
+        with self.assertRaises(ValueError):
+            marginal_crps(np.zeros((2, 5)), np.zeros((2,)))
+
+    def test_multi_case_multi_target_matches_per_case_per_target_loop(self):
+        rng = np.random.default_rng(11)
+        cases, draws, targets = 3, 6, 4
+        samples = rng.normal(size=(cases, draws, targets))
+        truth = rng.normal(size=(cases, targets))
+        batched = marginal_crps(samples, truth, estimator="empirical")
+        for c in range(cases):
+            for t in range(targets):
+                single = marginal_crps(
+                    samples[c : c + 1, :, t : t + 1],
+                    truth[c : c + 1, t : t + 1],
+                    estimator="empirical",
+                )
+                self.assertAlmostEqual(batched[c, t], single[0, 0], places=9)
+
+    def test_permutation_invariance_unweighted(self):
+        rng = np.random.default_rng(3)
+        x = rng.normal(size=8)
+        perm = rng.permutation(8)
+        samples_a = x.reshape(1, 8, 1)
+        samples_b = x[perm].reshape(1, 8, 1)
+        truth = np.array([[0.3]])
+        a = marginal_crps(samples_a, truth, estimator="empirical")
+        b = marginal_crps(samples_b, truth, estimator="empirical")
+        self.assertAlmostEqual(a[0, 0], b[0, 0], places=12)
+
+    def test_permutation_invariance_weighted_moves_with_draws(self):
+        rng = np.random.default_rng(4)
+        x = rng.normal(size=6)
+        w = rng.uniform(0.1, 1.0, size=6)
+        perm = rng.permutation(6)
+        samples_a = x.reshape(1, 6, 1)
+        samples_b = x[perm].reshape(1, 6, 1)
+        w_a = w.reshape(1, 6)
+        w_b = w[perm].reshape(1, 6)
+        truth = np.array([[0.0]])
+        a = marginal_crps(samples_a, truth, weights=w_a, estimator="empirical")
+        b = marginal_crps(samples_b, truth, weights=w_b, estimator="empirical")
+        self.assertAlmostEqual(a[0, 0], b[0, 0], places=9)
+
+
+class TestMarginalCrpsEquivariance(unittest.TestCase):
+    def test_shift_equivariance(self):
+        rng = np.random.default_rng(5)
+        samples = rng.normal(size=(2, 5, 2))
+        truth = rng.normal(size=(2, 2))
+        shift = 137.0
+        base = marginal_crps(samples, truth, estimator="empirical")
+        shifted = marginal_crps(samples + shift, truth + shift, estimator="empirical")
+        np.testing.assert_allclose(base, shifted, atol=1e-9)
+
+    def test_positive_scale_equivariance(self):
+        rng = np.random.default_rng(6)
+        samples = rng.normal(size=(2, 5, 2))
+        truth = rng.normal(size=(2, 2))
+        a = 3.5
+        base = marginal_crps(samples, truth, estimator="empirical")
+        scaled = marginal_crps(samples * a, truth * a, estimator="empirical")
+        np.testing.assert_allclose(scaled, base * a, atol=1e-9)
+
+
+class TestMarginalCrpsEdgeCasesAndValidation(unittest.TestCase):
+    def test_duplicate_draws_no_crash_and_correct(self):
+        samples = np.array([[[1.0], [1.0], [1.0], [1.0]]])  # all identical draws
+        truth = np.array([[2.0]])
+        result = marginal_crps(samples, truth, estimator="empirical")
+        # a point mass at 1 scored against truth 2: CRPS = |1-2| = 1
+        self.assertAlmostEqual(result[0, 0], 1.0, places=12)
+
+    def test_one_draw_valid_for_empirical(self):
+        samples = np.array([[[5.0]]])
+        truth = np.array([[2.0]])
+        result = marginal_crps(samples, truth, estimator="empirical")
+        self.assertAlmostEqual(result[0, 0], 3.0, places=12)
+
+    def test_one_draw_invalid_for_unbiased(self):
+        samples = np.array([[[5.0]]])
+        truth = np.array([[2.0]])
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, estimator="unbiased")
+
+    def test_unbiased_is_provably_nonnegative_not_possibly_negative(self):
+        # For a fixed truth value y, the triangle inequality ensures that the
+        # finite-sample iid U-estimator is non-negative in exact arithmetic.
+        # Include a tight cluster plus one distant point as a stress case.
+        samples = np.array([[[0.0], [0.0], [100.0]]])
+        truth = np.array([[0.0]])
+        result = marginal_crps(samples, truth, estimator="unbiased")
+        self.assertGreaterEqual(result[0, 0], 0.0)
+
+        rng = np.random.default_rng(2026)
+        for _ in range(200):
+            m = int(rng.integers(2, 9))
+            x = rng.normal(scale=rng.uniform(0.1, 5.0), size=m)
+            y = rng.normal(scale=rng.uniform(0.1, 10.0))
+            value = marginal_crps(x.reshape(1, m, 1), np.array([[y]]), estimator="unbiased")[0, 0]
+            self.assertGreaterEqual(
+                value, -1e-9, msg=f"unbiased CRPS went negative for x={x!r}, y={y!r}"
+            )
+
+    def test_unbiased_rejects_weights(self):
+        samples = np.zeros((1, 3, 1))
+        truth = np.zeros((1, 1))
+        weights = np.full((1, 3), 1.0 / 3.0)
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, weights=weights, estimator="unbiased")
+
+    def test_all_zero_weights_row_rejected(self):
+        samples = np.zeros((1, 3, 1))
+        truth = np.zeros((1, 1))
+        weights = np.zeros((1, 3))
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, weights=weights)
+
+    def test_negative_weights_rejected(self):
+        samples = np.zeros((1, 3, 1))
+        truth = np.zeros((1, 1))
+        weights = np.array([[0.5, -0.1, 0.6]])
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, weights=weights)
+
+    def test_nonfinite_weights_rejected(self):
+        samples = np.zeros((1, 3, 1))
+        truth = np.zeros((1, 1))
+        weights = np.array([[0.5, np.nan, 0.6]])
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, weights=weights)
+
+    def test_nonfinite_samples_rejected(self):
+        samples = np.array([[[0.0], [np.inf], [1.0]]])
+        truth = np.zeros((1, 1))
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth)
+
+    def test_empty_axes_rejected(self):
+        with self.assertRaises(ValueError):
+            marginal_crps(np.zeros((0, 3, 1)), np.zeros((0, 1)))
+        with self.assertRaises(ValueError):
+            marginal_crps(np.zeros((2, 0, 1)), np.zeros((2, 1)))
+        with self.assertRaises(ValueError):
+            marginal_crps(np.zeros((2, 3, 0)), np.zeros((2, 0)))
+
+    def test_invalid_estimator_name_rejected(self):
+        samples = np.zeros((1, 3, 1))
+        truth = np.zeros((1, 1))
+        with self.assertRaises(ValueError):
+            marginal_crps(samples, truth, estimator="not-a-real-estimator")
+
+
+# --------------------------------------------------------------------------
+# brier_score
+# --------------------------------------------------------------------------
+
+
+class TestBrierScore(unittest.TestCase):
+    def test_exact_binary_case(self):
+        probs = np.array([[0.7, 0.3], [0.2, 0.8]])
+        labels = np.array([0, 1])
+        result = brier_score(probs, labels)
+        # case 0: (0.7-1)^2 + (0.3-0)^2 = 0.09+0.09=0.18
+        # case 1: (0.2-0)^2 + (0.8-1)^2 = 0.04+0.04=0.08
+        np.testing.assert_allclose(result, [0.18, 0.08], atol=1e-12)
+
+    def test_binary_convention_is_twice_single_probability_convention(self):
+        p_positive = 0.63
+        probs = np.array([[1 - p_positive, p_positive]])
+        labels = np.array([1])
+        result = brier_score(probs, labels)
+        single_convention = (p_positive - 1.0) ** 2
+        self.assertAlmostEqual(result[0], 2.0 * single_convention, places=12)
+
+    def test_random_against_oracle(self):
+        rng = np.random.default_rng(42)
+        cases, classes = 6, 4
+        raw = rng.uniform(0.01, 1.0, size=(cases, classes))
+        probs = raw / raw.sum(axis=1, keepdims=True)
+        labels = rng.integers(0, classes, size=cases)
+        result = brier_score(probs, labels)
+        expected = _oracle_brier(probs, labels)
+        np.testing.assert_allclose(result, expected, atol=1e-9)
+
+    def test_rows_not_summing_to_one_rejected(self):
+        probs = np.array([[0.5, 0.4]])  # sums to 0.9
+        labels = np.array([0])
+        with self.assertRaises(ValueError):
+            brier_score(probs, labels)
+
+    def test_negative_probability_rejected(self):
+        probs = np.array([[1.2, -0.2]])
+        labels = np.array([0])
+        with self.assertRaises(ValueError):
+            brier_score(probs, labels)
+
+    def test_non_integer_labels_rejected(self):
+        probs = np.array([[0.5, 0.5]])
+        labels = np.array([0.0])  # float, not int
+        with self.assertRaises(TypeError):
+            brier_score(probs, labels)
+
+    def test_out_of_range_label_rejected(self):
+        probs = np.array([[0.5, 0.5]])
+        labels = np.array([2])  # only classes 0,1 exist
+        with self.assertRaises(ValueError):
+            brier_score(probs, labels)
+
+    def test_shape_mismatch_rejected(self):
+        probs = np.array([[0.5, 0.5], [0.3, 0.7]])
+        labels = np.array([0])  # only one label for two cases
+        with self.assertRaises(ValueError):
+            brier_score(probs, labels)
+
+    def test_empty_axes_rejected(self):
+        with self.assertRaises(ValueError):
+            brier_score(np.zeros((0, 2)), np.zeros((0,), dtype=int))
+
+
+# --------------------------------------------------------------------------
+# expected_action_losses / bayes_action_set
+# --------------------------------------------------------------------------
+
+
+class TestExpectedActionLossesAndBayesActionSet(unittest.TestCase):
+    def test_uniform_weights_matches_simple_mean(self):
+        losses = np.array([[1.0, 4.0], [3.0, 2.0], [5.0, 0.0]])
+        result = expected_action_losses(losses)
+        expected = losses.mean(axis=0)
+        np.testing.assert_allclose(result, expected, atol=1e-12)
+
+    def test_nonuniform_weights_against_oracle(self):
+        rng = np.random.default_rng(123)
+        draws, actions = 7, 3
+        losses = rng.normal(size=(draws, actions))
+        weights = rng.uniform(0.1, 2.0, size=draws)
+        result = expected_action_losses(losses, weights=weights)
+        expected = _oracle_expected_action_losses(losses, weights)
+        np.testing.assert_allclose(result, expected, atol=1e-9)
+
+    def test_single_best_action(self):
+        losses = np.array([[1.0, 2.0, 3.0], [1.0, 2.0, 3.0]])
+        min_risk, minimizers = bayes_action_set(losses)
+        self.assertAlmostEqual(min_risk, 1.0, places=12)
+        np.testing.assert_array_equal(minimizers, np.array([0]))
+
+    def test_exact_tie_returns_all_minimizers(self):
+        # actions 0 and 2 tie exactly at expected loss 2.0
+        losses = np.array([[2.0, 5.0, 2.0], [2.0, 5.0, 2.0], [2.0, 5.0, 2.0]])
+        min_risk, minimizers = bayes_action_set(losses)
+        self.assertAlmostEqual(min_risk, 2.0, places=12)
+        np.testing.assert_array_equal(np.sort(minimizers), np.array([0, 2]))
+
+    def test_exact_tie_with_nonuniform_weights(self):
+        # Constructed so actions 0 and 1 have identical expected loss under
+        # nonuniform weights, action 2 is strictly worse.
+        weights = np.array([0.2, 0.8])
+        losses = np.array(
+            [
+                [10.0, 0.0, 20.0],
+                [0.0, 2.5, 20.0],
+            ]
+        )
+        # action0: 0.2*10+0.8*0=2.0 ; action1: 0.2*0+0.8*2.5=2.0 ; action2: 20
+        min_risk, minimizers = bayes_action_set(losses, weights=weights)
+        self.assertAlmostEqual(min_risk, 2.0, places=12)
+        np.testing.assert_array_equal(np.sort(minimizers), np.array([0, 1]))
+
+    def test_losses_wrong_ndim_rejected(self):
+        with self.assertRaises(ValueError):
+            expected_action_losses(np.zeros(5))
+
+    def test_empty_axes_rejected(self):
+        with self.assertRaises(ValueError):
+            expected_action_losses(np.zeros((0, 3)))
+        with self.assertRaises(ValueError):
+            expected_action_losses(np.zeros((3, 0)))
+
+    def test_all_zero_weights_rejected(self):
+        losses = np.zeros((3, 2))
+        with self.assertRaises(ValueError):
+            expected_action_losses(losses, weights=np.zeros(3))
+
+    def test_huge_weights_are_rescaled_for_action_losses_and_ties(self):
+        losses = np.array([[1.0, 3.0], [1.0, 3.0]])
+        weights = np.array([1e308, 1e308])
+        np.testing.assert_allclose(expected_action_losses(losses, weights), [1.0, 3.0])
+        min_risk, minimizers = bayes_action_set(losses, weights)
+        self.assertEqual(min_risk, 1.0)
+        np.testing.assert_array_equal(minimizers, np.array([0]))
+
+    def test_negative_weights_rejected(self):
+        losses = np.zeros((3, 2))
+        with self.assertRaises(ValueError):
+            expected_action_losses(losses, weights=np.array([0.5, -0.1, 0.6]))
+
+    def test_nonfinite_losses_rejected(self):
+        losses = np.array([[1.0, np.nan], [2.0, 3.0]])
+        with self.assertRaises(ValueError):
+            expected_action_losses(losses)
+
+    def test_weight_shape_mismatch_rejected(self):
+        losses = np.zeros((3, 2))
+        with self.assertRaises(ValueError):
+            expected_action_losses(losses, weights=np.zeros(4))
+
+
+if __name__ == "__main__":
+    unittest.main()