Skip to content

feat: add energy_score() multivariate probabilistic metric - #5

Open
GiGiKoneti wants to merge 6 commits into
mllam:mainfrom
GiGiKoneti:feat/energy-score
Open

feat: add energy_score() multivariate probabilistic metric#5
GiGiKoneti wants to merge 6 commits into
mllam:mainfrom
GiGiKoneti:feat/energy-score

Conversation

@GiGiKoneti

Copy link
Copy Markdown

Describe your changes

adds energy_score() to mllam_verification/operations/statistics.py
right after spread_skill_ratio(). the energy score is a multivariate
proper scoring rule that generalises CRPS to vector-valued forecasts —
useful for checking spatial coherence across ensemble members, not just
marginal correctness at individual grid points.

implemented from scratch since scores>=1.2.0 doesn't include it:

ES = E||X - y||₂ - 0.5 · E||X - X'||₂

where X, X' are independent ensemble draws and y is the observation.
used the xarray rename trick (ensemble_member_prime) for clean O(M²)
pairwise computation without python loops. mirrored the
update_cell_methods() call and overall structure from crps() exactly
so it fits the existing pattern.

no new dependencies — formula implemented using xarray ops already
present in the codebase.

Issue Link

related to #3 (probabilistic evaluation coverage) — energy score is
the natural next metric after crps and spread_skill_ratio from PR #4

Type of change

  • 🐛 Bug fix (non-breaking change that fixes an issue)
  • ✨ New feature (non-breaking change that adds functionality)
  • 💥 Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • 📖 Documentation (Addition or improvements to documentation)

Checklist before requesting a review

  • My branch is up-to-date with the target branch
  • I have performed a self-review of my code
  • For any new/modified functions/classes I have added docstrings
    that clearly describe its purpose, expected inputs and returned values
  • I have placed in-line comments to clarify the intent of any
    hard-to-understand passages of my code
  • I have updated the documentation to cover introduced code changes
  • I have added tests that prove my fix is effective or that my feature works
  • I have given the PR a name that clearly describes the change,
    written in imperative form
  • I have requested a reviewer and an assignee

Checklist for reviewers

  • the code is readable
  • the code is well tested
  • the code is documented (including return types and parameters)
  • the code is easy to maintain

Author checklist after completed review

  • I have added a line to the CHANGELOG describing this change,
    in a section reflecting type of change:
    - added: energy_score() multivariate proper scoring rule

Checklist for assignee

  • PR is up to date with the base branch
  • the tests pass
  • author has added an entry to the changelog

…semble evaluation

Adds probabilistic ensemble evaluation metrics to mllam-verification,
as proposed in issue mllam#3 and approved by @mfroelund.

New functions in statistics.py:
- crps(): wraps scores.probability.crps_for_ensemble via
  compute_pipeline_statistic, following the same pattern as rmse()
  and mae(). Uses the fair (unbiased) estimator. Accepts any
  ensemble_member_dim name.
- spread_skill_ratio(): computes ensemble spread / RMSE of ensemble
  mean. SSR = 1.0 indicates perfect calibration. SSR < 1.0 indicates
  underdispersion.

New function in plot.py:
- plot_rank_histogram(): wraps scores.plotdata.rank_histogram to
  produce a Talagrand diagram. Includes a reference line for perfect
  calibration.

New test fixture in conftest.py:
- da_ensemble_prediction_2d_utc: 10-member ensemble DataArray built
  from the existing deterministic prediction fixture.

All functions follow the existing architecture exactly:
compute_pipeline_statistic backbone, Google-style docstrings,
cell_methods CF-convention attribute, 90-char line length.

No new dependencies introduced. All functions use the existing
scores>=1.2.0 dependency already pinned in pyproject.toml.

Closes mllam#3
…ove imports

- Remove unused groupby parameter from crps() and spread_skill_ratio()
- Swap argument order to (ds_reference, ds_prediction) matching plot
  calling convention and mae() signature
- Add preserve_dims support to spread_skill_ratio() for hovmoller plots
- Move matplotlib/plot_rank_histogram imports to top of test_plot.py
- Add crps and spread_skill_ratio to timeseries/hovmoller test
  parametrizations
- Add da_ensemble_prediction_2d_elapsed fixture for elapsed-time tests
- Update test_statistics.py to match new argument order

All 38 tests pass. Pre-commit hooks (isort, black, flake8, mypy) clean.
Implements the Energy Score multivariate proper scoring
rule from scratch (not available in scores>=1.2.0).

Formula: ES = E||X-Y||₂ - 0.5·E||X-X'||₂

Follows exact same signature and convention as crps()
and spread_skill_ratio(). Four unit tests added to
TestEnergyScore class in test_statistics.py:
- returns correct type
- ensemble_member dim collapsed in output
- perfect ensemble produces ES=0 (mathematical ground truth)
- non-negative for any input (ES >= 0 by definition)

Cite: Gneiting & Raftery (2007), JASA 102(477)
      Rasp et al. (2024), WeatherBench 2, JAMES

Relevant to neural-lam Issue #62 and GSoC Project 3
probabilistic evaluation pipeline.
@GiGiKoneti

Copy link
Copy Markdown
Author

Note: This PR depends on #4 and is rebased on top of it. The diff currently includes #4's changes because it hasn't merged into main yet. Once #4 is merged, this PR's diff will automatically reduce to just the energy_score() addition and its tests.

@mfroelund

Copy link
Copy Markdown
Collaborator

I just looked through the PR, and your implementation looks great. I'll wait for #4 to get merged before I submit my review.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants