diff --git a/docs/crps_estimators.md b/docs/background/crps_estimators.md similarity index 100% rename from docs/crps_estimators.md rename to docs/background/crps_estimators.md diff --git a/docs/forecast_dists.md b/docs/background/forecast_dists.md similarity index 100% rename from docs/forecast_dists.md rename to docs/background/forecast_dists.md diff --git a/docs/theory.md b/docs/background/theory.md similarity index 100% rename from docs/theory.md rename to docs/background/theory.md diff --git a/docs/weighted_scores.md b/docs/background/weighted_scores.md similarity index 100% rename from docs/weighted_scores.md rename to docs/background/weighted_scores.md diff --git a/docs/contributing.md b/docs/contributing.md deleted file mode 100644 index 3294e06..0000000 --- a/docs/contributing.md +++ /dev/null @@ -1,48 +0,0 @@ -# Contributing - -We welcome contributions! You can help improve the library in many ways: - -- Report issues or propose enhancements using on the [GitHub issue tracker](https://github.com/frazane/scoringrules/issues) -- Improve or extend the codebase -- Improve or extend the documentation - -## Getting started - -Fork the repository on GitHub and clone it to your computer: - -``` -git clone https://github.com//scoringrules.git -``` - -We use [uv](https://docs.astral.sh/uv/) for project management and packaging. Install it with - -``` -curl -LsSf https://astral.sh/uv/install.sh | sh -``` - -Then, you can install the library and all dependencies (including development dependencies) and install the pre-commit hooks: - -``` -uv install -uv run pre-commit install -``` - -From here you can work on your changes! Once you're satisfied with your changes, and followed the additional instructions below, push everything to your repository and open a pull request on GitHub. - - -### Contributing to the codebase -Don't forget to include new tests if necessary, then make sure that all tests are passing with - -``` -uv run pytest tests/ -``` - -### Contributing to the documentation - -You can work on the documentation by modifying files in `docs/`. The most convenient way to do it is to run - -``` -uvx --with-requirements docs/requirements.txt sphinx-autobuild docs/ docs/_build/ -``` - -and open the locally hosted documentation on your browser. It will be updated automatically every time you make changes and save. If you edit or add pieces of LaTex math, please make sure they are rendered correctly. diff --git a/docs/index.md b/docs/index.md index 739f9c0..a1174c5 100644 --- a/docs/index.md +++ b/docs/index.md @@ -17,7 +17,6 @@ The scoring rules available in `scoringrules` include, but are not limited to, t - Logarithmic Score - (Discrete) Ranked Probability Score - Continuous Ranked Probability Score (CRPS) -- Dawid-Sebastiani Score - Energy Score - Variogram Score - Gaussian Kernel Score @@ -25,27 +24,27 @@ The scoring rules available in `scoringrules` include, but are not limited to, t - Threshold-Weighted Energy Score - Quantile Score - Interval Score +- Weighted Interval Score - - +Functionality is available for forecasts in the form of samples (i.e. ensemble forecasts) as well as popular univariate parametric distributions. ```{toctree} :hidden: :caption: Background -theory.md -forecast_dists.md -crps_estimators.md -weighted_scores.md +background/theory.md +background/forecast_dists.md +background/crps_estimators.md +background/weighted_scores.md ``` ```{toctree} :hidden: :caption: Library -user_guide.md -contributing.md -reference.md +library/user_guide.md +library/contributing.md +library/reference.md ``` diff --git a/docs/library/contributing.md b/docs/library/contributing.md new file mode 100644 index 0000000..8784c24 --- /dev/null +++ b/docs/library/contributing.md @@ -0,0 +1,49 @@ +# Contributing + +We welcome contributions! You can help improve the library in many ways. For example, by: + +- Reporting issues or proposing enhancements via the [GitHub issue tracker](https://github.com/frazane/scoringrules/issues) +- Improving or extending the codebase +- Improving or extending the documentation + +## Getting started + +To make changes to the library, fork the repository on GitHub and clone it to your computer: + +``` +git clone https://github.com//scoringrules.git +``` + +We use [uv](https://docs.astral.sh/uv/) for project management and packaging. Install it by running + +``` +curl -LsSf https://astral.sh/uv/install.sh | sh +``` + +You should then be able to install the library and all dependencies (including development dependencies), and install the pre-commit hooks: + +``` +uv install +uv run pre-commit install +``` + +From here, you can work on your changes. Once you're satisfied with your changes, and have followed the additional instructions below, push everything to your repository and open a pull request on GitHub. + + +### Contributing to the codebase + +Don't forget to include new tests if necessary. Make sure that all tests are passing with + +``` +uv run pytest tests/ +``` + +### Contributing to the documentation + +You can work on the documentation by modifying files in `docs/`. The most convenient way to do this is to run + +``` +uvx --with-requirements docs/requirements.txt sphinx-autobuild docs/ docs/_build/ +``` + +and open the locally hosted documentation on your browser. It will be updated automatically every time you make changes and save. If you edit or add pieces of LaTeX math, please make sure they are rendered correctly. diff --git a/docs/library/generated/scoringrules.brier_score.rst b/docs/library/generated/scoringrules.brier_score.rst new file mode 100644 index 0000000..c821784 --- /dev/null +++ b/docs/library/generated/scoringrules.brier_score.rst @@ -0,0 +1,6 @@ +scoringrules.brier\_score +========================= + +.. currentmodule:: scoringrules + +.. autofunction:: brier_score diff --git a/docs/library/generated/scoringrules.clogs_ensemble.rst b/docs/library/generated/scoringrules.clogs_ensemble.rst new file mode 100644 index 0000000..fcb2b79 --- /dev/null +++ b/docs/library/generated/scoringrules.clogs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.clogs\_ensemble +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: clogs_ensemble diff --git a/docs/library/generated/scoringrules.crps_2pexponential.rst b/docs/library/generated/scoringrules.crps_2pexponential.rst new file mode 100644 index 0000000..742b5a8 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_2pexponential.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_2pexponential +================================ + +.. currentmodule:: scoringrules + +.. autofunction:: crps_2pexponential diff --git a/docs/library/generated/scoringrules.crps_2pnormal.rst b/docs/library/generated/scoringrules.crps_2pnormal.rst new file mode 100644 index 0000000..56e42fa --- /dev/null +++ b/docs/library/generated/scoringrules.crps_2pnormal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_2pnormal +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_2pnormal diff --git a/docs/library/generated/scoringrules.crps_beta.rst b/docs/library/generated/scoringrules.crps_beta.rst new file mode 100644 index 0000000..fd41a3e --- /dev/null +++ b/docs/library/generated/scoringrules.crps_beta.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_beta +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: crps_beta diff --git a/docs/library/generated/scoringrules.crps_binomial.rst b/docs/library/generated/scoringrules.crps_binomial.rst new file mode 100644 index 0000000..cfa1350 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_binomial.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_binomial +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_binomial diff --git a/docs/library/generated/scoringrules.crps_clogistic.rst b/docs/library/generated/scoringrules.crps_clogistic.rst new file mode 100644 index 0000000..b5a658d --- /dev/null +++ b/docs/library/generated/scoringrules.crps_clogistic.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_clogistic +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: crps_clogistic diff --git a/docs/library/generated/scoringrules.crps_cnormal.rst b/docs/library/generated/scoringrules.crps_cnormal.rst new file mode 100644 index 0000000..0739fbc --- /dev/null +++ b/docs/library/generated/scoringrules.crps_cnormal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_cnormal +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_cnormal diff --git a/docs/library/generated/scoringrules.crps_ct.rst b/docs/library/generated/scoringrules.crps_ct.rst new file mode 100644 index 0000000..70fb9ba --- /dev/null +++ b/docs/library/generated/scoringrules.crps_ct.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_ct +===================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_ct diff --git a/docs/library/generated/scoringrules.crps_ensemble.rst b/docs/library/generated/scoringrules.crps_ensemble.rst new file mode 100644 index 0000000..557690f --- /dev/null +++ b/docs/library/generated/scoringrules.crps_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_ensemble diff --git a/docs/library/generated/scoringrules.crps_exponential.rst b/docs/library/generated/scoringrules.crps_exponential.rst new file mode 100644 index 0000000..67c096b --- /dev/null +++ b/docs/library/generated/scoringrules.crps_exponential.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_exponential +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_exponential diff --git a/docs/library/generated/scoringrules.crps_exponentialM.rst b/docs/library/generated/scoringrules.crps_exponentialM.rst new file mode 100644 index 0000000..a13250d --- /dev/null +++ b/docs/library/generated/scoringrules.crps_exponentialM.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_exponentialM +=============================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_exponentialM diff --git a/docs/library/generated/scoringrules.crps_gamma.rst b/docs/library/generated/scoringrules.crps_gamma.rst new file mode 100644 index 0000000..ff552fc --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gamma.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gamma +======================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gamma diff --git a/docs/library/generated/scoringrules.crps_gev.rst b/docs/library/generated/scoringrules.crps_gev.rst new file mode 100644 index 0000000..262103e --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gev.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gev +====================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gev diff --git a/docs/library/generated/scoringrules.crps_gpd.rst b/docs/library/generated/scoringrules.crps_gpd.rst new file mode 100644 index 0000000..ff78341 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gpd.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gpd +====================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gpd diff --git a/docs/library/generated/scoringrules.crps_gtclogistic.rst b/docs/library/generated/scoringrules.crps_gtclogistic.rst new file mode 100644 index 0000000..3f2f82d --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gtclogistic.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gtclogistic +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gtclogistic diff --git a/docs/library/generated/scoringrules.crps_gtcnormal.rst b/docs/library/generated/scoringrules.crps_gtcnormal.rst new file mode 100644 index 0000000..3e75a79 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gtcnormal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gtcnormal +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gtcnormal diff --git a/docs/library/generated/scoringrules.crps_gtct.rst b/docs/library/generated/scoringrules.crps_gtct.rst new file mode 100644 index 0000000..89735e6 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_gtct.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_gtct +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: crps_gtct diff --git a/docs/library/generated/scoringrules.crps_hypergeometric.rst b/docs/library/generated/scoringrules.crps_hypergeometric.rst new file mode 100644 index 0000000..6c7cc1b --- /dev/null +++ b/docs/library/generated/scoringrules.crps_hypergeometric.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_hypergeometric +================================= + +.. currentmodule:: scoringrules + +.. autofunction:: crps_hypergeometric diff --git a/docs/library/generated/scoringrules.crps_laplace.rst b/docs/library/generated/scoringrules.crps_laplace.rst new file mode 100644 index 0000000..481596e --- /dev/null +++ b/docs/library/generated/scoringrules.crps_laplace.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_laplace +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_laplace diff --git a/docs/library/generated/scoringrules.crps_logistic.rst b/docs/library/generated/scoringrules.crps_logistic.rst new file mode 100644 index 0000000..687b73a --- /dev/null +++ b/docs/library/generated/scoringrules.crps_logistic.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_logistic +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_logistic diff --git a/docs/library/generated/scoringrules.crps_loglaplace.rst b/docs/library/generated/scoringrules.crps_loglaplace.rst new file mode 100644 index 0000000..d6a0bee --- /dev/null +++ b/docs/library/generated/scoringrules.crps_loglaplace.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_loglaplace +============================= + +.. currentmodule:: scoringrules + +.. autofunction:: crps_loglaplace diff --git a/docs/library/generated/scoringrules.crps_loglogistic.rst b/docs/library/generated/scoringrules.crps_loglogistic.rst new file mode 100644 index 0000000..2107a2e --- /dev/null +++ b/docs/library/generated/scoringrules.crps_loglogistic.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_loglogistic +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_loglogistic diff --git a/docs/library/generated/scoringrules.crps_lognormal.rst b/docs/library/generated/scoringrules.crps_lognormal.rst new file mode 100644 index 0000000..fd49df7 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_lognormal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_lognormal +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: crps_lognormal diff --git a/docs/library/generated/scoringrules.crps_mixnorm.rst b/docs/library/generated/scoringrules.crps_mixnorm.rst new file mode 100644 index 0000000..8fba00b --- /dev/null +++ b/docs/library/generated/scoringrules.crps_mixnorm.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_mixnorm +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_mixnorm diff --git a/docs/library/generated/scoringrules.crps_negbinom.rst b/docs/library/generated/scoringrules.crps_negbinom.rst new file mode 100644 index 0000000..a1d12d2 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_negbinom.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_negbinom +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_negbinom diff --git a/docs/library/generated/scoringrules.crps_normal.rst b/docs/library/generated/scoringrules.crps_normal.rst new file mode 100644 index 0000000..a670094 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_normal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_normal +========================= + +.. currentmodule:: scoringrules + +.. autofunction:: crps_normal diff --git a/docs/library/generated/scoringrules.crps_poisson.rst b/docs/library/generated/scoringrules.crps_poisson.rst new file mode 100644 index 0000000..aca7e0e --- /dev/null +++ b/docs/library/generated/scoringrules.crps_poisson.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_poisson +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_poisson diff --git a/docs/library/generated/scoringrules.crps_quantile.rst b/docs/library/generated/scoringrules.crps_quantile.rst new file mode 100644 index 0000000..1db2b18 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_quantile.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_quantile +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_quantile diff --git a/docs/library/generated/scoringrules.crps_t.rst b/docs/library/generated/scoringrules.crps_t.rst new file mode 100644 index 0000000..c9ad9fc --- /dev/null +++ b/docs/library/generated/scoringrules.crps_t.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_t +==================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_t diff --git a/docs/library/generated/scoringrules.crps_tlogistic.rst b/docs/library/generated/scoringrules.crps_tlogistic.rst new file mode 100644 index 0000000..04699f9 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_tlogistic.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_tlogistic +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: crps_tlogistic diff --git a/docs/library/generated/scoringrules.crps_tnormal.rst b/docs/library/generated/scoringrules.crps_tnormal.rst new file mode 100644 index 0000000..939e73b --- /dev/null +++ b/docs/library/generated/scoringrules.crps_tnormal.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_tnormal +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_tnormal diff --git a/docs/library/generated/scoringrules.crps_tt.rst b/docs/library/generated/scoringrules.crps_tt.rst new file mode 100644 index 0000000..1d2dc90 --- /dev/null +++ b/docs/library/generated/scoringrules.crps_tt.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_tt +===================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_tt diff --git a/docs/library/generated/scoringrules.crps_uniform.rst b/docs/library/generated/scoringrules.crps_uniform.rst new file mode 100644 index 0000000..1ed6d9f --- /dev/null +++ b/docs/library/generated/scoringrules.crps_uniform.rst @@ -0,0 +1,6 @@ +scoringrules.crps\_uniform +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: crps_uniform diff --git a/docs/library/generated/scoringrules.es_ensemble.rst b/docs/library/generated/scoringrules.es_ensemble.rst new file mode 100644 index 0000000..efbd76a --- /dev/null +++ b/docs/library/generated/scoringrules.es_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.es\_ensemble +========================= + +.. currentmodule:: scoringrules + +.. autofunction:: es_ensemble diff --git a/docs/library/generated/scoringrules.gksmv_ensemble.rst b/docs/library/generated/scoringrules.gksmv_ensemble.rst new file mode 100644 index 0000000..ae5873f --- /dev/null +++ b/docs/library/generated/scoringrules.gksmv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.gksmv\_ensemble +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: gksmv_ensemble diff --git a/docs/library/generated/scoringrules.gksuv_ensemble.rst b/docs/library/generated/scoringrules.gksuv_ensemble.rst new file mode 100644 index 0000000..1792921 --- /dev/null +++ b/docs/library/generated/scoringrules.gksuv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.gksuv\_ensemble +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: gksuv_ensemble diff --git a/docs/library/generated/scoringrules.interval_score.rst b/docs/library/generated/scoringrules.interval_score.rst new file mode 100644 index 0000000..48f50c5 --- /dev/null +++ b/docs/library/generated/scoringrules.interval_score.rst @@ -0,0 +1,6 @@ +scoringrules.interval\_score +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: interval_score diff --git a/docs/library/generated/scoringrules.log_score.rst b/docs/library/generated/scoringrules.log_score.rst new file mode 100644 index 0000000..831a334 --- /dev/null +++ b/docs/library/generated/scoringrules.log_score.rst @@ -0,0 +1,6 @@ +scoringrules.log\_score +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: log_score diff --git a/docs/library/generated/scoringrules.logs_2pexponential.rst b/docs/library/generated/scoringrules.logs_2pexponential.rst new file mode 100644 index 0000000..6639250 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_2pexponential.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_2pexponential +================================ + +.. currentmodule:: scoringrules + +.. autofunction:: logs_2pexponential diff --git a/docs/library/generated/scoringrules.logs_2pnormal.rst b/docs/library/generated/scoringrules.logs_2pnormal.rst new file mode 100644 index 0000000..78308b7 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_2pnormal.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_2pnormal +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_2pnormal diff --git a/docs/library/generated/scoringrules.logs_beta.rst b/docs/library/generated/scoringrules.logs_beta.rst new file mode 100644 index 0000000..fc7cf26 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_beta.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_beta +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: logs_beta diff --git a/docs/library/generated/scoringrules.logs_binomial.rst b/docs/library/generated/scoringrules.logs_binomial.rst new file mode 100644 index 0000000..726d986 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_binomial.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_binomial +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_binomial diff --git a/docs/library/generated/scoringrules.logs_ensemble.rst b/docs/library/generated/scoringrules.logs_ensemble.rst new file mode 100644 index 0000000..666be79 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_ensemble diff --git a/docs/library/generated/scoringrules.logs_exponential.rst b/docs/library/generated/scoringrules.logs_exponential.rst new file mode 100644 index 0000000..cbd8559 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_exponential.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_exponential +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_exponential diff --git a/docs/library/generated/scoringrules.logs_exponential2.rst b/docs/library/generated/scoringrules.logs_exponential2.rst new file mode 100644 index 0000000..8b9b636 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_exponential2.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_exponential2 +=============================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_exponential2 diff --git a/docs/library/generated/scoringrules.logs_gamma.rst b/docs/library/generated/scoringrules.logs_gamma.rst new file mode 100644 index 0000000..ee15b3d --- /dev/null +++ b/docs/library/generated/scoringrules.logs_gamma.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_gamma +======================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_gamma diff --git a/docs/library/generated/scoringrules.logs_gev.rst b/docs/library/generated/scoringrules.logs_gev.rst new file mode 100644 index 0000000..0284bc0 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_gev.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_gev +====================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_gev diff --git a/docs/library/generated/scoringrules.logs_gpd.rst b/docs/library/generated/scoringrules.logs_gpd.rst new file mode 100644 index 0000000..9b8e0e5 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_gpd.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_gpd +====================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_gpd diff --git a/docs/library/generated/scoringrules.logs_hypergeometric.rst b/docs/library/generated/scoringrules.logs_hypergeometric.rst new file mode 100644 index 0000000..46095fb --- /dev/null +++ b/docs/library/generated/scoringrules.logs_hypergeometric.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_hypergeometric +================================= + +.. currentmodule:: scoringrules + +.. autofunction:: logs_hypergeometric diff --git a/docs/library/generated/scoringrules.logs_laplace.rst b/docs/library/generated/scoringrules.logs_laplace.rst new file mode 100644 index 0000000..74538c6 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_laplace.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_laplace +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_laplace diff --git a/docs/library/generated/scoringrules.logs_logistic.rst b/docs/library/generated/scoringrules.logs_logistic.rst new file mode 100644 index 0000000..f8d0639 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_logistic.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_logistic +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_logistic diff --git a/docs/library/generated/scoringrules.logs_loglaplace.rst b/docs/library/generated/scoringrules.logs_loglaplace.rst new file mode 100644 index 0000000..3cf835b --- /dev/null +++ b/docs/library/generated/scoringrules.logs_loglaplace.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_loglaplace +============================= + +.. currentmodule:: scoringrules + +.. autofunction:: logs_loglaplace diff --git a/docs/library/generated/scoringrules.logs_loglogistic.rst b/docs/library/generated/scoringrules.logs_loglogistic.rst new file mode 100644 index 0000000..a0ea488 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_loglogistic.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_loglogistic +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_loglogistic diff --git a/docs/library/generated/scoringrules.logs_lognormal.rst b/docs/library/generated/scoringrules.logs_lognormal.rst new file mode 100644 index 0000000..ef016c0 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_lognormal.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_lognormal +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: logs_lognormal diff --git a/docs/library/generated/scoringrules.logs_mixnorm.rst b/docs/library/generated/scoringrules.logs_mixnorm.rst new file mode 100644 index 0000000..f56a9d7 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_mixnorm.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_mixnorm +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_mixnorm diff --git a/docs/library/generated/scoringrules.logs_negbinom.rst b/docs/library/generated/scoringrules.logs_negbinom.rst new file mode 100644 index 0000000..8324c66 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_negbinom.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_negbinom +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_negbinom diff --git a/docs/library/generated/scoringrules.logs_normal.rst b/docs/library/generated/scoringrules.logs_normal.rst new file mode 100644 index 0000000..beb64a9 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_normal.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_normal +========================= + +.. currentmodule:: scoringrules + +.. autofunction:: logs_normal diff --git a/docs/library/generated/scoringrules.logs_poisson.rst b/docs/library/generated/scoringrules.logs_poisson.rst new file mode 100644 index 0000000..98d6ea2 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_poisson.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_poisson +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_poisson diff --git a/docs/library/generated/scoringrules.logs_t.rst b/docs/library/generated/scoringrules.logs_t.rst new file mode 100644 index 0000000..9b11053 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_t.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_t +==================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_t diff --git a/docs/library/generated/scoringrules.logs_tlogistic.rst b/docs/library/generated/scoringrules.logs_tlogistic.rst new file mode 100644 index 0000000..e74a1c6 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_tlogistic.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_tlogistic +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: logs_tlogistic diff --git a/docs/library/generated/scoringrules.logs_tnormal.rst b/docs/library/generated/scoringrules.logs_tnormal.rst new file mode 100644 index 0000000..61df3dd --- /dev/null +++ b/docs/library/generated/scoringrules.logs_tnormal.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_tnormal +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_tnormal diff --git a/docs/library/generated/scoringrules.logs_tt.rst b/docs/library/generated/scoringrules.logs_tt.rst new file mode 100644 index 0000000..cf2a300 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_tt.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_tt +===================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_tt diff --git a/docs/library/generated/scoringrules.logs_uniform.rst b/docs/library/generated/scoringrules.logs_uniform.rst new file mode 100644 index 0000000..03fad84 --- /dev/null +++ b/docs/library/generated/scoringrules.logs_uniform.rst @@ -0,0 +1,6 @@ +scoringrules.logs\_uniform +========================== + +.. currentmodule:: scoringrules + +.. autofunction:: logs_uniform diff --git a/docs/library/generated/scoringrules.owcrps_ensemble.rst b/docs/library/generated/scoringrules.owcrps_ensemble.rst new file mode 100644 index 0000000..860679b --- /dev/null +++ b/docs/library/generated/scoringrules.owcrps_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.owcrps\_ensemble +============================= + +.. currentmodule:: scoringrules + +.. autofunction:: owcrps_ensemble diff --git a/docs/library/generated/scoringrules.owes_ensemble.rst b/docs/library/generated/scoringrules.owes_ensemble.rst new file mode 100644 index 0000000..9f46986 --- /dev/null +++ b/docs/library/generated/scoringrules.owes_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.owes\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: owes_ensemble diff --git a/docs/library/generated/scoringrules.owgksmv_ensemble.rst b/docs/library/generated/scoringrules.owgksmv_ensemble.rst new file mode 100644 index 0000000..8277dab --- /dev/null +++ b/docs/library/generated/scoringrules.owgksmv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.owgksmv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: owgksmv_ensemble diff --git a/docs/library/generated/scoringrules.owgksuv_ensemble.rst b/docs/library/generated/scoringrules.owgksuv_ensemble.rst new file mode 100644 index 0000000..62a7130 --- /dev/null +++ b/docs/library/generated/scoringrules.owgksuv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.owgksuv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: owgksuv_ensemble diff --git a/docs/library/generated/scoringrules.owvs_ensemble.rst b/docs/library/generated/scoringrules.owvs_ensemble.rst new file mode 100644 index 0000000..38a1765 --- /dev/null +++ b/docs/library/generated/scoringrules.owvs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.owvs\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: owvs_ensemble diff --git a/docs/library/generated/scoringrules.quantile_score.rst b/docs/library/generated/scoringrules.quantile_score.rst new file mode 100644 index 0000000..560bb96 --- /dev/null +++ b/docs/library/generated/scoringrules.quantile_score.rst @@ -0,0 +1,6 @@ +scoringrules.quantile\_score +============================ + +.. currentmodule:: scoringrules + +.. autofunction:: quantile_score diff --git a/docs/library/generated/scoringrules.register_backend.rst b/docs/library/generated/scoringrules.register_backend.rst new file mode 100644 index 0000000..0820837 --- /dev/null +++ b/docs/library/generated/scoringrules.register_backend.rst @@ -0,0 +1,6 @@ +scoringrules.register\_backend +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: register_backend diff --git a/docs/library/generated/scoringrules.rls_score.rst b/docs/library/generated/scoringrules.rls_score.rst new file mode 100644 index 0000000..1cc6638 --- /dev/null +++ b/docs/library/generated/scoringrules.rls_score.rst @@ -0,0 +1,6 @@ +scoringrules.rls\_score +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: rls_score diff --git a/docs/library/generated/scoringrules.rps_score.rst b/docs/library/generated/scoringrules.rps_score.rst new file mode 100644 index 0000000..ac8a19e --- /dev/null +++ b/docs/library/generated/scoringrules.rps_score.rst @@ -0,0 +1,6 @@ +scoringrules.rps\_score +======================= + +.. currentmodule:: scoringrules + +.. autofunction:: rps_score diff --git a/docs/library/generated/scoringrules.twcrps_ensemble.rst b/docs/library/generated/scoringrules.twcrps_ensemble.rst new file mode 100644 index 0000000..ae01854 --- /dev/null +++ b/docs/library/generated/scoringrules.twcrps_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.twcrps\_ensemble +============================= + +.. currentmodule:: scoringrules + +.. autofunction:: twcrps_ensemble diff --git a/docs/library/generated/scoringrules.twes_ensemble.rst b/docs/library/generated/scoringrules.twes_ensemble.rst new file mode 100644 index 0000000..dd247f4 --- /dev/null +++ b/docs/library/generated/scoringrules.twes_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.twes\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: twes_ensemble diff --git a/docs/library/generated/scoringrules.twgksmv_ensemble.rst b/docs/library/generated/scoringrules.twgksmv_ensemble.rst new file mode 100644 index 0000000..5118306 --- /dev/null +++ b/docs/library/generated/scoringrules.twgksmv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.twgksmv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: twgksmv_ensemble diff --git a/docs/library/generated/scoringrules.twgksuv_ensemble.rst b/docs/library/generated/scoringrules.twgksuv_ensemble.rst new file mode 100644 index 0000000..b8e0848 --- /dev/null +++ b/docs/library/generated/scoringrules.twgksuv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.twgksuv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: twgksuv_ensemble diff --git a/docs/library/generated/scoringrules.twvs_ensemble.rst b/docs/library/generated/scoringrules.twvs_ensemble.rst new file mode 100644 index 0000000..05c0f94 --- /dev/null +++ b/docs/library/generated/scoringrules.twvs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.twvs\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: twvs_ensemble diff --git a/docs/library/generated/scoringrules.vrcrps_ensemble.rst b/docs/library/generated/scoringrules.vrcrps_ensemble.rst new file mode 100644 index 0000000..c9d842d --- /dev/null +++ b/docs/library/generated/scoringrules.vrcrps_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vrcrps\_ensemble +============================= + +.. currentmodule:: scoringrules + +.. autofunction:: vrcrps_ensemble diff --git a/docs/library/generated/scoringrules.vres_ensemble.rst b/docs/library/generated/scoringrules.vres_ensemble.rst new file mode 100644 index 0000000..5e72fda --- /dev/null +++ b/docs/library/generated/scoringrules.vres_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vres\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: vres_ensemble diff --git a/docs/library/generated/scoringrules.vrgksmv_ensemble.rst b/docs/library/generated/scoringrules.vrgksmv_ensemble.rst new file mode 100644 index 0000000..326018e --- /dev/null +++ b/docs/library/generated/scoringrules.vrgksmv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vrgksmv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: vrgksmv_ensemble diff --git a/docs/library/generated/scoringrules.vrgksuv_ensemble.rst b/docs/library/generated/scoringrules.vrgksuv_ensemble.rst new file mode 100644 index 0000000..d7bae3f --- /dev/null +++ b/docs/library/generated/scoringrules.vrgksuv_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vrgksuv\_ensemble +============================== + +.. currentmodule:: scoringrules + +.. autofunction:: vrgksuv_ensemble diff --git a/docs/library/generated/scoringrules.vrvs_ensemble.rst b/docs/library/generated/scoringrules.vrvs_ensemble.rst new file mode 100644 index 0000000..a622194 --- /dev/null +++ b/docs/library/generated/scoringrules.vrvs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vrvs\_ensemble +=========================== + +.. currentmodule:: scoringrules + +.. autofunction:: vrvs_ensemble diff --git a/docs/library/generated/scoringrules.vs_ensemble.rst b/docs/library/generated/scoringrules.vs_ensemble.rst new file mode 100644 index 0000000..a37cd2e --- /dev/null +++ b/docs/library/generated/scoringrules.vs_ensemble.rst @@ -0,0 +1,6 @@ +scoringrules.vs\_ensemble +========================= + +.. currentmodule:: scoringrules + +.. autofunction:: vs_ensemble diff --git a/docs/library/generated/scoringrules.weighted_interval_score.rst b/docs/library/generated/scoringrules.weighted_interval_score.rst new file mode 100644 index 0000000..37dc886 --- /dev/null +++ b/docs/library/generated/scoringrules.weighted_interval_score.rst @@ -0,0 +1,6 @@ +scoringrules.weighted\_interval\_score +====================================== + +.. currentmodule:: scoringrules + +.. autofunction:: weighted_interval_score diff --git a/docs/reference.md b/docs/library/reference.md similarity index 82% rename from docs/reference.md rename to docs/library/reference.md index b3274fd..8cad92e 100644 --- a/docs/reference.md +++ b/docs/library/reference.md @@ -1,13 +1,33 @@ # API reference This page provides a summary of scoringrules' API. All functions are -available in the top-level namespace of the package and are here -organized by category. +available in the top-level namespace of the package and are organised here by category. ```{eval-rst} .. currentmodule:: scoringrules -Ensemble forecasts +Categorical forecasts +===================== +.. autosummary:: + :toctree: generated + + brier_score + rps_score + log_score + rls_score + + +Consistent scoring functions +============================ +.. autosummary:: + :toctree: generated + + quantile_score + interval_score + weighted_interval_score + + +Scoring rules for ensemble (sample) forecasts ================== Univariate @@ -17,14 +37,8 @@ Univariate :toctree: generated crps_ensemble - twcrps_ensemble - owcrps_ensemble - vrcrps_ensemble + logs_ensemble gksuv_ensemble - twgksuv_ensemble - owgksuv_ensemble - vrgksuv_ensemble - crps_quantile Multivariate -------------------- @@ -33,22 +47,44 @@ Multivariate :toctree: generated es_ensemble + vs_ensemble + gksmv_ensemble + +Weighted scoring rules +-------------------- + +.. autosummary:: + :toctree: generated + + twcrps_ensemble + owcrps_ensemble + vrcrps_ensemble + + clogs_ensemble + + twgksuv_ensemble + owgksuv_ensemble + vrgksuv_ensemble + owes_ensemble twes_ensemble vres_ensemble - vs_ensemble owvs_ensemble twvs_ensemble vrvs_ensemble - gksmv_ensemble twgksmv_ensemble owgksmv_ensemble vrgksmv_ensemble -Parametric distributions forecasts + +Scoring rules for parametric forecast distributions ==================================== + +CRPS +-------------------- + .. autosummary:: :toctree: generated @@ -83,6 +119,14 @@ Parametric distributions forecasts crps_quantile crps_t crps_uniform + + +Log score +-------------------- + +.. autosummary:: + :toctree: generated + logs_beta logs_binomial logs_ensemble @@ -109,25 +153,6 @@ Parametric distributions forecasts logs_tt logs_uniform -Consistent scoring functions -============================ -.. autosummary:: - :toctree: generated - - interval_score - weighted_interval_score - - -Categorical forecasts -===================== -.. autosummary:: - :toctree: generated - - brier_score - rps_score - log_score - rls_score - Backends ======== diff --git a/docs/library/user_guide.md b/docs/library/user_guide.md new file mode 100644 index 0000000..cf6e821 --- /dev/null +++ b/docs/library/user_guide.md @@ -0,0 +1,95 @@ +# User guide + +## First steps +Start by importing the library using +```python +import scoringrules as sr +``` + +The library API is simple: all metrics are available under the main namespace. Some examples are given below. + +```python +import numpy as np + +## scalar data + +# Brier score (observation = 1, forecast = 0.1) +sr.brier_score(1, 0.1) + +# CRPS (observation = 0.1, forecast = normal distribution with mean 1.2 and sd 0.3) +sr.crps_normal(0.1, 1.2, 0.3) + + +## array data + +# Brier score +obs = np.random.uniform(0, 1, 100) +fct = np.random.binomial(1, 0.5, 100) +sr.brier_score(obs, fct) + +# CRPS (forecast = lognormal distribution) +obs = np.random.lognormal(0, 1, 100) +mu = np.random.randn(100) +sig = np.random.uniform(0.5, 1.5, 100) +sr.crps_lognormal(obs, mu, sig) + + +## ensemble metrics (univariate) +obs = np.random.randn(100) +fct = obs[:, None] + np.random.randn(100, 21) * 0.1 + +sr.crps_ensemble(obs, fct) # CRPS +sr.gksuv_ensemble(obs, fct) # univariate Gaussian kernel score + + +## ensemble metrics (multivariate) +obs = np.random.randn(100, 3) +fct = obs[:, None] + np.random.randn(100, 21, 3) * 0.1 + +sr.energy_score(obs, fct) # energy score +sr.variogram_score(obs, fct) # variogram score +``` + +For the ensemble metrics, the forecast should have one more dimension than the observations, containing the ensemble (i.e. sample) members. In the univariate case, this ensemble dimension is assumed to be the last axis of the input forecast array, though this can be changed manually by specifying the `m_axis` argument (default is `m_axis=-1`). In the multivariate case, the ensemble dimension is assumed to be the second last axis, with the number of variables the last axis. These can similarly be changed manually by specifying the `m_axis` and `v_axis` arguments respectively (default is `m_axis=-2` and `v_axis=-1`). + +## Backends + +The scoringrules library supports multiple backends: `numpy` (accelerated with `numba`), `torch`, `tensorflow`, and `jax`. By default, the `numpy` and `numba` backends will be registered when importing the library. You can see the list of registered backends by running + +```python +print(sr.backends) +# {'numpy': , +# 'numba': } +``` + +and the currently active backend, used by default in all metrics, by running + +```python +print(sr.backends.active) +# +``` + +The default backend can be changed manually using `sr.backends.set_active()`. For example, the following code sets the active backend to `numba`. + +```python +sr.backends.set_active("numba") +print(sr.backends.active) +# +``` +Alternatively, the `backend` argument to the score functions can be used to override the default choice. For example, + +```python +sr.crps_normal(0.1, 1.2, 0.3, backend="numba") +``` + +To register a new backend, for example `torch`, simply use + +```python +sr.register_backend("torch") +``` + +You can now use `torch` to compute scores, either by setting it as the default backend, or by specifying it manually when running the desired function: + +```python +sr.crps_normal(0.1, 1.2, 0.3, backend="torch") +``` diff --git a/docs/user_guide.md b/docs/user_guide.md deleted file mode 100644 index b200f99..0000000 --- a/docs/user_guide.md +++ /dev/null @@ -1,75 +0,0 @@ -# User guide - -## First steps -Start by importing the library in your code with -```python -import scoringrules as sr -``` - -the library API is simple: all metrics are available under the main namespace. Let's look at some examples: - -```python -import numpy as np - -# on scalars -sr.brier_score(0.1, 1) -sr.crps_normal(0.1, 1.2, 0.3) - -# on arrays -sr.brier_score(np.random.uniform(0, 1, 100), np.random.binomial(1, 0.5, 100)) -sr.crps_lognormal(np.random.lognormal(0, 1, 100), np.random.randn(100), np.random.uniform(0.5, 1.5, 100)) - -# ensemble metrics -obs = np.random.randn(100) -fct = obs[:,None] + np.random.randn(100, 21) * 0.1 - -sr.crps_ensemble(obs, fct) -sr.error_spread_score(obs, fct) - -# multivariate ensemble metrics -obs = np.random.randn(100,3) -fct = obs[:,None] + np.random.randn(100, 21, 3) * 0.1 - -sr.energy_score(obs, fct) -sr.variogram_score(obs, fct) -``` - -For the univariate ensemble metrics, the ensemble dimension is on the last axis unless you specify otherwise with the `axis` argument. For the multivariate ensemble metrics, the ensemble dimension and the variable dimension are on the second last and last axis respectively, unless specified otherwise with `m_axis` and `v_axis`. - -## Backends -Scoringrules supports multiple backends. By default, the `numpy` and `numba` backends will be registered when importing the library. You can see the list of registered backends with - -```python -print(sr.backends) -# {'numpy': , -# 'numba': } -``` - -and the currently active backend, used by default in all metrics, can be seen with - -```python -print(sr.backends.active) -# -``` - -The default backend can also be changed with - -```python -sr.backends.set_active("numba") -print(sr.backends.active) -# -``` -When computing a metric, the `backend` argument can be used to override the default choice. - - -To register a new backend, for example `torch`, simply use - -```python -sr.register_backend("torch") -``` - -You can now use `torch` to compute metrics, either by setting it as the default backend or by specifying it on a specific metric: - -```python -sr.crps_normal(0.1, 1.0, 0.0, backend="torch") -``` diff --git a/scoringrules/_brier.py b/scoringrules/_brier.py index 554f0f2..59f884d 100644 --- a/scoringrules/_brier.py +++ b/scoringrules/_brier.py @@ -15,29 +15,41 @@ def brier_score( backend: "Backend" = None, ) -> "Array": r""" - Compute the Brier Score (BS). + Brier Score - The BS is formulated as + The Brier Score is defined as .. math:: - BS(f, y) = (f - y)^2, + \text{BS}(F, y) = (F - y)^2, - where :math:`f \in [0, 1]` is the predicted probability of an event and :math:`y \in \{0, 1\}` the actual outcome. + where :math:`F \in [0, 1]` is the predicted probability of an event and :math:`y \in \{0, 1\}` + is the outcome [1]_. Parameters ---------- obs : array_like - Observed outcome, either 0 or 1. + Observed outcomes, either 0 or 1. fct : array_like - Forecasted probabilities between 0 and 1. + Forecast probabilities, between 0 and 1. backend : str - The name of the backend used for computations. Defaults to 'numpy'. + The name of the backend used for computations. Default is 'numpy'. Returns ------- - brier_score : array_like + score : array_like The computed Brier Score. + References + ---------- + .. [1] Brier, G. W. (1950). + Verification of forecasts expressed in terms of probability. + Monthly Weather Review, 78, 1-3. + + Examples + -------- + >>> import scoringrules as sr + >>> sr.brier_score(1, 0.2) + 0.64000 """ return brier.brier_score(obs=obs, fct=fct, backend=backend) @@ -48,46 +60,76 @@ def rps_score( /, k_axis: int = -1, *, + onehot: bool = False, backend: "Backend" = None, ) -> "Array": r""" - Compute the (Discrete) Ranked Probability Score (RPS). + (Discrete) Ranked Probability Score (RPS) Suppose the outcome corresponds to one of :math:`K` ordered categories. The RPS is defined as .. math:: - RPS(f, y) = \sum_{k=1}^{K}(\tilde{f}_{k} - \tilde{y}_{k})^2, + \text{RPS}(F, y) = \sum_{k=1}^{K}(\tilde{F}_{k} - \tilde{y}_{k})^2, - where :math:`f \in [0, 1]^{K}` is a vector of length :math:`K` containing forecast probabilities - that each of the :math:`K` categories will occur, and :math:`y \in \{0, 1\}^{K}` is a vector of - length :math:`K`, with the :math:`k`-th element equal to one if the :math:`k`-th category occurs. We - have :math:`\sum_{k=1}^{K} y_{k} = \sum_{k=1}^{K} f_{k} = 1`, and, for :math:`k = 1, \dots, K`, - :math:`\tilde{y}_{k} = \sum_{i=1}^{k} y_{i}` and :math:`\tilde{f}_{k} = \sum_{i=1}^{k} f_{i}`. + where :math:`F \in [0, 1]^{K}` is a vector of length :math:`K`, containing forecast probabilities + that each of the :math:`K` categories will occur, with :math:`\sum_{k=1}^{K} F_{k} = 1` and + :math:`\tilde{F}_{k} = \sum_{i=1}^{k} F_{i}` for all :math:`k = 1, \dots, K`, and where + :math:`y \in \{1, \dots, K\}` is the category that occurs, with :math:`\tilde{y}_{k} = 1\{y \le i\}` + for all :math:`k = 1, \dots, K` [1]_. + + The outcome can alternatively be interpreted as a vector :math:`y \in \{0, 1\}^K` of length :math:`K`, with the + :math:`k`-th element equal to one if the :math:`k`-th category occurs, and zero otherwise. + Using this one-hot encoding, the RPS is defined analogously to as above, but with + :math:`\tilde{y}_{k} = \sum_{i=1}^{k} y_{i}`. Parameters ---------- obs : array_like - Array of 0's and 1's corresponding to unobserved and observed categories - forecasts : + Category that occurs. Or array of 0's and 1's corresponding to unobserved and + observed categories if `onehot=True`. + fct : array Array of forecast probabilities for each category. k_axis: int - The axis corresponding to the categories. Default is the last axis. + The axis of `obs` and `fct` corresponding to the categories. Default is the last axis. + onehot: bool + Boolean indicating whether the observation is the category that occurs or a onehot + encoded vector of 0's and 1's. Default is False. backend : str - The name of the backend used for computations. Defaults to 'numpy'. + The name of the backend used for computations. Default is 'numpy'. Returns ------- - score: + score: array_like The computed Ranked Probability Score. + References + ---------- + .. [1] Epstein, E. S. (1969). + A scoring system for probability forecasts of ranked categories. + Journal of Applied Meteorology, 8, 985-987. + https://www.jstor.org/stable/26174707. + + Examples + -------- + >>> import scoringrules as sr + >>> import numpy as np + >>> fct = np.array([0.1, 0.2, 0.3, 0.4]) + >>> obs = 3 + >>> sr.rps_score(obs, fct) + 0.25999999999999995 + >>> obs = np.array([0, 0, 1, 0]) + >>> sr.rps_score(obs, fct, onehot=True) + 0.25999999999999995 """ B = backends.active if backend is None else backends[backend] fct = B.asarray(fct) if k_axis != -1: fct = B.moveaxis(fct, k_axis, -1) + if onehot: + obs = B.moveaxis(obs, k_axis, -1) - return brier.rps_score(obs=obs, fct=fct, backend=backend) + return brier.rps_score(obs=obs, fct=fct, onehot=onehot, backend=backend) def log_score( diff --git a/scoringrules/_crps.py b/scoringrules/_crps.py index 135ffe5..1550e90 100644 --- a/scoringrules/_crps.py +++ b/scoringrules/_crps.py @@ -13,6 +13,7 @@ def crps_ensemble( /, m_axis: int = -1, *, + ens_w: "Array" = None, sorted_ensemble: bool = False, estimator: str = "pwm", backend: "Backend" = None, @@ -50,6 +51,9 @@ def crps_ensemble( represented by the last axis. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array, shape (..., m) + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. sorted_ensemble : bool Boolean indicating whether the ensemble members are already in ascending order. Default is False. @@ -102,23 +106,43 @@ def crps_ensemble( if m_axis != -1: fct = B.moveaxis(fct, m_axis, -1) - if not sorted_ensemble and estimator not in [ - "nrg", - "akr", - "akr_circperm", - "fair", - ]: + sort_ensemble = not sorted_ensemble and estimator in ["qd", "pwm", "int"] + + if sort_ensemble: + fct_uns = fct fct = B.sort(fct, axis=-1) - if backend == "numba": - if estimator not in crps.estimator_gufuncs: - raise ValueError( - f"{estimator} is not a valid estimator. " - f"Must be one of {crps.estimator_gufuncs.keys()}" - ) - return crps.estimator_gufuncs[estimator](obs, fct) + if ens_w is None: + if backend == "numba": + if estimator not in crps.estimator_gufuncs: + raise ValueError( + f"{estimator} is not a valid estimator. " + f"Must be one of {crps.estimator_gufuncs.keys()}" + ) + return crps.estimator_gufuncs[estimator](obs, fct) + + return crps.ensemble(obs, fct, estimator, backend=backend) - return crps.ensemble(obs, fct, estimator, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + if sort_ensemble: + ind = B.argsort(fct_uns, axis=-1) + ens_w = B.gather(ens_w, ind, axis=-1) + + if backend == "numba": + if estimator not in crps.estimator_gufuncs: + raise ValueError( + f"{estimator} is not a valid estimator. " + f"Must be one of {crps.estimator_gufuncs.keys()}" + ) + return crps.estimator_gufuncs_w[estimator](obs, fct, ens_w) + + return crps.ensemble_w(obs, fct, ens_w, estimator, backend=backend) def twcrps_ensemble( @@ -129,6 +153,7 @@ def twcrps_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, v_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, estimator: str = "pwm", sorted_ensemble: bool = False, @@ -162,6 +187,9 @@ def twcrps_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. v_func : callable, array_like -> array_like Chaining function used to emphasise particular outcomes. For example, a function that only considers values above a certain threshold :math:`t` by projecting forecasts and observations @@ -216,6 +244,7 @@ def v_func(x): obs, fct, m_axis=m_axis, + ens_w=ens_w, sorted_ensemble=sorted_ensemble, estimator=estimator, backend=backend, @@ -230,8 +259,8 @@ def owcrps_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, w_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, - estimator: tp.Literal["nrg"] = "nrg", backend: "Backend" = None, ) -> "Array": r"""Estimate the outcome-weighted CRPS (owCRPS) for a finite ensemble. @@ -257,7 +286,7 @@ def owcrps_ensemble( ---------- obs : array_like The observed values. - fct : array_like + fct : array The predicted forecast ensemble, where the ensemble dimension is by default represented by the last axis. a : float @@ -268,6 +297,9 @@ def owcrps_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. backend : str, optional @@ -308,11 +340,6 @@ def owcrps_ensemble( B = backends.active if backend is None else backends[backend] obs, fct = map(B.asarray, (obs, fct)) - if estimator != "nrg": - raise ValueError( - "Only the energy form of the estimator is available " - "for the outcome-weighted CRPS." - ) if m_axis != -1: fct = B.moveaxis(fct, m_axis, -1) @@ -322,14 +349,31 @@ def w_func(x): return ((a <= x) & (x <= b)) * 1.0 obs_weights, fct_weights = map(w_func, (obs, fct)) - obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights)) + if B.any(obs_weights < 0) or B.any(fct_weights < 0): + raise ValueError("`w_func` returns negative values") - if backend == "numba": - return crps.estimator_gufuncs["ow" + estimator]( - obs, fct, obs_weights, fct_weights - ) + if ens_w is None: + if backend == "numba": + return crps.estimator_gufuncs["ownrg"](obs, fct, obs_weights, fct_weights) - return crps.ow_ensemble(obs, fct, obs_weights, fct_weights, backend=backend) + return crps.ow_ensemble(obs, fct, obs_weights, fct_weights, backend=backend) + + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + + if backend == "numba": + return crps.estimator_gufuncs_w["ownrg"]( + obs, fct, obs_weights, fct_weights, ens_w + ) + + return crps.ow_ensemble_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) def vrcrps_ensemble( @@ -340,8 +384,8 @@ def vrcrps_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, w_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, - estimator: tp.Literal["nrg"] = "nrg", backend: "Backend" = None, ) -> "Array": r"""Estimate the vertically re-scaled CRPS (vrCRPS) for a finite ensemble. @@ -365,7 +409,7 @@ def vrcrps_ensemble( ---------- obs : array_like The observed values. - fct : array_like + fct : array The predicted forecast ensemble, where the ensemble dimension is by default represented by the last axis. a : float @@ -376,6 +420,9 @@ def vrcrps_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. backend : str, optional @@ -415,11 +462,6 @@ def vrcrps_ensemble( B = backends.active if backend is None else backends[backend] obs, fct = map(B.asarray, (obs, fct)) - if estimator != "nrg": - raise ValueError( - "Only the energy form of the estimator is available " - "for the outcome-weighted CRPS." - ) if m_axis != -1: fct = B.moveaxis(fct, m_axis, -1) @@ -429,21 +471,38 @@ def w_func(x): return ((a <= x) & (x <= b)) * 1.0 obs_weights, fct_weights = map(w_func, (obs, fct)) - obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights)) + if B.any(obs_weights < 0) or B.any(fct_weights < 0): + raise ValueError("`w_func` returns negative values") - if backend == "numba": - return crps.estimator_gufuncs["vr" + estimator]( - obs, fct, obs_weights, fct_weights - ) + if ens_w is None: + if backend == "numba": + return crps.estimator_gufuncs["vrnrg"](obs, fct, obs_weights, fct_weights) + + return crps.vr_ensemble(obs, fct, obs_weights, fct_weights, backend=backend) + + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + + if backend == "numba": + return crps.estimator_gufuncs_w["vrnrg"]( + obs, fct, obs_weights, fct_weights, ens_w + ) - return crps.vr_ensemble(obs, fct, obs_weights, fct_weights, backend=backend) + return crps.vr_ensemble_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) def crps_quantile( obs: "ArrayLike", fct: "Array", - alpha: "Array", /, + alpha: "Array", m_axis: int = -1, *, backend: "Backend" = None, @@ -491,11 +550,14 @@ def crps_quantile( Journal of Econometrics, 237(2), 105221. Available at https://arxiv.org/abs/2102.00968. - # TODO: add example + # TODO: add example, change reference """ B = backends.active if backend is None else backends[backend] obs, fct, alpha = map(B.asarray, (obs, fct, alpha)) + if B.any(alpha <= 0) or B.any(alpha >= 1): + raise ValueError("`alpha` contains entries that are not between 0 and 1.") + if m_axis != -1: fct = B.moveaxis(fct, m_axis, -1) @@ -510,13 +572,14 @@ def crps_quantile( def crps_beta( obs: "ArrayLike", + /, a: "ArrayLike", b: "ArrayLike", - /, lower: "ArrayLike" = 0.0, upper: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the beta distribution. @@ -550,6 +613,9 @@ def crps_beta( Upper bound of the forecast beta distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -569,16 +635,30 @@ def crps_beta( >>> sr.crps_beta(0.3, 0.7, 1.1) 0.08501024366637236 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + a, b, lower, upper = map(B.asarray, (a, b, lower, upper)) + if B.any(a <= 0): + raise ValueError( + "`a` contains non-positive entries. The shape parameters of the Beta distribution must be positive." + ) + if B.any(b <= 0): + raise ValueError( + "`b` contains non-positive entries. The shape parameters of the Beta distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.beta(obs, a, b, lower, upper, backend=backend) def crps_binomial( obs: "ArrayLike", + /, n: "ArrayLike", prob: "ArrayLike", - /, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the binomial distribution. @@ -601,6 +681,9 @@ def crps_binomial( Probability parameter of the forecast binomial distribution as a float or array of floats. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -620,15 +703,21 @@ def crps_binomial( >>> sr.crps_binomial(4, 10, 0.5) 0.5955772399902344 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + prob = B.asarray(prob) + if B.any(prob < 0) or B.any(prob > 1): + raise ValueError("`prob` contains values outside the range [0, 1].") return crps.binomial(obs, n, prob, backend=backend) def crps_exponential( obs: "ArrayLike", - rate: "ArrayLike", /, + rate: "ArrayLike", *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the exponential distribution. @@ -647,6 +736,9 @@ def crps_exponential( Rate parameter of the forecast exponential distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -669,17 +761,25 @@ def crps_exponential( >>> sr.crps_exponential(np.array([0.8, 0.9]), np.array([3.0, 2.0])) array([0.36047864, 0.31529889]) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + rate = B.asarray(rate) + if B.any(rate <= 0): + raise ValueError( + "`rate` contains non-positive entries. The rate parameter of the exponential distribution must be positive." + ) return crps.exponential(obs, rate, backend=backend) def crps_exponentialM( obs: "ArrayLike", /, - mass: "ArrayLike" = 0.0, location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, + mass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the standard exponential distribution with a point mass at the boundary. @@ -713,6 +813,9 @@ def crps_exponentialM( Scale parameter of the forecast exponential distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -732,17 +835,27 @@ def crps_exponentialM( >>> sr.crps_exponentialM(0.4, 0.2, 0.0, 1.0) 0.19251207365702294 """ - return crps.exponentialM(obs, mass, location, scale, backend=backend) + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, mass = map(B.asarray, (scale, mass)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter must be positive." + ) + if B.any(mass < 0) or B.any(mass > 1): + raise ValueError("`mass` contains entries outside the range [0, 1].") + return crps.exponentialM(obs, location, scale, mass, backend=backend) def crps_2pexponential( obs: "ArrayLike", + /, scale1: "ArrayLike", scale2: "ArrayLike", - location: "ArrayLike", - /, + location: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the two-piece exponential distribution. @@ -769,6 +882,9 @@ def crps_2pexponential( Location parameter of the forecast two-piece exponential distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -788,17 +904,29 @@ def crps_2pexponential( >>> sr.crps_2pexponential(0.8, 3.0, 1.4, 0.0) array(1.18038524) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale1, scale2 = map(B.asarray, (scale1, scale2)) + if B.any(scale1 <= 0): + raise ValueError( + "`scale1` contains non-positive entries. The scale parameters of the two-piece exponential distribution must be positive." + ) + if B.any(scale2 <= 0): + raise ValueError( + "`scale2` contains non-positive entries. The scale parameters of the two-piece exponential distribution must be positive." + ) return crps.twopexponential(obs, scale1, scale2, location, backend=backend) def crps_gamma( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", rate: "ArrayLike | None" = None, *, scale: "ArrayLike | None" = None, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the gamma distribution. @@ -827,6 +955,9 @@ def crps_gamma( Either ``rate`` or ``scale`` must be provided. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -859,17 +990,30 @@ def crps_gamma( if rate is None: rate = 1.0 / scale + if check_pars: + B = backends.active if backend is None else backends[backend] + shape, rate = map(B.asarray, (shape, rate)) + if B.any(shape <= 0): + raise ValueError( + "`shape` contains non-positive entries. The shape parameter of the gamma distribution must be positive." + ) + if B.any(rate <= 0): + raise ValueError( + "`rate` or `scale` contains non-positive entries. The rate and scale parameters of the gamma distribution must be positive." + ) + return crps.gamma(obs, shape, rate, backend=backend) def crps_gev( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the generalised extreme value (GEV) distribution. @@ -893,6 +1037,9 @@ def crps_gev( Scale parameter of the forecast GEV distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -954,18 +1101,30 @@ def crps_gev( >>> sr.crps_gev(0.3, 0.1) 0.2924712413052034 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, shape = map(B.asarray, (scale, shape)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the GEV distribution must be positive." + ) + if B.any(shape >= 1): + raise ValueError( + "`shape` contains entries larger than 1. The CRPS for the GEV distribution is only valid for shape values less than 1." + ) return crps.gev(obs, shape, location, scale, backend=backend) def crps_gpd( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, mass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the generalised pareto distribution (GPD). @@ -998,6 +1157,9 @@ def crps_gpd( Mass parameter at the lower boundary of the forecast GPD distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1017,20 +1179,36 @@ def crps_gpd( >>> sr.crps_gpd(0.3, 0.9) 0.6849331901197213 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, shape, mass = map(B.asarray, (scale, shape, mass)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the GPD distribution must be positive. `nan` is returned in these places." + ) + if B.any(shape >= 1): + raise ValueError( + "`shape` contains entries larger than 1. The CRPS for the GPD distribution is only valid for shape values less than 1. `nan` is returned in these places." + ) + if B.any(mass < 0) or B.any(mass > 1): + raise ValueError( + "`mass` contains entries outside the range [0, 1]. `nan` is returned in these places." + ) return crps.gpd(obs, shape, location, scale, mass, backend=backend) def crps_gtclogistic( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), lmass: "ArrayLike" = 0.0, umass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the generalised truncated and censored logistic distribution. @@ -1077,6 +1255,9 @@ def crps_gtclogistic( Point mass assigned to the lower boundary of the forecast distribution. umass : array_like Point mass assigned to the upper boundary of the forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1089,6 +1270,23 @@ def crps_gtclogistic( >>> sr.crps_gtclogistic(0.0, 0.1, 0.4, -1.0, 1.0, 0.1, 0.1) 0.1658713056903939 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lmass, umass, lower, upper = map( + B.asarray, (scale, lmass, umass, lower, upper) + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the generalised logistic distribution must be positive." + ) + if B.any(lmass < 0) or B.any(lmass > 1): + raise ValueError("`lmass` contains entries outside the range [0, 1].") + if B.any(umass < 0) or B.any(umass > 1): + raise ValueError("`umass` contains entries outside the range [0, 1].") + if B.any(umass + lmass >= 1): + raise ValueError("The sum of `umass` and `lmass` should be smaller than 1.") + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtclogistic( obs, location, @@ -1103,13 +1301,14 @@ def crps_gtclogistic( def crps_tlogistic( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the truncated logistic distribution. @@ -1128,6 +1327,9 @@ def crps_tlogistic( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1140,20 +1342,37 @@ def crps_tlogistic( >>> sr.crps_tlogistic(0.0, 0.1, 0.4, -1.0, 1.0) 0.12714830546327846 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated logistic distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtclogistic( - obs, location, scale, lower, upper, 0.0, 0.0, backend=backend + obs, + location, + scale, + lower, + upper, + 0.0, + 0.0, + backend=backend, ) def crps_clogistic( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the censored logistic distribution. @@ -1172,6 +1391,9 @@ def crps_clogistic( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1184,6 +1406,15 @@ def crps_clogistic( >>> sr.crps_clogistic(0.0, 0.1, 0.4, -1.0, 1.0) 0.15805632276434345 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the censored logistic distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") lmass = stats._logis_cdf((lower - location) / scale) umass = 1 - stats._logis_cdf((upper - location) / scale) return crps.gtclogistic( @@ -1200,15 +1431,16 @@ def crps_clogistic( def crps_gtcnormal( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), lmass: "ArrayLike" = 0.0, umass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the generalised truncated and censored normal distribution. @@ -1242,6 +1474,23 @@ def crps_gtcnormal( >>> sr.crps_gtcnormal(0.0, 0.1, 0.4, -1.0, 1.0, 0.1, 0.1) 0.1351100832878575 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lmass, umass, lower, upper = map( + B.asarray, (scale, lmass, umass, lower, upper) + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the generalised normal distribution must be positive." + ) + if B.any(lmass < 0) or B.any(lmass > 1): + raise ValueError("`lmass` contains entries outside the range [0, 1].") + if B.any(umass < 0) or B.any(umass > 1): + raise ValueError("`umass` contains entries outside the range [0, 1].") + if B.any(umass + lmass >= 1): + raise ValueError("The sum of `umass` and `lmass` should be smaller than 1.") + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtcnormal( obs, location, @@ -1256,13 +1505,14 @@ def crps_gtcnormal( def crps_tnormal( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the truncated normal distribution. @@ -1281,6 +1531,9 @@ def crps_tnormal( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1293,18 +1546,28 @@ def crps_tnormal( >>> sr.crps_tnormal(0.0, 0.1, 0.4, -1.0, 1.0) 0.10070146718008832 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated normal distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtcnormal(obs, location, scale, lower, upper, 0.0, 0.0, backend=backend) def crps_cnormal( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the censored normal distribution. @@ -1323,6 +1586,9 @@ def crps_cnormal( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1335,6 +1601,15 @@ def crps_cnormal( >>> sr.crps_cnormal(0.0, 0.1, 0.4, -1.0, 1.0) 0.10338851213123085 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the censored normal distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") lmass = stats._norm_cdf((lower - location) / scale) umass = 1 - stats._norm_cdf((upper - location) / scale) return crps.gtcnormal( @@ -1351,8 +1626,8 @@ def crps_cnormal( def crps_gtct( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), @@ -1361,6 +1636,7 @@ def crps_gtct( umass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the generalised truncated and censored t distribution. @@ -1405,6 +1681,9 @@ def crps_gtct( Point mass assigned to the lower boundary of the forecast distribution. umass : array_like Point mass assigned to the upper boundary of the forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1424,6 +1703,27 @@ def crps_gtct( >>> sr.crps_gtct(0.0, 2.0, 0.1, 0.4, -1.0, 1.0, 0.1, 0.1) 0.13997789333289662 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lmass, umass, lower, upper = map( + B.asarray, (scale, lmass, umass, lower, upper) + ) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the generalised t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the generalised t distribution must be positive." + ) + if B.any(lmass < 0) or B.any(lmass > 1): + raise ValueError("`lmass` contains entries outside the range [0, 1].") + if B.any(umass < 0) or B.any(umass > 1): + raise ValueError("`umass` contains entries outside the range [0, 1].") + if B.any(umass + lmass >= 1): + raise ValueError("The sum of `umass` and `lmass` should be smaller than 1.") + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtct( obs, df, @@ -1439,14 +1739,15 @@ def crps_gtct( def crps_tt( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the truncated t distribution. @@ -1467,6 +1768,9 @@ def crps_tt( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1479,6 +1783,19 @@ def crps_tt( >>> sr.crps_tt(0.0, 2.0, 0.1, 0.4, -1.0, 1.0) 0.10323007471747117 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the truncated t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated t distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return crps.gtct( obs, df, @@ -1494,14 +1811,15 @@ def crps_tt( def crps_ct( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the censored t distribution. @@ -1522,6 +1840,9 @@ def crps_ct( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1534,6 +1855,19 @@ def crps_ct( >>> sr.crps_ct(0.0, 2.0, 0.1, 0.4, -1.0, 1.0) 0.12672580744453948 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the censored t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the censored t distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") lmass = stats._t_cdf((lower - location) / scale, df) umass = 1 - stats._t_cdf((upper - location) / scale, df) return crps.gtct( @@ -1551,12 +1885,13 @@ def crps_ct( def crps_hypergeometric( obs: "ArrayLike", + /, m: "ArrayLike", n: "ArrayLike", k: "ArrayLike", - /, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the hypergeometric distribution. @@ -1582,6 +1917,9 @@ def crps_hypergeometric( Number of draws, without replacement. Must be in 0, 1, ..., m + n. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1601,6 +1939,7 @@ def crps_hypergeometric( >>> sr.crps_hypergeometric(5, 7, 13, 12) 0.44697415547610597 """ + # TODO: add check that m,n,k are integers return crps.hypergeometric(obs, m, n, k, backend=backend) @@ -1611,6 +1950,7 @@ def crps_laplace( scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the laplace distribution. @@ -1633,6 +1973,9 @@ def crps_laplace( Scale parameter of the forecast laplace distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1652,16 +1995,24 @@ def crps_laplace( >>> sr.crps_laplace(0.3, 0.1, 0.2) 0.12357588823428847 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale = B.asarray(scale) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the laplace must be positive." + ) return crps.laplace(obs, location, scale, backend=backend) def crps_logistic( obs: "ArrayLike", - mu: "ArrayLike", - sigma: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the logistic distribution. @@ -1677,15 +2028,18 @@ def crps_logistic( ---------- obs : array_like Observed values. - mu: array_like + location: array_like Location parameter of the forecast logistic distribution. - sigma: array_like + scale: array_like Scale parameter of the forecast logistic distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- crps : array_like - The CRPS for the Logistic(mu, sigma) forecasts given the observations. + The CRPS for the Logistic(location, scale) forecasts given the observations. References ---------- @@ -1700,15 +2054,24 @@ def crps_logistic( >>> sr.crps_logistic(0.0, 0.4, 0.1) 0.3036299855835619 """ - return crps.logistic(obs, mu, sigma, backend=backend) + if check_pars: + B = backends.active if backend is None else backends[backend] + scale = B.asarray(scale) + if B.any(scale <= 0): + raise ValueError( + "`sigma` contains non-positive entries. The scale parameter of the logistic distribution must be positive." + ) + return crps.logistic(obs, location, scale, backend=backend) def crps_loglaplace( obs: "ArrayLike", - locationlog: "ArrayLike", - scalelog: "ArrayLike", + /, + locationlog: "ArrayLike" = 0.0, + scalelog: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the log-Laplace distribution. @@ -1741,6 +2104,9 @@ def crps_loglaplace( Scale parameter of the forecast log-laplace distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1760,14 +2126,24 @@ def crps_loglaplace( >>> sr.crps_loglaplace(3.0, 0.1, 0.9) 1.162020513653791 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scalelog = B.asarray(scalelog) + if B.any(scalelog <= 0) or B.any(scalelog >= 1): + raise ValueError( + "`scalelog` contains entries outside of the range (0, 1). The scale parameter of the log-laplace distribution must be between 0 and 1." + ) return crps.loglaplace(obs, locationlog, scalelog, backend=backend) def crps_loglogistic( obs: "ArrayLike", - mulog: "ArrayLike", - sigmalog: "ArrayLike", + /, + mulog: "ArrayLike" = 0.0, + sigmalog: "ArrayLike" = 1.0, + *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the log-logistic distribution. @@ -1800,7 +2176,9 @@ def crps_loglogistic( Scale parameter of the log-logistic distribution. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. - + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1820,14 +2198,24 @@ def crps_loglogistic( >>> sr.crps_loglogistic(3.0, 0.1, 0.9) 1.1329527730161177 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigmalog = B.asarray(sigmalog) + if B.any(sigmalog <= 0) or B.any(sigmalog >= 1): + raise ValueError( + "`sigmalog` contains entries outside of the range (0, 1). The scale parameter of the log-logistic distribution must be between 0 and 1." + ) return crps.loglogistic(obs, mulog, sigmalog, backend=backend) def crps_lognormal( obs: "ArrayLike", - mulog: "ArrayLike", - sigmalog: "ArrayLike", + /, + mulog: "ArrayLike" = 0.0, + sigmalog: "ArrayLike" = 1.0, + *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the lognormal distribution. @@ -1852,6 +2240,9 @@ def crps_lognormal( Mean of the normal underlying distribution. sigmalog : array_like Standard deviation of the underlying normal distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1870,18 +2261,26 @@ def crps_lognormal( >>> sr.crps_lognormal(0.1, 0.4, 0.0) 1.3918246976412703 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigmalog = B.asarray(sigmalog) + if B.any(sigmalog <= 0): + raise ValueError( + "`sigmalog` contains non-positive entries. The scale parameter of the log-normal distribution must be positive." + ) return crps.lognormal(obs, mulog, sigmalog, backend=backend) def crps_mixnorm( obs: "ArrayLike", - m: "ArrayLike", - s: "ArrayLike", /, + m: "ArrayLike" = 0.0, + s: "ArrayLike" = 1.0, w: "ArrayLike" = None, m_axis: "ArrayLike" = -1, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for a mixture of normal distributions. @@ -1907,6 +2306,9 @@ def crps_mixnorm( The axis corresponding to the mixture components. Default is the last axis. backend : str, optional The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1934,6 +2336,15 @@ def crps_mixnorm( w = B.zeros(m.shape) + 1 / M else: w = B.asarray(w) + w = w / B.sum(w, axis=m_axis, keepdims=True) + + if check_pars: + if B.any(s <= 0): + raise ValueError( + "`s` contains non-positive entries. The scale parameters of the normal distributions should be positive." + ) + if B.any(w < 0): + raise ValueError("`w` contains negative entries") if m_axis != -1: m = B.moveaxis(m, m_axis, -1) @@ -1945,12 +2356,13 @@ def crps_mixnorm( def crps_negbinom( obs: "ArrayLike", - n: "ArrayLike", /, + n: "ArrayLike", prob: "ArrayLike | None" = None, *, mu: "ArrayLike | None" = None, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the negative binomial distribution. @@ -1973,6 +2385,9 @@ def crps_negbinom( Probability parameter of the forecast negative binomial distribution. mu: array_like Mean of the forecast negative binomial distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2004,16 +2419,23 @@ def crps_negbinom( if prob is None: prob = n / (n + mu) + if check_pars: + B = backends.active if backend is None else backends[backend] + prob = B.asarray(prob) + if B.any(prob < 0) or B.any(prob > 1): + raise ValueError("`prob` contains values outside the range [0, 1].") + return crps.negbinom(obs, n, prob, backend=backend) def crps_normal( obs: "ArrayLike", - mu: "ArrayLike", - sigma: "ArrayLike", /, + mu: "ArrayLike" = 0.0, + sigma: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the normal distribution. @@ -2033,6 +2455,9 @@ def crps_normal( Mean of the forecast normal distribution. sigma: array_like Standard deviation of the forecast normal distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2051,17 +2476,25 @@ def crps_normal( >>> sr.crps_normal(0.0, 0.1, 0.4) 0.10339992515976162 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigma = B.asarray(sigma) + if B.any(sigma <= 0): + raise ValueError( + "`sigma` contains non-positive entries. The standard deviation of the normal distribution must be positive." + ) return crps.normal(obs, mu, sigma, backend=backend) def crps_2pnormal( obs: "ArrayLike", - scale1: "ArrayLike", - scale2: "ArrayLike", - location: "ArrayLike", /, + scale1: "ArrayLike" = 1.0, + scale2: "ArrayLike" = 1.0, + location: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the two-piece normal distribution. @@ -2086,6 +2519,9 @@ def crps_2pnormal( Scale parameter of the upper half of the forecast two-piece normal distribution. mu: array_like Location parameter of the forecast two-piece normal distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2107,6 +2543,15 @@ def crps_2pnormal( """ B = backends.active if backend is None else backends[backend] obs, scale1, scale2, location = map(B.asarray, (obs, scale1, scale2, location)) + if check_pars: + if B.any(scale1 <= 0): + raise ValueError( + "`scale1` contains non-positive entries. The scale parameters of the two-piece normal distribution must be positive." + ) + if B.any(scale2 <= 0): + raise ValueError( + "`scale2` contains non-positive entries. The scale parameters of the two-piece normal distribution must be positive." + ) lower = float("-inf") upper = 0.0 lmass = 0.0 @@ -2128,10 +2573,11 @@ def crps_2pnormal( def crps_poisson( obs: "ArrayLike", - mean: "ArrayLike", /, + mean: "ArrayLike", *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the Poisson distribution. @@ -2151,6 +2597,9 @@ def crps_poisson( The observed values. mean : array_like Mean parameter of the forecast poisson distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2169,17 +2618,25 @@ def crps_poisson( >>> sr.crps_poisson(1, 2) 0.4991650450203817 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + mean = B.asarray(mean) + if B.any(mean <= 0): + raise ValueError( + "`mean` contains non-positive entries. The mean parameter of the Poisson distribution must be positive." + ) return crps.poisson(obs, mean, backend=backend) def crps_t( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the student's t distribution. @@ -2204,8 +2661,11 @@ def crps_t( Degrees of freedom parameter of the forecast t distribution. location : array_like Location parameter of the forecast t distribution. - sigma : array_like + scale : array_like Scale parameter of the forecast t distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2225,18 +2685,30 @@ def crps_t( >>> sr.crps_t(0.0, 0.1, 0.4, 0.1) 0.07687151141732129 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + df, scale = map(B.asarray, (df, scale)) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the t distribution must be positive." + ) return crps.t(obs, df, location, scale, backend=backend) def crps_uniform( obs: "ArrayLike", - min: "ArrayLike", - max: "ArrayLike", /, + min: "ArrayLike" = 0.0, + max: "ArrayLike" = 1.0, lmass: "ArrayLike" = 0.0, umass: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the closed form of the CRPS for the uniform distribution. @@ -2264,6 +2736,9 @@ def crps_uniform( Point mass on the lower bound of the forecast uniform distribution. umass : array_like Point mass on the upper bound of the forecast uniform distribution. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -2283,6 +2758,17 @@ def crps_uniform( >>> sr.crps_uniform(0.4, 0.0, 1.0, 0.0, 0.0) 0.09333333333333332 """ + if check_pars: + B = backends.active if backend is None else backends[backend] + lmass, umass, min, max = map(B.asarray, (lmass, umass, min, max)) + if B.any(lmass < 0) or B.any(lmass > 1): + raise ValueError("`lmass` contains entries outside the range [0, 1].") + if B.any(umass < 0) or B.any(umass > 1): + raise ValueError("`umass` contains entries outside the range [0, 1].") + if B.any(umass + lmass >= 1): + raise ValueError("The sum of `umass` and `lmass` should be smaller than 1.") + if B.any(min >= max): + raise ValueError("`min` is not always smaller than `max`.") return crps.uniform(obs, min, max, lmass, umass, backend=backend) diff --git a/scoringrules/_energy.py b/scoringrules/_energy.py index b4fd990..7732099 100644 --- a/scoringrules/_energy.py +++ b/scoringrules/_energy.py @@ -15,6 +15,7 @@ def es_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the Energy Score for a finite multivariate ensemble. @@ -39,6 +40,9 @@ def es_ensemble( v_axis : int The axis corresponding to the variables dimension on the forecasts array (or the observations array with an extra dimension on `m_axis`). Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -59,13 +63,23 @@ def es_ensemble( :ref:`theory.multivariate` Some theoretical background on scoring rules for multivariate forecasts. """ - backend = backend if backend is not None else backends._active + B = backends.active if backend is None else backends[backend] obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) - if backend == "numba": - return energy._energy_score_gufunc(obs, fct) + if ens_w is None: + if backend == "numba": + return energy._energy_score_gufunc(obs, fct) + + return energy.nrg(obs, fct, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if backend == "numba": + return energy._energy_score_gufunc_w(obs, fct, ens_w) - return energy.nrg(obs, fct, backend=backend) + return energy.nrg_w(obs, fct, ens_w, backend=backend) def twes_ensemble( @@ -76,6 +90,7 @@ def twes_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the Threshold-Weighted Energy Score (twES) for a finite multivariate ensemble. @@ -105,6 +120,9 @@ def twes_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -114,7 +132,9 @@ def twes_ensemble( The computed Threshold-Weighted Energy Score. """ obs, fct = map(v_func, (obs, fct)) - return es_ensemble(obs, fct, m_axis=m_axis, v_axis=v_axis, backend=backend) + return es_ensemble( + obs, fct, m_axis=m_axis, v_axis=v_axis, ens_w=ens_w, backend=backend + ) def owes_ensemble( @@ -125,6 +145,7 @@ def owes_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the Outcome-Weighted Energy Score (owES) for a finite multivariate ensemble. @@ -157,6 +178,9 @@ def owes_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of ints The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -166,16 +190,29 @@ def owes_ensemble( The computed Outcome-Weighted Energy Score. """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) fct_weights = B.apply_along_axis(w_func, fct, -1) obs_weights = B.apply_along_axis(w_func, obs, -1) - if B.name == "numba": - return energy._owenergy_score_gufunc(obs, fct, obs_weights, fct_weights) + if ens_w is None: + if B.name == "numba": + return energy._owenergy_score_gufunc(obs, fct, obs_weights, fct_weights) + + return energy.ownrg(obs, fct, obs_weights, fct_weights, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if B.name == "numba": + return energy._owenergy_score_gufunc_w( + obs, fct, obs_weights, fct_weights, ens_w + ) - return energy.ownrg(obs, fct, obs_weights, fct_weights, backend=backend) + return energy.ownrg_w( + obs, fct, obs_weights, fct_weights, ens_w=ens_w, backend=backend + ) def vres_ensemble( @@ -183,9 +220,10 @@ def vres_ensemble( fct: "Array", w_func: tp.Callable[["ArrayLike"], "ArrayLike"], /, - *, m_axis: int = -2, v_axis: int = -1, + *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the Vertically Re-scaled Energy Score (vrES) for a finite multivariate ensemble. @@ -219,6 +257,9 @@ def vres_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -228,13 +269,26 @@ def vres_ensemble( The computed Vertically Re-scaled Energy Score. """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) fct_weights = B.apply_along_axis(w_func, fct, -1) obs_weights = B.apply_along_axis(w_func, obs, -1) - if backend == "numba": - return energy._vrenergy_score_gufunc(obs, fct, obs_weights, fct_weights) - - return energy.vrnrg(obs, fct, obs_weights, fct_weights, backend=backend) + if ens_w is None: + if backend == "numba": + return energy._vrenergy_score_gufunc(obs, fct, obs_weights, fct_weights) + + return energy.vrnrg(obs, fct, obs_weights, fct_weights, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if backend == "numba": + return energy._vrenergy_score_gufunc_w( + obs, fct, obs_weights, fct_weights, ens_w + ) + + return energy.vrnrg_w( + obs, fct, obs_weights, fct_weights, ens_w=ens_w, backend=backend + ) diff --git a/scoringrules/_kernels.py b/scoringrules/_kernels.py index 3ea14c2..3596b7b 100644 --- a/scoringrules/_kernels.py +++ b/scoringrules/_kernels.py @@ -14,6 +14,7 @@ def gksuv_ensemble( /, m_axis: int = -1, *, + ens_w: "Array" = None, estimator: str = "nrg", backend: "Backend" = None, ) -> "Array": @@ -37,11 +38,14 @@ def gksuv_ensemble( ---------- obs : array_like The observed values. - fct : array_like + fct : array The predicted forecast ensemble, where the ensemble dimension is by default represented by the last axis. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. estimator : str Indicates the estimator to be used. backend : str @@ -60,11 +64,14 @@ def gksuv_ensemble( B = backends.active if backend is None else backends[backend] obs, fct = map(B.asarray, (obs, fct)) + if m_axis != -1: + fct = B.moveaxis(fct, m_axis, -1) + if backend == "numba": - if estimator not in kernels.estimator_gufuncs: + if estimator not in kernels.estimator_gufuncs_uv: raise ValueError( f"{estimator} is not a valid estimator. " - f"Must be one of {kernels.estimator_gufuncs.keys()}" + f"Must be one of {kernels.estimator_gufuncs_uv.keys()}" ) else: if estimator not in ["fair", "nrg"]: @@ -73,13 +80,26 @@ def gksuv_ensemble( f"Must be one of ['fair', 'nrg']" ) - if m_axis != -1: - fct = B.moveaxis(fct, m_axis, -1) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_uv[estimator](obs, fct) - if backend == "numba": - return kernels.estimator_gufuncs[estimator](obs, fct) + return kernels.ensemble_uv(obs, fct, estimator=estimator, backend=backend) - return kernels.ensemble_uv(obs, fct, estimator, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + + if backend == "numba": + return kernels.estimator_gufuncs_uv_w[estimator](obs, fct, ens_w) + + return kernels.ensemble_uv_w( + obs, fct, ens_w=ens_w, estimator=estimator, backend=backend + ) def twgksuv_ensemble( @@ -90,6 +110,7 @@ def twgksuv_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, v_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, estimator: str = "nrg", backend: "Backend" = None, @@ -126,6 +147,9 @@ def twgksuv_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. v_func : callable, array_like -> array_like Chaining function used to emphasise particular outcomes. For example, a function that only considers values above a certain threshold :math:`t` by projecting forecasts and observations @@ -160,6 +184,7 @@ def v_func(x): obs, fct, m_axis=m_axis, + ens_w=ens_w, estimator=estimator, backend=backend, ) @@ -173,6 +198,7 @@ def owgksuv_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, w_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, backend: "Backend" = None, ) -> "Array": @@ -209,6 +235,9 @@ def owgksuv_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. backend : str @@ -242,11 +271,35 @@ def w_func(x): obs_weights, fct_weights = map(w_func, (obs, fct)) obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights)) + if B.any(obs_weights < 0) or B.any(fct_weights < 0): + raise ValueError("`w_func` returns negative values") - if backend == "numba": - return kernels.estimator_gufuncs["ow"](obs, fct, obs_weights, fct_weights) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_uv["ow"]( + obs, fct, obs_weights, fct_weights + ) + + return kernels.ow_ensemble_uv( + obs, fct, obs_weights, fct_weights, backend=backend + ) - return kernels.ow_ensemble_uv(obs, fct, obs_weights, fct_weights, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + + if backend == "numba": + return kernels.estimator_gufuncs_uv_w["ow"]( + obs, fct, obs_weights, fct_weights, ens_w + ) + + return kernels.ow_ensemble_uv_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) def vrgksuv_ensemble( @@ -257,6 +310,7 @@ def vrgksuv_ensemble( b: float = float("inf"), m_axis: int = -1, *, + ens_w: "Array" = None, w_func: tp.Callable[["ArrayLike"], "ArrayLike"] = None, backend: "Backend" = None, ) -> "Array": @@ -292,6 +346,9 @@ def vrgksuv_ensemble( to values in the range [a, b]. m_axis : int The axis corresponding to the ensemble. Default is the last axis. + ens_w : array + Weights assigned to the ensemble members. Array with the same shape as fct. + Default is equal weighting. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. backend : str @@ -325,11 +382,35 @@ def w_func(x): obs_weights, fct_weights = map(w_func, (obs, fct)) obs_weights, fct_weights = map(B.asarray, (obs_weights, fct_weights)) + if B.any(obs_weights < 0) or B.any(fct_weights < 0): + raise ValueError("`w_func` returns negative values") - if backend == "numba": - return kernels.estimator_gufuncs["vr"](obs, fct, obs_weights, fct_weights) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_uv["vr"]( + obs, fct, obs_weights, fct_weights + ) - return kernels.vr_ensemble_uv(obs, fct, obs_weights, fct_weights, backend=backend) + return kernels.vr_ensemble_uv( + obs, fct, obs_weights, fct_weights, backend=backend + ) + + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=m_axis, keepdims=True) + if m_axis != -1: + ens_w = B.moveaxis(ens_w, m_axis, -1) + + if backend == "numba": + return kernels.estimator_gufuncs_uv_w["vr"]( + obs, fct, obs_weights, fct_weights, ens_w + ) + + return kernels.vr_ensemble_uv_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) def gksmv_ensemble( @@ -339,6 +420,7 @@ def gksmv_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, estimator: str = "nrg", backend: "Backend" = None, ) -> "Array": @@ -373,6 +455,9 @@ def gksmv_ensemble( v_axis : int The axis corresponding to the variables dimension on the forecasts array (or the observations array with an extra dimension on `m_axis`). Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. estimator : str Indicates the estimator to be used. backend : str @@ -383,7 +468,7 @@ def gksmv_ensemble( score : array_like The GKS between the forecast ensemble and obs. """ - backend = backend if backend is not None else backends._active + B = backends.active if backend is None else backends[backend] obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) if estimator not in kernels.estimator_gufuncs_mv: @@ -392,10 +477,22 @@ def gksmv_ensemble( f"Must be one of {kernels.estimator_gufuncs_mv.keys()}" ) - if backend == "numba": - return kernels.estimator_gufuncs_mv[estimator](obs, fct) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_mv[estimator](obs, fct) - return kernels.ensemble_mv(obs, fct, estimator, backend=backend) + return kernels.ensemble_mv(obs, fct, estimator=estimator, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if backend == "numba": + return kernels.estimator_gufuncs_mv_w[estimator](obs, fct, ens_w) + + return kernels.ensemble_mv_w( + obs, fct, ens_w, estimator=estimator, backend=backend + ) def twgksmv_ensemble( @@ -406,6 +503,8 @@ def twgksmv_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, + estimator: str = "nrg", backend: "Backend" = None, ) -> "Array": r"""Compute the Threshold-Weighted Gaussian Kernel Score (twGKS) for a finite multivariate ensemble. @@ -438,6 +537,9 @@ def twgksmv_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of ints The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -447,7 +549,15 @@ def twgksmv_ensemble( The computed Threshold-Weighted Gaussian Kernel Score. """ obs, fct = map(v_func, (obs, fct)) - return gksmv_ensemble(obs, fct, m_axis=m_axis, v_axis=v_axis, backend=backend) + return gksmv_ensemble( + obs, + fct, + m_axis=m_axis, + v_axis=v_axis, + ens_w=ens_w, + estimator=estimator, + backend=backend, + ) def owgksmv_ensemble( @@ -458,6 +568,7 @@ def owgksmv_ensemble( m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the multivariate Outcome-Weighted Gaussian Kernel Score (owGKS) for a finite ensemble. @@ -505,6 +616,9 @@ def owgksmv_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of ints The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -514,16 +628,33 @@ def owgksmv_ensemble( The computed Outcome-Weighted GKS. """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) fct_weights = B.apply_along_axis(w_func, fct, -1) obs_weights = B.apply_along_axis(w_func, obs, -1) - if B.name == "numba": - return kernels.estimator_gufuncs_mv["ow"](obs, fct, obs_weights, fct_weights) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_mv["ow"]( + obs, fct, obs_weights, fct_weights + ) + + return kernels.ow_ensemble_mv( + obs, fct, obs_weights, fct_weights, backend=backend + ) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if backend == "numba": + return kernels.estimator_gufuncs_mv_w["ow"]( + obs, fct, obs_weights, fct_weights, ens_w + ) - return kernels.ow_ensemble_mv(obs, fct, obs_weights, fct_weights, backend=backend) + return kernels.ow_ensemble_mv_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) def vrgksmv_ensemble( @@ -531,9 +662,10 @@ def vrgksmv_ensemble( fct: "Array", w_func: tp.Callable[["ArrayLike"], "ArrayLike"], /, - *, m_axis: int = -2, v_axis: int = -1, + *, + ens_w: "Array" = None, backend: "Backend" = None, ) -> "Array": r"""Compute the Vertically Re-scaled Gaussian Kernel Score (vrGKS) for a finite multivariate ensemble. @@ -567,6 +699,9 @@ def vrgksmv_ensemble( The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int or tuple of ints The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. backend: str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -576,16 +711,33 @@ def vrgksmv_ensemble( The computed Vertically Re-scaled Gaussian Kernel Score. """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) fct_weights = B.apply_along_axis(w_func, fct, -1) obs_weights = B.apply_along_axis(w_func, obs, -1) - if B.name == "numba": - return kernels.estimator_gufuncs_mv["vr"](obs, fct, obs_weights, fct_weights) + if ens_w is None: + if backend == "numba": + return kernels.estimator_gufuncs_mv["vr"]( + obs, fct, obs_weights, fct_weights + ) + + return kernels.vr_ensemble_mv( + obs, fct, obs_weights, fct_weights, backend=backend + ) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + if backend == "numba": + return kernels.estimator_gufuncs_mv_w["vr"]( + obs, fct, obs_weights, fct_weights, ens_w + ) - return kernels.vr_ensemble_mv(obs, fct, obs_weights, fct_weights, backend=backend) + return kernels.vr_ensemble_mv_w( + obs, fct, obs_weights, fct_weights, ens_w, backend=backend + ) __all__ = [ diff --git a/scoringrules/_logs.py b/scoringrules/_logs.py index 6b5f96d..8ac32aa 100644 --- a/scoringrules/_logs.py +++ b/scoringrules/_logs.py @@ -156,13 +156,14 @@ def clogs_ensemble( def logs_beta( obs: "ArrayLike", + /, a: "ArrayLike", b: "ArrayLike", - /, lower: "ArrayLike" = 0.0, upper: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the beta distribution. @@ -181,7 +182,10 @@ def logs_beta( upper : array_like Upper bound of the forecast beta distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -193,16 +197,30 @@ def logs_beta( >>> import scoringrules as sr >>> sr.logs_beta(0.3, 0.7, 1.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + a, b, lower, upper = map(B.asarray, (a, b, lower, upper)) + if B.any(a <= 0): + raise ValueError( + "`a` contains non-positive entries. The shape parameters of the Beta distribution must be positive." + ) + if B.any(b <= 0): + raise ValueError( + "`b` contains non-positive entries. The shape parameters of the Beta distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return logarithmic.beta(obs, a, b, lower, upper, backend=backend) def logs_binomial( obs: "ArrayLike", + /, n: "ArrayLike", prob: "ArrayLike", - /, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the binomial distribution. @@ -217,7 +235,10 @@ def logs_binomial( prob : array_like Probability parameter of the forecast binomial distribution as a float or array of floats. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -229,15 +250,21 @@ def logs_binomial( >>> import scoringrules as sr >>> sr.logs_binomial(4, 10, 0.5) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + prob = B.asarray(prob) + if B.any(prob < 0) or B.any(prob > 1): + raise ValueError("`prob` contains values outside the range [0, 1].") return logarithmic.binomial(obs, n, prob, backend=backend) def logs_exponential( obs: "ArrayLike", - rate: "ArrayLike", /, + rate: "ArrayLike", *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the exponential distribution. @@ -250,7 +277,10 @@ def logs_exponential( rate : array_like Rate parameter of the forecast exponential distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -262,6 +292,13 @@ def logs_exponential( >>> import scoringrules as sr >>> sr.logs_exponential(0.8, 3.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + rate = B.asarray(rate) + if B.any(rate <= 0): + raise ValueError( + "`rate` contains non-positive entries. The rate parameter of the exponential distribution must be positive." + ) return logarithmic.exponential(obs, rate, backend=backend) @@ -272,6 +309,7 @@ def logs_exponential2( scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the exponential distribution with location and scale parameters. @@ -286,7 +324,10 @@ def logs_exponential2( scale : array_like Scale parameter of the forecast exponential distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -298,17 +339,25 @@ def logs_exponential2( >>> import scoringrules as sr >>> sr.logs_exponential2(0.2, 0.0, 1.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale = B.asarray(scale) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the exponential distribution must be positive." + ) return logarithmic.exponential2(obs, location, scale, backend=backend) def logs_2pexponential( obs: "ArrayLike", + /, scale1: "ArrayLike", scale2: "ArrayLike", - location: "ArrayLike", - /, + location: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the two-piece exponential distribution. @@ -325,7 +374,10 @@ def logs_2pexponential( location : array_like Location parameter of the forecast two-piece exponential distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -337,17 +389,29 @@ def logs_2pexponential( >>> import scoringrules as sr >>> sr.logs_2pexponential(0.8, 3.0, 1.4, 0.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale1, scale2 = map(B.asarray, (scale1, scale2)) + if B.any(scale1 <= 0): + raise ValueError( + "`scale1` contains non-positive entries. The scale parameters of the two-piece exponential distribution must be positive." + ) + if B.any(scale2 <= 0): + raise ValueError( + "`scale2` contains non-positive entries. The scale parameters of the two-piece exponential distribution must be positive." + ) return logarithmic.twopexponential(obs, scale1, scale2, location, backend=backend) def logs_gamma( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", rate: "ArrayLike | None" = None, *, scale: "ArrayLike | None" = None, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the gamma distribution. @@ -363,6 +427,11 @@ def logs_gamma( Rate parameter of the forecast gamma distribution. scale : array_like Scale parameter of the forecast gamma distribution, where `scale = 1 / rate`. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -387,17 +456,30 @@ def logs_gamma( if rate is None: rate = 1.0 / scale + if check_pars: + B = backends.active if backend is None else backends[backend] + shape, rate = map(B.asarray, (shape, rate)) + if B.any(shape <= 0): + raise ValueError( + "`shape` contains non-positive entries. The shape parameter of the gamma distribution must be positive." + ) + if B.any(rate <= 0): + raise ValueError( + "`rate` or `scale` contains non-positive entries. The rate and scale parameters of the gamma distribution must be positive." + ) + return logarithmic.gamma(obs, shape, rate, backend=backend) def logs_gev( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the generalised extreme value (GEV) distribution. @@ -414,7 +496,10 @@ def logs_gev( scale : array_like Scale parameter of the forecast GEV distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -426,17 +511,25 @@ def logs_gev( >>> import scoringrules as sr >>> sr.logs_gev(0.3, 0.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, shape = map(B.asarray, (scale, shape)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the GEV distribution must be positive." + ) return logarithmic.gev(obs, shape, location, scale, backend=backend) def logs_gpd( obs: "ArrayLike", - shape: "ArrayLike", /, + shape: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the generalised Pareto distribution (GPD). @@ -454,7 +547,10 @@ def logs_gpd( scale : array_like Scale parameter of the forecast GPD distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -466,17 +562,25 @@ def logs_gpd( >>> import scoringrules as sr >>> sr.logs_gpd(0.3, 0.9) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, shape = map(B.asarray, (scale, shape)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the GPD distribution must be positive. `nan` is returned in these places." + ) return logarithmic.gpd(obs, shape, location, scale, backend=backend) def logs_hypergeometric( obs: "ArrayLike", + /, m: "ArrayLike", n: "ArrayLike", k: "ArrayLike", - /, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the hypergeometric distribution. @@ -493,7 +597,10 @@ def logs_hypergeometric( k : array_like Number of draws, without replacement. Must be in 0, 1, ..., m + n. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -505,16 +612,18 @@ def logs_hypergeometric( >>> import scoringrules as sr >>> sr.logs_hypergeometric(5, 7, 13, 12) """ + # TODO: add check that m,n,k are integers return logarithmic.hypergeometric(obs, m, n, k, backend=backend) def logs_laplace( obs: "ArrayLike", + /, location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, - /, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the Laplace distribution. @@ -530,6 +639,11 @@ def logs_laplace( scale : array_like Scale parameter of the forecast laplace distribution. The LS between obs and Laplace(location, scale). + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -541,15 +655,24 @@ def logs_laplace( >>> import scoringrules as sr >>> sr.logs_laplace(0.3, 0.1, 0.2) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale = B.asarray(scale) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the laplace must be positive." + ) return logarithmic.laplace(obs, location, scale, backend=backend) def logs_loglaplace( obs: "ArrayLike", - locationlog: "ArrayLike", - scalelog: "ArrayLike", + /, + locationlog: "ArrayLike" = 0.0, + scalelog: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the log-Laplace distribution. @@ -563,6 +686,11 @@ def logs_loglaplace( Location parameter of the forecast log-laplace distribution. scalelog : array_like Scale parameter of the forecast log-laplace distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -574,16 +702,24 @@ def logs_loglaplace( >>> import scoringrules as sr >>> sr.logs_loglaplace(3.0, 0.1, 0.9) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scalelog = B.asarray(scalelog) + if B.any(scalelog <= 0) or B.any(scalelog >= 1): + raise ValueError( + "`scalelog` contains entries outside of the range (0, 1). The scale parameter of the log-laplace distribution must be between 0 and 1." + ) return logarithmic.loglaplace(obs, locationlog, scalelog, backend=backend) def logs_logistic( obs: "ArrayLike", - mu: "ArrayLike", - sigma: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the logistic distribution. @@ -597,6 +733,11 @@ def logs_logistic( Location parameter of the forecast logistic distribution. sigma : array_like Scale parameter of the forecast logistic distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -608,14 +749,24 @@ def logs_logistic( >>> import scoringrules as sr >>> sr.logs_logistic(0.0, 0.4, 0.1) """ - return logarithmic.logistic(obs, mu, sigma, backend=backend) + if check_pars: + B = backends.active if backend is None else backends[backend] + scale = B.asarray(scale) + if B.any(scale <= 0): + raise ValueError( + "`sigma` contains non-positive entries. The scale parameter of the logistic distribution must be positive." + ) + return logarithmic.logistic(obs, location, scale, backend=backend) def logs_loglogistic( obs: "ArrayLike", - mulog: "ArrayLike", - sigmalog: "ArrayLike", + /, + mulog: "ArrayLike" = 0.0, + sigmalog: "ArrayLike" = 1.0, + *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the log-logistic distribution. @@ -630,7 +781,10 @@ def logs_loglogistic( sigmalog : array_like Scale parameter of the log-logistic distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -642,14 +796,24 @@ def logs_loglogistic( >>> import scoringrules as sr >>> sr.logs_loglogistic(3.0, 0.1, 0.9) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigmalog = B.asarray(sigmalog) + if B.any(sigmalog <= 0) or B.any(sigmalog >= 1): + raise ValueError( + "`sigmalog` contains entries outside of the range (0, 1). The scale parameter of the log-logistic distribution must be between 0 and 1." + ) return logarithmic.loglogistic(obs, mulog, sigmalog, backend=backend) def logs_lognormal( obs: "ArrayLike", - mulog: "ArrayLike", - sigmalog: "ArrayLike", + /, + mulog: "ArrayLike" = 0.0, + sigmalog: "ArrayLike" = 1.0, + *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the log-normal distribution. @@ -663,6 +827,11 @@ def logs_lognormal( Mean of the normal underlying distribution. sigmalog : array_like Standard deviation of the underlying normal distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -674,18 +843,26 @@ def logs_lognormal( >>> import scoringrules as sr >>> sr.logs_lognormal(0.0, 0.4, 0.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigmalog = B.asarray(sigmalog) + if B.any(sigmalog <= 0): + raise ValueError( + "`sigmalog` contains non-positive entries. The scale parameter of the log-normal distribution must be positive." + ) return logarithmic.lognormal(obs, mulog, sigmalog, backend=backend) def logs_mixnorm( obs: "ArrayLike", - m: "ArrayLike", - s: "ArrayLike", /, + m: "ArrayLike" = 0.0, + s: "ArrayLike" = 1.0, w: "ArrayLike" = None, mc_axis: "ArrayLike" = -1, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score for a mixture of normal distributions. @@ -704,7 +881,10 @@ def logs_mixnorm( mc_axis : int The axis corresponding to the mixture components. Default is the last axis. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -724,6 +904,15 @@ def logs_mixnorm( w = B.zeros(m.shape) + 1 / M else: w = B.asarray(w) + w = w / B.sum(w, axis=mc_axis, keepdims=True) + + if check_pars: + if B.any(s <= 0): + raise ValueError( + "`s` contains non-positive entries. The scale parameters of the normal distributions should be positive." + ) + if B.any(w < 0): + raise ValueError("`w` contains negative entries") if mc_axis != -1: m = B.moveaxis(m, mc_axis, -1) @@ -735,12 +924,13 @@ def logs_mixnorm( def logs_negbinom( obs: "ArrayLike", - n: "ArrayLike", /, + n: "ArrayLike", prob: "ArrayLike | None" = None, *, mu: "ArrayLike | None" = None, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the negative binomial distribution. @@ -756,6 +946,11 @@ def logs_negbinom( Probability parameter of the forecast negative binomial distribution. mu : array_like Mean of the forecast negative binomial distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -780,17 +975,24 @@ def logs_negbinom( if prob is None: prob = n / (n + mu) + if check_pars: + B = backends.active if backend is None else backends[backend] + prob = B.asarray(prob) + if B.any(prob < 0) or B.any(prob > 1): + raise ValueError("`prob` contains values outside the range [0, 1].") + return logarithmic.negbinom(obs, n, prob, backend=backend) def logs_normal( obs: "ArrayLike", - mu: "ArrayLike", - sigma: "ArrayLike", /, + mu: "ArrayLike" = 0.0, + sigma: "ArrayLike" = 1.0, *, negative: bool = True, backend: "Backend" = None, + check_pars: bool = False, ) -> "Array": r"""Compute the logarithmic score (LS) for the normal distribution. @@ -804,8 +1006,11 @@ def logs_normal( Mean of the forecast normal distribution. sigma : array_like Standard deviation of the forecast normal distribution. - backend : str, optional - The backend used for computations. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -817,17 +1022,25 @@ def logs_normal( >>> import scoringrules as sr >>> sr.logs_normal(0.0, 0.4, 0.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + sigma = B.asarray(sigma) + if B.any(sigma <= 0): + raise ValueError( + "`sigma` contains non-positive entries. The standard deviation of the normal distribution must be positive." + ) return logarithmic.normal(obs, mu, sigma, negative=negative, backend=backend) def logs_2pnormal( obs: "ArrayLike", - scale1: "ArrayLike", - scale2: "ArrayLike", - location: "ArrayLike", /, + scale1: "ArrayLike" = 1.0, + scale2: "ArrayLike" = 1.0, + location: "ArrayLike" = 0.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the two-piece normal distribution. @@ -844,7 +1057,10 @@ def logs_2pnormal( location : array_like Location parameter of the forecast two-piece normal distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -855,15 +1071,27 @@ def logs_2pnormal( >>> import scoringrules as sr >>> sr.logs_2pnormal(0.0, 0.4, 2.0, 0.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale1, scale2 = map(B.asarray, (scale1, scale2)) + if B.any(scale1 <= 0): + raise ValueError( + "`scale1` contains non-positive entries. The scale parameters of the two-piece normal distribution must be positive." + ) + if B.any(scale2 <= 0): + raise ValueError( + "`scale2` contains non-positive entries. The scale parameters of the two-piece normal distribution must be positive." + ) return logarithmic.twopnormal(obs, scale1, scale2, location, backend=backend) def logs_poisson( obs: "ArrayLike", - mean: "ArrayLike", /, + mean: "ArrayLike", *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the Poisson distribution. @@ -876,7 +1104,10 @@ def logs_poisson( mean : array_like Mean parameter of the forecast poisson distribution. backend : str - The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -888,17 +1119,25 @@ def logs_poisson( >>> import scoringrules as sr >>> sr.logs_poisson(1, 2) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + mean = B.asarray(mean) + if B.any(mean <= 0): + raise ValueError( + "`mean` contains non-positive entries. The mean parameter of the Poisson distribution must be positive." + ) return logarithmic.poisson(obs, mean, backend=backend) def logs_t( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the Student's t distribution. @@ -912,8 +1151,13 @@ def logs_t( Degrees of freedom parameter of the forecast t distribution. location : array_like Location parameter of the forecast t distribution. - sigma : array_like + scale : array_like Scale parameter of the forecast t distribution. + backend : str, optional + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -925,18 +1169,30 @@ def logs_t( >>> import scoringrules as sr >>> sr.logs_t(0.0, 0.1, 0.4, 0.1) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + df, scale = map(B.asarray, (df, scale)) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the t distribution must be positive." + ) return logarithmic.t(obs, df, location, scale, backend=backend) def logs_tlogistic( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the truncated logistic distribution. @@ -954,6 +1210,11 @@ def logs_tlogistic( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -965,18 +1226,28 @@ def logs_tlogistic( >>> import scoringrules as sr >>> sr.logs_tlogistic(0.0, 0.1, 0.4, -1.0, 1.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated logistic distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return logarithmic.tlogistic(obs, location, scale, lower, upper, backend=backend) def logs_tnormal( obs: "ArrayLike", - location: "ArrayLike", - scale: "ArrayLike", /, + location: "ArrayLike" = 0.0, + scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the truncated normal distribution. @@ -994,6 +1265,11 @@ def logs_tnormal( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1005,19 +1281,29 @@ def logs_tnormal( >>> import scoringrules as sr >>> sr.logs_tnormal(0.0, 0.1, 0.4, -1.0, 1.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated normal distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return logarithmic.tnormal(obs, location, scale, lower, upper, backend=backend) def logs_tt( obs: "ArrayLike", - df: "ArrayLike", /, + df: "ArrayLike", location: "ArrayLike" = 0.0, scale: "ArrayLike" = 1.0, lower: "ArrayLike" = float("-inf"), upper: "ArrayLike" = float("inf"), *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the truncated Student's t distribution. @@ -1037,6 +1323,11 @@ def logs_tt( Lower boundary of the truncated forecast distribution. upper : array_like Upper boundary of the truncated forecast distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- @@ -1048,16 +1339,30 @@ def logs_tt( >>> import scoringrules as sr >>> sr.logs_tt(0.0, 2.0, 0.1, 0.4, -1.0, 1.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + scale, lower, upper = map(B.asarray, (scale, lower, upper)) + if B.any(df <= 0): + raise ValueError( + "`df` contains non-positive entries. The degrees of freedom parameter of the truncated t distribution must be positive." + ) + if B.any(scale <= 0): + raise ValueError( + "`scale` contains non-positive entries. The scale parameter of the truncated t distribution must be positive." + ) + if B.any(lower >= upper): + raise ValueError("`lower` is not always smaller than `upper`.") return logarithmic.tt(obs, df, location, scale, lower, upper, backend=backend) def logs_uniform( obs: "ArrayLike", - min: "ArrayLike", - max: "ArrayLike", /, + min: "ArrayLike" = 0.0, + max: "ArrayLike" = 1.0, *, backend: "Backend" = None, + check_pars: bool = False, ) -> "ArrayLike": r"""Compute the logarithmic score (LS) for the uniform distribution. @@ -1071,17 +1376,27 @@ def logs_uniform( Lower bound of the forecast uniform distribution. max : array_like Upper bound of the forecast uniform distribution. + backend : str + The name of the backend used for computations. Defaults to ``numba`` if available, else ``numpy``. + check_pars: bool + Boolean indicating whether distribution parameter checks should be carried out prior to implementation. + Default is False. Returns ------- score : array_like - The LS between U(min, max, lmass, umass) and obs. + The LS between U(min, max) and obs. Examples -------- >>> import scoringrules as sr >>> sr.logs_uniform(0.4, 0.0, 1.0) """ + if check_pars: + B = backends.active if backend is None else backends[backend] + min, max = map(B.asarray, (min, max)) + if B.any(min >= max): + raise ValueError("`min` is not always smaller than `max`.") return logarithmic.uniform(obs, min, max, backend=backend) diff --git a/scoringrules/_variogram.py b/scoringrules/_variogram.py index b60db8e..0d88a9d 100644 --- a/scoringrules/_variogram.py +++ b/scoringrules/_variogram.py @@ -12,9 +12,11 @@ def vs_ensemble( obs: "Array", fct: "Array", /, + w: "Array" = None, m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, p: float = 1.0, backend: "Backend" = None, ) -> "Array": @@ -23,7 +25,7 @@ def vs_ensemble( For a :math:`D`-variate ensemble the Variogram Score [1]_ is defined as: .. math:: - \text{VS}_{p}(F_{ens}, \mathbf{y})= \sum_{i=1}^{d} \sum_{j=1}^{d} + \text{VS}_{p}(F_{ens}, \mathbf{y})= \sum_{i=1}^{d} \sum_{j=1}^{d} w_{i,j} \left( \frac{1}{M} \sum_{m=1}^{M} | x_{m,i} - x_{m,j} |^{p} - | y_{i} - y_{j} |^{p} \right)^{2}, where :math:`\mathbf{X}` and :math:`\mathbf{X'}` are independently sampled ensembles from from :math:`F`. @@ -36,12 +38,18 @@ def vs_ensemble( fct : array_like The predicted forecast ensemble, where the ensemble dimension is by default represented by the second last axis and the variables dimension by the last axis. - p : float - The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. + w : array_like + The weights assigned to pairs of dimensions. Must be of shape (..., D, D), where + D is the dimension, so that the weights are in the last two axes. m_axis : int The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. + p : float + The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. backend: str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -66,12 +74,29 @@ def vs_ensemble( >>> sr.vs_ensemble(obs, fct) array([ 8.65630139, 6.84693866, 19.52993307]) """ + B = backends.active if backend is None else backends[backend] obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) - if backend == "numba": - return variogram._variogram_score_gufunc(obs, fct, p) + if w is None: + D = fct.shape[-1] + w = B.zeros(obs.shape + (D,)) + 1.0 + + if ens_w is None: + if backend == "numba": + return variogram._variogram_score_gufunc(obs, fct, w, p) + + return variogram.vs(obs, fct, w, p, backend=backend) + + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + + if backend == "numba": + return variogram._variogram_score_gufunc_w(obs, fct, w, ens_w, p) - return variogram.vs(obs, fct, p, backend=backend) + return variogram.vs_w(obs, fct, w, ens_w, p, backend=backend) def twvs_ensemble( @@ -79,9 +104,11 @@ def twvs_ensemble( fct: "Array", v_func: tp.Callable, /, + w: "Array" = None, m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, p: float = 1.0, backend: "Backend" = None, ) -> "Array": @@ -103,14 +130,20 @@ def twvs_ensemble( fct : array_like The predicted forecast ensemble, where the ensemble dimension is by default represented by the second last axis and the variables dimension by the last axis. - p : float - The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. + w : array_like + The weights assigned to pairs of dimensions. Must be of shape (..., D, D), where + D is the dimension, so that the weights are in the last two axes. v_func : callable, array_like -> array_like Chaining function used to emphasise particular outcomes. m_axis : int The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. + p : float + The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -137,7 +170,7 @@ def twvs_ensemble( array([5.94996894, 4.72029765, 6.08947229]) """ obs, fct = map(v_func, (obs, fct)) - return vs_ensemble(obs, fct, m_axis, v_axis, p=p, backend=backend) + return vs_ensemble(obs, fct, w, m_axis, v_axis, ens_w=ens_w, p=p, backend=backend) def owvs_ensemble( @@ -145,9 +178,11 @@ def owvs_ensemble( fct: "Array", w_func: tp.Callable, /, + w: "Array" = None, m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, p: float = 1.0, backend: "Backend" = None, ) -> "Array": @@ -174,14 +209,20 @@ def owvs_ensemble( fct : array_like The predicted forecast ensemble, where the ensemble dimension is by default represented by the second last axis and the variables dimension by the last axis. - p : float - The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. + w : array_like + The weights assigned to pairs of dimensions. Must be of shape (..., D, D), where + D is the dimension, so that the weights are in the last two axes. m_axis : int The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. + p : float + The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -201,18 +242,39 @@ def owvs_ensemble( array([ 9.86816636, 6.75532522, 19.59353723]) """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) obs_weights = B.apply_along_axis(w_func, obs, -1) fct_weights = B.apply_along_axis(w_func, fct, -1) - if backend == "numba": - return variogram._owvariogram_score_gufunc( - obs, fct, p, obs_weights, fct_weights + if w is None: + D = fct.shape[-1] + w = B.zeros(obs.shape + (D,)) + 1.0 + + if ens_w is None: + if backend == "numba": + return variogram._owvariogram_score_gufunc( + obs, fct, w, obs_weights, fct_weights, p + ) + + return variogram.owvs( + obs, fct, w, obs_weights, fct_weights, p=p, backend=backend ) - return variogram.owvs(obs, fct, obs_weights, fct_weights, p=p, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + + if backend == "numba": + return variogram._owvariogram_score_gufunc_w( + obs, fct, w, obs_weights, fct_weights, ens_w, p + ) + + return variogram.owvs_w( + obs, fct, w, obs_weights, fct_weights, ens_w=ens_w, p=p, backend=backend + ) def vrvs_ensemble( @@ -220,9 +282,11 @@ def vrvs_ensemble( fct: "Array", w_func: tp.Callable, /, + w: "Array" = None, m_axis: int = -2, v_axis: int = -1, *, + ens_w: "Array" = None, p: float = 1.0, backend: "Backend" = None, ) -> "Array": @@ -251,14 +315,20 @@ def vrvs_ensemble( fct : array_like The predicted forecast ensemble, where the ensemble dimension is by default represented by the second last axis and the variables dimension by the last axis. - p : float - The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. w_func : callable, array_like -> array_like Weight function used to emphasise particular outcomes. + w : array_like + The weights assigned to pairs of dimensions. Must be of shape (..., D, D), where + D is the dimension, so that the weights are in the last two axes. m_axis : int The axis corresponding to the ensemble dimension. Defaults to -2. v_axis : int The axis corresponding to the variables dimension. Defaults to -1. + ens_w : array_like + Weights assigned to the ensemble members. Array with one less dimension than fct (without the v_axis dimension). + Default is equal weighting. + p : float + The order of the Variogram Score. Typical values are 0.5, 1.0 or 2.0. Defaults to 1.0. backend : str The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'. @@ -278,15 +348,36 @@ def vrvs_ensemble( array([46.48256493, 57.90759816, 92.37153472]) """ B = backends.active if backend is None else backends[backend] - obs, fct = multivariate_array_check(obs, fct, m_axis, v_axis, backend=backend) obs_weights = B.apply_along_axis(w_func, obs, -1) fct_weights = B.apply_along_axis(w_func, fct, -1) - if backend == "numba": - return variogram._vrvariogram_score_gufunc( - obs, fct, p, obs_weights, fct_weights + if w is None: + D = fct.shape[-1] + w = B.zeros(obs.shape + (D,)) + 1.0 + + if ens_w is None: + if backend == "numba": + return variogram._vrvariogram_score_gufunc( + obs, fct, w, obs_weights, fct_weights, p + ) + + return variogram.vrvs( + obs, fct, w, obs_weights, fct_weights, p=p, backend=backend ) - return variogram.vrvs(obs, fct, obs_weights, fct_weights, p=p, backend=backend) + else: + ens_w = B.asarray(ens_w) + if B.any(ens_w < 0): + raise ValueError("`ens_w` contains negative entries") + ens_w = ens_w / B.sum(ens_w, axis=-1, keepdims=True) + + if backend == "numba": + return variogram._vrvariogram_score_gufunc_w( + obs, fct, w, obs_weights, fct_weights, ens_w, p + ) + + return variogram.vrvs_w( + obs, fct, w, obs_weights, fct_weights, ens_w=ens_w, p=p, backend=backend + ) diff --git a/scoringrules/backend/base.py b/scoringrules/backend/base.py index 10d6200..64f0828 100644 --- a/scoringrules/backend/base.py +++ b/scoringrules/backend/base.py @@ -201,6 +201,17 @@ def any( ) -> "Array": """Test whether any input array element evaluates to ``True`` along a specified axis.""" + @abc.abstractmethod + def argsort( + self, + x: "Array", + /, + *, + axis: int = -1, + descending: bool = False, + ) -> "Array": + """Return the indices of a sorted copy of an input array ``x``.""" + @abc.abstractmethod def sort( self, @@ -297,3 +308,7 @@ def size(self, x: "Array") -> int: @abc.abstractmethod def indices(self, x: "Array") -> int: """Return an array representing the indices of a grid.""" + + @abc.abstractmethod + def gather(self, x: "Array", ind: "Array", axis: int) -> "Array": + """Reorder an array ``x`` depending on a template ``ind`` across an axis ``axis``.""" diff --git a/scoringrules/backend/jax.py b/scoringrules/backend/jax.py index 80332cb..2e5f872 100644 --- a/scoringrules/backend/jax.py +++ b/scoringrules/backend/jax.py @@ -194,6 +194,16 @@ def sort( out = jnp.sort(x, axis=axis) # TODO: this is slow! why? return -out if descending else out + def argsort( + self, + x: "Array", + /, + *, + axis: int = -1, + descending: bool = False, + ) -> "Array": + return jnp.argsort(x, axis=axis, descending=descending) + def norm(self, x: "Array", axis: int | tuple[int, ...] | None = None) -> "Array": return jnp.linalg.norm(x, axis=axis) @@ -269,6 +279,9 @@ def size(self, x: "Array") -> int: def indices(self, dimensions: tuple) -> "Array": return jnp.indices(dimensions) + def gather(self, x: "Array", ind: "Array", axis: int) -> "Array": + return jnp.take_along_axis(x, ind, axis=axis) + if __name__ == "__main__": B = JaxBackend() diff --git a/scoringrules/backend/numpy.py b/scoringrules/backend/numpy.py index 1fce50f..c1f9b08 100644 --- a/scoringrules/backend/numpy.py +++ b/scoringrules/backend/numpy.py @@ -198,6 +198,17 @@ def sort( out = np.sort(x, axis=axis) return -out if descending else out + def argsort( + self, + x: "NDArray", + /, + *, + axis: int = -1, + descending: bool = False, + ) -> "NDArray": + x = -x if descending else x + return np.argsort(x, axis=axis) + def norm( self, x: "NDArray", axis: int | tuple[int, ...] | None = None ) -> "NDArray": @@ -265,6 +276,9 @@ def size(self, x: "NDArray") -> int: def indices(self, dimensions: tuple) -> "NDArray": return np.indices(dimensions) + def gather(self, x: "NDArray", ind: "NDArray", axis: int) -> "NDArray": + return np.take_along_axis(x, ind, axis=axis) + class NumbaBackend(NumpyBackend): """Numba backend.""" diff --git a/scoringrules/backend/tensorflow.py b/scoringrules/backend/tensorflow.py index 9143607..453bfd6 100644 --- a/scoringrules/backend/tensorflow.py +++ b/scoringrules/backend/tensorflow.py @@ -225,6 +225,17 @@ def sort( direction = "DESCENDING" if descending else "ASCENDING" return tf.sort(x, axis=axis, direction=direction) + def argsort( + self, + x: "Tensor", + /, + *, + axis: int = -1, + descending: bool = False, + ) -> "Tensor": + direction = "DESCENDING" if descending else "ASCENDING" + return tf.argsort(x, axis=axis, direction=direction) + def norm(self, x: "Tensor", axis: int | tuple[int, ...] | None = None) -> "Tensor": return tf.norm(x, axis=axis) @@ -304,6 +315,10 @@ def indices(self, dimensions: tuple) -> "Tensor": indices = tf.stack(index_grids) return indices + def gather(self, x: "Tensor", ind: "Tensor", axis: int) -> "Tensor": + d = len(x.shape) + return tf.gather(x, ind, axis=axis, batch_dims=d) + if __name__ == "__main__": B = TensorflowBackend() diff --git a/scoringrules/backend/torch.py b/scoringrules/backend/torch.py index 3c7528b..1ff0e59 100644 --- a/scoringrules/backend/torch.py +++ b/scoringrules/backend/torch.py @@ -204,6 +204,17 @@ def sort( ) -> "Tensor": return torch.sort(x, stable=stable, dim=axis, descending=descending)[0] + def argsort( + self, + x: "Tensor", + /, + *, + axis: int = -1, + descending: bool = False, + stable: bool = True, + ) -> "Tensor": + return torch.argsort(x, stable=stable, dim=axis, descending=descending) + def norm(self, x: "Tensor", axis: int | tuple[int, ...] | None = None) -> "Tensor": return torch.norm(x, dim=axis) @@ -286,3 +297,6 @@ def indices(self, dimensions: tuple) -> "Tensor": index_grids = torch.meshgrid(*ranges, indexing="ij") indices = torch.stack(index_grids) return indices + + def gather(self, x: "Tensor", ind: "Tensor", axis: int) -> "Tensor": + return torch.gather(x, index=ind, dim=axis) diff --git a/scoringrules/core/brier.py b/scoringrules/core/brier.py index e66fbae..a721519 100644 --- a/scoringrules/core/brier.py +++ b/scoringrules/core/brier.py @@ -21,12 +21,13 @@ def brier_score( if not set(v.item() for v in B.unique_values(obs)) <= {0, 1}: raise ValueError("Observations must be 0, 1, or NaN.") - return B.asarray((fct - obs) ** 2) + return (fct - obs) ** 2 def rps_score( obs: "ArrayLike", fct: "ArrayLike", + onehot: bool = False, backend: "Backend" = None, ) -> "Array": """Compute the Ranked Probability Score for ordinal categorical forecasts.""" @@ -36,12 +37,11 @@ def rps_score( if B.any(fct < 0.0) or B.any(fct > 1.0 + EPSILON): raise ValueError("Forecast probabilities must be between 0 and 1.") - categories = B.arange(1, fct.shape[-1] + 1) - obs_one_hot = B.where(B.expand_dims(obs, -1) == categories, 1, 0) + if not onehot: + categories = B.arange(1, fct.shape[-1] + 1) + obs = B.where(B.expand_dims(obs, -1) == categories, 1, 0) - return B.sum( - (B.cumsum(fct, axis=-1) - B.cumsum(obs_one_hot, axis=-1)) ** 2, axis=-1 - ) + return B.sum((B.cumsum(fct, axis=-1) - B.cumsum(obs, axis=-1)) ** 2, axis=-1) def log_score(obs: "ArrayLike", fct: "ArrayLike", backend: "Backend" = None) -> "Array": diff --git a/scoringrules/core/crps/__init__.py b/scoringrules/core/crps/__init__.py index 74f563a..d13f562 100644 --- a/scoringrules/core/crps/__init__.py +++ b/scoringrules/core/crps/__init__.py @@ -1,4 +1,5 @@ -from ._approx import ensemble, ow_ensemble, quantile_pinball, vr_ensemble +from ._approx import ensemble, ow_ensemble, vr_ensemble, quantile_pinball +from ._approx_w import ensemble_w, ow_ensemble_w, vr_ensemble_w from ._closed import ( beta, binomial, @@ -27,14 +28,19 @@ try: from ._gufuncs import estimator_gufuncs, quantile_pinball_gufunc + from ._gufuncs_w import estimator_gufuncs_w except ImportError: estimator_gufuncs = None + estimator_gufuncs_w = None quantile_pinball_gufunc = None __all__ = [ "ensemble", + "ensemble_w", "ow_ensemble", + "ow_ensemble_w", "vr_ensemble", + "vr_ensemble_w", "beta", "binomial", "exponential", @@ -59,6 +65,7 @@ "t", "uniform", "estimator_gufuncs", + "estimator_gufuncs_w", "quantile_pinball", "quantile_pinball_gufunc", ] diff --git a/scoringrules/core/crps/_approx.py b/scoringrules/core/crps/_approx.py index 8a40abd..18564dc 100644 --- a/scoringrules/core/crps/_approx.py +++ b/scoringrules/core/crps/_approx.py @@ -19,6 +19,12 @@ def ensemble( out = _crps_ensemble_pwm(obs, fct, backend=backend) elif estimator == "fair": out = _crps_ensemble_fair(obs, fct, backend=backend) + elif estimator == "qd": + out = _crps_ensemble_qd(obs, fct, backend=backend) + elif estimator == "akr": + out = _crps_ensemble_akr(obs, fct, backend=backend) + elif estimator == "akr_circperm": + out = _crps_ensemble_akr_circperm(obs, fct, backend=backend) else: raise ValueError( f"{estimator} can only be used with `numpy` " @@ -33,7 +39,7 @@ def _crps_ensemble_fair( ) -> "Array": """Fair version of the CRPS estimator based on the energy form.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-1] + M = fct.shape[-1] e_1 = B.sum(B.abs(obs[..., None] - fct), axis=-1) / M e_2 = B.sum( B.abs(fct[..., None] - fct[..., None, :]), @@ -65,6 +71,42 @@ def _crps_ensemble_pwm( return expected_diff + β_0 - 2.0 * β_1 +def _crps_ensemble_qd(obs: "Array", fct: "Array", backend: "Backend" = None) -> "Array": + """CRPS estimator based on the quantile decomposition form.""" + B = backends.active if backend is None else backends[backend] + M: int = fct.shape[-1] + alpha = B.arange(1, M + 1) - 0.5 + below = (fct <= obs[..., None]) * alpha * (obs[..., None] - fct) + above = (fct > obs[..., None]) * (M - alpha) * (fct - obs[..., None]) + out = B.sum(below + above, axis=-1) / (M**2) + return 2 * out + + +def _crps_ensemble_akr( + obs: "Array", fct: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the approximate kernel representation.""" + B = backends.active if backend is None else backends[backend] + M = fct.shape[-1] + e_1 = B.mean(B.abs(obs[..., None] - fct), axis=-1) + ind = [(i + 1) % M for i in range(M)] + e_2 = B.mean(B.abs(fct[..., ind] - fct)[..., ind], axis=-1) + return e_1 - 0.5 * e_2 + + +def _crps_ensemble_akr_circperm( + obs: "Array", fct: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the AKR with cyclic permutation.""" + B = backends.active if backend is None else backends[backend] + M = fct.shape[-1] + e_1 = B.mean(B.abs(obs[..., None] - fct), axis=-1) + shift = int((M - 1) / 2) + ind = [(i + shift) % M for i in range(M)] + e_2 = B.mean(B.abs(fct[..., ind] - fct), axis=-1) + return e_1 - 0.5 * e_2 + + def quantile_pinball( obs: "Array", fct: "Array", alpha: "Array", backend: "Backend" = None ) -> "Array": @@ -82,16 +124,15 @@ def ow_ensemble( fw: "Array", backend: "Backend" = None, ) -> "Array": - """Outcome-Weighted CRPS estimator based on the energy form.""" + """Outcome-Weighted CRPS for an ensemble forecast.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-1] wbar = B.mean(fw, axis=-1) - e_1 = B.sum(B.abs(obs[..., None] - fct) * fw, axis=-1) * ow / (M * wbar) - e_2 = B.sum( + e_1 = B.mean(B.abs(obs[..., None] - fct) * fw, axis=-1) * ow / wbar + e_2 = B.mean( B.abs(fct[..., None] - fct[..., None, :]) * fw[..., None] * fw[..., None, :], axis=(-1, -2), ) - e_2 *= ow / (M**2 * wbar**2) + e_2 *= ow / (wbar**2) return e_1 - 0.5 * e_2 @@ -102,15 +143,13 @@ def vr_ensemble( fw: "Array", backend: "Backend" = None, ) -> "Array": - """Vertically Re-scaled CRPS estimator based on the energy form.""" + """Vertically Re-scaled CRPS for an ensemble forecast.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-1] - e_1 = B.sum(B.abs(obs[..., None] - fct) * fw, axis=-1) * ow / M - e_2 = B.sum( - B.abs(B.expand_dims(fct, axis=-1) - B.expand_dims(fct, axis=-2)) - * (B.expand_dims(fw, axis=-1) * B.expand_dims(fw, axis=-2)), + e_1 = B.mean(B.abs(obs[..., None] - fct) * fw, axis=-1) * ow + e_2 = B.mean( + B.abs(fct[..., None] - fct[..., None, :]) * fw[..., None] * fw[..., None, :], axis=(-1, -2), - ) / (M**2) + ) e_3 = B.mean(B.abs(fct) * fw, axis=-1) - B.abs(obs) * ow e_3 *= B.mean(fw, axis=1) - ow return e_1 - 0.5 * e_2 + e_3 diff --git a/scoringrules/core/crps/_approx_w.py b/scoringrules/core/crps/_approx_w.py new file mode 100644 index 0000000..d342da2 --- /dev/null +++ b/scoringrules/core/crps/_approx_w.py @@ -0,0 +1,159 @@ +import typing as tp + +from scoringrules.backend import backends + +if tp.TYPE_CHECKING: + from scoringrules.core.typing import Array, ArrayLike, Backend + + +def ensemble_w( + obs: "ArrayLike", + fct: "Array", + ens_w: "Array" = None, + estimator: str = "pwm", + backend: "Backend" = None, +) -> "Array": + """Compute the CRPS for a finite weighted ensemble.""" + if estimator == "nrg": + out = _crps_ensemble_nrg_w(obs, fct, ens_w, backend=backend) + elif estimator == "pwm": + out = _crps_ensemble_pwm_w(obs, fct, ens_w, backend=backend) + elif estimator == "fair": + out = _crps_ensemble_fair_w(obs, fct, ens_w, backend=backend) + elif estimator == "qd": + out = _crps_ensemble_qd_w(obs, fct, ens_w, backend=backend) + elif estimator == "akr": + out = _crps_ensemble_akr_w(obs, fct, ens_w, backend=backend) + elif estimator == "akr_circperm": + out = _crps_ensemble_akr_circperm_w(obs, fct, ens_w, backend=backend) + else: + raise ValueError( + f"{estimator} can only be used with `numpy` " + "backend and needs `numba` to be installed" + ) + + return out + + +def _crps_ensemble_fair_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """Fair version of the CRPS estimator based on the energy form.""" + B = backends.active if backend is None else backends[backend] + e_1 = B.sum(B.abs(obs[..., None] - fct) * w, axis=-1) + e_2 = B.sum( + B.abs(fct[..., None] - fct[..., None, :]) * w[..., None] * w[..., None, :], + axis=(-1, -2), + ) / (1 - B.sum(w * w, axis=-1)) + return e_1 - 0.5 * e_2 + + +def _crps_ensemble_nrg_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the energy form.""" + B = backends.active if backend is None else backends[backend] + e_1 = B.sum(B.abs(obs[..., None] - fct) * w, axis=-1) + e_2 = B.sum( + B.abs(fct[..., None] - fct[..., None, :]) * w[..., None] * w[..., None, :], + (-1, -2), + ) + return e_1 - 0.5 * e_2 + + +def _crps_ensemble_pwm_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the probability weighted moment (PWM) form.""" + B = backends.active if backend is None else backends[backend] + w_sum = B.cumsum(w, axis=-1) + expected_diff = B.sum(B.abs(obs[..., None] - fct) * w, axis=-1) + β_0 = B.sum(fct * w, axis=-1) + β_1 = B.sum(fct * w * (w_sum - w), axis=-1) / (1 - B.sum(w * w, axis=-1)) + return expected_diff + β_0 - 2.0 * β_1 + + +def _crps_ensemble_qd_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the quantile score decomposition.""" + B = backends.active if backend is None else backends[backend] + w_sum = B.cumsum(w, axis=-1) + a = w_sum - 0.5 * w + dif = fct - obs[..., None] + c = B.where(dif > 0, 1 - a, -a) + s = B.sum(w * c * dif, axis=-1) + return 2 * s + + +def _crps_ensemble_akr_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the approximate kernel representation.""" + B = backends.active if backend is None else backends[backend] + M = fct.shape[-1] + e_1 = B.sum(B.abs(obs[..., None] - fct) * w, axis=-1) + ind = [(i + 1) % M for i in range(M)] + e_2 = B.sum(B.abs(fct[..., ind] - fct) * w[..., ind], axis=-1) + return e_1 - 0.5 * e_2 + + +def _crps_ensemble_akr_circperm_w( + obs: "Array", fct: "Array", w: "Array", backend: "Backend" = None +) -> "Array": + """CRPS estimator based on the AKR with cyclic permutation.""" + B = backends.active if backend is None else backends[backend] + M = fct.shape[-1] + e_1 = B.sum(B.abs(obs[..., None] - fct) * w, axis=-1) + shift = int((M - 1) / 2) + ind = [(i + shift) % M for i in range(M)] + e_2 = B.sum(B.abs(fct[..., ind] - fct) * w[..., ind], axis=-1) + return e_1 - 0.5 * e_2 + + +def ow_ensemble_w( + obs: "Array", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Outcome-Weighted CRPS for an ensemble forecast.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(ens_w * fw, axis=-1) + e_1 = B.sum(ens_w * B.abs(obs[..., None] - fct) * fw, axis=-1) * ow / wbar + e_2 = B.sum( + ens_w[..., None] + * ens_w[..., None, :] + * B.abs(fct[..., None] - fct[..., None, :]) + * fw[..., None] + * fw[..., None, :], + axis=(-1, -2), + ) + e_2 *= ow / (wbar**2) + return e_1 - 0.5 * e_2 + + +def vr_ensemble_w( + obs: "Array", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Vertically Re-scaled CRPS for an ensemble forecast.""" + B = backends.active if backend is None else backends[backend] + e_1 = B.sum(ens_w * B.abs(obs[..., None] - fct) * fw, axis=-1) * ow + e_2 = B.sum( + ens_w[..., None] + * ens_w[..., None, :] + * B.abs(fct[..., None] - fct[..., None, :]) + * fw[..., None] + * fw[..., None, :], + axis=(-1, -2), + ) + e_3 = B.sum(ens_w * B.abs(fct) * fw, axis=-1) - B.abs(obs) * ow + e_3 *= B.sum(ens_w * fw, axis=1) - ow + return e_1 - 0.5 * e_2 + e_3 diff --git a/scoringrules/core/crps/_closed.py b/scoringrules/core/crps/_closed.py index 4fc7822..cd98084 100644 --- a/scoringrules/core/crps/_closed.py +++ b/scoringrules/core/crps/_closed.py @@ -35,28 +35,19 @@ def beta( """Compute the CRPS for the beta distribution.""" B = backends.active if backend is None else backends[backend] obs, a, b, lower, upper = map(B.asarray, (obs, a, b, lower, upper)) + a = B.where(a > 0.0, a, B.nan) + b = B.where(b > 0.0, b, B.nan) + lower = B.where(lower < upper, lower, B.nan) - if _is_scalar_value(lower, 0.0) and _is_scalar_value(upper, 1.0): - special_limits = False - else: - if B.any(lower >= upper): - raise ValueError("lower must be less than upper") - special_limits = True - - if special_limits: - obs = (obs - lower) / (upper - lower) + obs = (obs - lower) / (upper - lower) + obs_std = B.minimum(B.maximum(obs, 0.0), 1.0) - I_ab = B.betainc(a, b, obs) - I_a1b = B.betainc(a + 1, b, obs) - F_ab = B.minimum(B.maximum(I_ab, 0), 1) - F_a1b = B.minimum(B.maximum(I_a1b, 0), 1) + F_ab = B.betainc(a, b, obs_std) + F_a1b = B.betainc(a + 1, b, obs_std) bet_rat = 2 * B.beta(2 * a, 2 * b) / (a * B.beta(a, b) ** 2) s = obs * (2 * F_ab - 1) + (a / (a + b)) * (1 - 2 * F_a1b - bet_rat) - if special_limits: - s = s * (upper - lower) - - return s + return s * (upper - lower) def binomial( @@ -133,26 +124,18 @@ def exponential( def exponentialM( obs: "ArrayLike", - mass: "ArrayLike", location: "ArrayLike", scale: "ArrayLike", + mass: "ArrayLike", backend: "Backend" = None, ) -> "Array": """Compute the CRPS for the standard exponential distribution with a point mass at the boundary.""" B = backends.active if backend is None else backends[backend] obs, location, scale, mass = map(B.asarray, (obs, location, scale, mass)) - - if not _is_scalar_value(location, 0.0): - obs -= location - - a = 1.0 if _is_scalar_value(mass, 0.0) else 1 - mass + obs -= location + a = 1.0 - mass s = B.abs(obs) - - if _is_scalar_value(scale, 1.0): - s -= a * (2 * _exp_cdf(obs, 1.0, backend=backend) - 0.5 * a) - else: - s -= scale * a * (2 * _exp_cdf(obs, 1 / scale, backend=backend) - 0.5 * a) - + s -= scale * a * (2 * _exp_cdf(obs, 1 / scale, backend=backend) - 0.5 * a) return s @@ -266,6 +249,7 @@ def gpd( ) shape = B.where(shape < 1.0, shape, B.nan) mass = B.where((mass >= 0.0) & (mass <= 1.0), mass, B.nan) + scale = B.where(scale > 0.0, scale, B.nan) ω = (obs - location) / scale F_xi = _gpd_cdf(ω, shape, backend=backend) s = ( @@ -536,7 +520,7 @@ def logistic( sigma: "ArrayLike", backend: "Backend" = None, ) -> "Array": - """Compute the CRPS for the normal distribution.""" + """Compute the CRPS for the logistic distribution.""" B = backends.active if backend is None else backends[backend] mu, sigma, obs = map(B.asarray, (mu, sigma, obs)) ω = (obs - mu) / sigma @@ -665,7 +649,7 @@ def normal( sigma: "ArrayLike", backend: "Backend" = None, ) -> "Array": - """Compute the CRPS for the logistic distribution.""" + """Compute the CRPS for the normal distribution.""" B = backends.active if backend is None else backends[backend] mu, sigma, obs = map(B.asarray, (mu, sigma, obs)) ω = (obs - mu) / sigma diff --git a/scoringrules/core/crps/_gufuncs.py b/scoringrules/core/crps/_gufuncs.py index 1e2ed1d..2d54d80 100644 --- a/scoringrules/core/crps/_gufuncs.py +++ b/scoringrules/core/crps/_gufuncs.py @@ -1,7 +1,5 @@ -import math - import numpy as np -from numba import guvectorize, njit, vectorize +from numba import guvectorize from scoringrules.core.utils import lazy_gufunc_wrapper_uv @@ -99,11 +97,9 @@ def _crps_ensemble_nrg_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray) """CRPS estimator based on the energy form.""" obs = obs[0] M = fct.shape[-1] - if np.isnan(obs): out[0] = np.nan return - e_1 = 0 e_2 = 0 @@ -137,9 +133,9 @@ def _crps_ensemble_fair_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray for i in range(M): e_1 += abs(fct[i] - obs) for j in range(i + 1, M): - e_2 += 2 * abs(fct[j] - fct[i]) + e_2 += abs(fct[j] - fct[i]) - out[0] = e_1 / M - 0.5 * e_2 / (M * (M - 1)) + out[0] = e_1 / M - e_2 / (M * (M - 1)) @guvectorize("(),(n)->()") @@ -174,6 +170,7 @@ def _crps_ensemble_akr_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray) i = M - 1 e_1 += abs(forecast - obs) e_2 += abs(forecast - fct[i - 1]) + out[0] = e_1 / M - 0.5 * 1 / M * e_2 @@ -202,11 +199,9 @@ def _owcrps_ensemble_nrg_gufunc( ): """Outcome-weighted CRPS estimator based on the energy form.""" M = fct.shape[-1] - if np.isnan(obs): out[0] = np.nan return - e_1 = 0.0 e_2 = 0.0 @@ -230,11 +225,9 @@ def _vrcrps_ensemble_nrg_gufunc( ): """Vertically re-scaled CRPS estimator based on the energy form.""" M = fct.shape[-1] - if np.isnan(obs): out[0] = np.nan return - e_1 = 0.0 e_2 = 0.0 @@ -250,51 +243,6 @@ def _vrcrps_ensemble_nrg_gufunc( out[0] = e_1 / M - 0.5 * e_2 / (M**2) + (wabs_x - wabs_y) * (wbar - ow) -@njit(["float32(float32)", "float64(float64)"]) -def _norm_cdf(x: float) -> float: - """Cumulative distribution function for the standard normal distribution.""" - out: float = (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0 - return out - - -@njit(["float32(float32)", "float64(float64)"]) -def _norm_pdf(x: float) -> float: - """Probability density function for the standard normal distribution.""" - out: float = (1 / math.sqrt(2 * math.pi)) * math.exp(-(x**2) / 2) - return out - - -@njit(["float32(float32)", "float64(float64)"]) -def _logis_cdf(x: float) -> float: - """Cumulative distribution function for the standard logistic distribution.""" - out: float = 1.0 / (1.0 + math.exp(-x)) - return out - - -@vectorize(["float32(float32, float32, float32)", "float64(float64, float64, float64)"]) -def _crps_normal_ufunc(obs: float, mu: float, sigma: float) -> float: - ω = (obs - mu) / sigma - out: float = sigma * (ω * (2 * _norm_cdf(ω) - 1) + 2 * _norm_pdf(ω) - INV_SQRT_PI) - return out - - -@vectorize(["float32(float32, float32, float32)", "float64(float64, float64, float64)"]) -def _crps_lognormal_ufunc(obs: float, mulog: float, sigmalog: float) -> float: - ω = (np.log(obs) - mulog) / sigmalog - ex = 2 * np.exp(mulog + sigmalog**2 / 2) - out: float = obs * (2 * _norm_cdf(ω) - 1) - ex * ( - _norm_cdf(ω - sigmalog) + _norm_cdf(sigmalog / np.sqrt(2)) - 1 - ) - return out - - -@vectorize(["float32(float32, float32, float32)", "float64(float64, float64, float64)"]) -def _crps_logistic_ufunc(obs: float, mu: float, sigma: float) -> float: - ω = (obs - mu) / sigma - out: float = sigma * (ω - 2 * np.log(_logis_cdf(ω)) - 1) - return out - - estimator_gufuncs = { "akr_circperm": lazy_gufunc_wrapper_uv(_crps_ensemble_akr_circperm_gufunc), "akr": lazy_gufunc_wrapper_uv(_crps_ensemble_akr_gufunc), @@ -315,8 +263,5 @@ def _crps_logistic_ufunc(obs: float, mu: float, sigma: float) -> float: "_crps_ensemble_nrg_gufunc", "_crps_ensemble_pwm_gufunc", "_crps_ensemble_qd_gufunc", - "_crps_normal_ufunc", - "_crps_lognormal_ufunc", - "_crps_logistic_ufunc", "quantile_pinball_gufunc", ] diff --git a/scoringrules/core/crps/_gufuncs_w.py b/scoringrules/core/crps/_gufuncs_w.py new file mode 100644 index 0000000..f4ccaad --- /dev/null +++ b/scoringrules/core/crps/_gufuncs_w.py @@ -0,0 +1,253 @@ +import numpy as np +from numba import guvectorize + +from scoringrules.core.utils import lazy_gufunc_wrapper_uv + +INV_SQRT_PI = 1 / np.sqrt(np.pi) +EPSILON = 1e-6 + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_int_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimator based on the integral form.""" + + if np.isnan(obs): + out[0] = np.nan + return + + obs_cdf = 0 + forecast_cdf = 0.0 + prev_forecast = 0.0 + integral = 0.0 + + for n, forecast in enumerate(fct): + if np.isnan(forecast): + if n == 0: + integral = np.nan + forecast = prev_forecast # noqa: PLW2901 + break + + if obs_cdf == 0 and obs < forecast: + # this correctly handles the transition point of the obs CDF + integral += (obs - prev_forecast) * forecast_cdf**2 + integral += (forecast - obs) * (forecast_cdf - 1) ** 2 + obs_cdf = 1 + else: + integral += (forecast_cdf - obs_cdf) ** 2 * (forecast - prev_forecast) + + forecast_cdf += w[n] + prev_forecast = forecast + + if obs_cdf == 0: + integral += obs - forecast + + out[0] = integral + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_qd_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimator based on the quantile decomposition form.""" + + if np.isnan(obs): + out[0] = np.nan + return + + obs_cdf = 0.0 + integral = 0.0 + a = np.cumsum(w) - w / 2 + + for i, forecast in enumerate(fct): + if obs < forecast: + obs_cdf = 1.0 + + integral += w[i] * (obs_cdf - a[i]) * (forecast - obs) + + out[0] = 2 * integral + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_nrg_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimator based on the energy form.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0 + e_2 = 0 + + for i in range(M): + e_1 += abs(fct[i] - obs) * w[i] + for j in range(i + 1, M): + e_2 += abs(fct[j] - fct[i]) * w[j] * w[i] + + out[0] = e_1 - e_2 + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_fair_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """Fair version of the CRPS estimator based on the energy form.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0 + e_2 = 0 + + for i in range(M): + e_1 += abs(fct[i] - obs) * w[i] + for j in range(i + 1, M): + e_2 += abs(fct[j] - fct[i]) * w[j] * w[i] + + fair_c = 1 - np.sum(w**2) + + out[0] = e_1 - e_2 / fair_c + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_pwm_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimator based on the probability weighted moment (PWM) form.""" + + if np.isnan(obs): + out[0] = np.nan + return + + expected_diff = 0.0 + β_0 = 0.0 + β_1 = 0.0 + + w_sum = np.cumsum(w) + + for i, forecast in enumerate(fct): + expected_diff += np.abs(forecast - obs) * w[i] + β_0 += forecast * w[i] * (1.0 - w[i]) + β_1 += forecast * w[i] * (w_sum[i] - w[i]) + + out[0] = expected_diff + β_0 - 2 * β_1 + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_akr_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimaton based on the approximate kernel representation.""" + M = fct.shape[-1] + e_1 = 0 + e_2 = 0 + for i, forecast in enumerate(fct): + if i == 0: + i = M - 1 + e_1 += abs(forecast - obs) * w[i] + e_2 += abs(forecast - fct[i - 1]) * w[i] + out[0] = e_1 - 0.5 * e_2 + + +@guvectorize("(),(n),(n)->()") +def _crps_ensemble_akr_circperm_w_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """CRPS estimaton based on the AKR with cyclic permutation.""" + M = fct.shape[-1] + e_1 = 0.0 + e_2 = 0.0 + for i, forecast in enumerate(fct): + sigma_i = int((i + 1 + ((M - 1) / 2)) % M) + e_1 += abs(forecast - obs) * w[i] + e_2 += abs(forecast - fct[sigma_i]) * w[i] + out[0] = e_1 - 0.5 * e_2 + + +@guvectorize("(),(n),(),(n),(n)->()") +def _owcrps_ensemble_nrg_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + ens_w: np.ndarray, + out: np.ndarray, +): + """Outcome-weighted CRPS estimator based on the energy form.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0.0 + e_2 = 0.0 + + for i in range(M): + e_1 += ens_w[i] * abs(fct[i] - obs) * fw[i] * ow + for j in range(i + 1, M): + e_2 += 2 * ens_w[i] * ens_w[j] * abs(fct[i] - fct[j]) * fw[i] * fw[j] * ow + + wbar = np.sum(ens_w * fw) + + out[0] = e_1 / wbar - 0.5 * e_2 / (wbar**2) + + +@guvectorize("(),(n),(),(n),(n)->()") +def _vrcrps_ensemble_nrg_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + ens_w: np.ndarray, + out: np.ndarray, +): + """Vertically re-scaled CRPS estimator based on the energy form.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0.0 + e_2 = 0.0 + + for i in range(M): + e_1 += ens_w[i] * abs(fct[i] - obs) * fw[i] * ow + for j in range(i + 1, M): + e_2 += 2 * ens_w[i] * ens_w[j] * abs(fct[i] - fct[j]) * fw[i] * fw[j] + + wbar = np.sum(ens_w * fw) + wabs_x = np.sum(ens_w * np.abs(fct) * fw) + wabs_y = abs(obs) * ow + + out[0] = e_1 - 0.5 * e_2 + (wabs_x - wabs_y) * (wbar - ow) + + +estimator_gufuncs_w = { + "akr_circperm": lazy_gufunc_wrapper_uv(_crps_ensemble_akr_circperm_w_gufunc), + "akr": lazy_gufunc_wrapper_uv(_crps_ensemble_akr_w_gufunc), + "fair": lazy_gufunc_wrapper_uv(_crps_ensemble_fair_w_gufunc), + "int": lazy_gufunc_wrapper_uv(_crps_ensemble_int_w_gufunc), + "nrg": lazy_gufunc_wrapper_uv(_crps_ensemble_nrg_w_gufunc), + "pwm": lazy_gufunc_wrapper_uv(_crps_ensemble_pwm_w_gufunc), + "qd": lazy_gufunc_wrapper_uv(_crps_ensemble_qd_w_gufunc), + "ownrg": lazy_gufunc_wrapper_uv(_owcrps_ensemble_nrg_w_gufunc), + "vrnrg": lazy_gufunc_wrapper_uv(_vrcrps_ensemble_nrg_w_gufunc), +} + +__all__ = [ + "_crps_ensemble_akr_circperm_w_gufunc", + "_crps_ensemble_akr_w_gufunc", + "_crps_ensemble_fair_w_gufunc", + "_crps_ensemble_int_w_gufunc", + "_crps_ensemble_nrg_w_gufunc", + "_crps_ensemble_pwm_w_gufunc", + "_crps_ensemble_qd_w_gufunc", +] diff --git a/scoringrules/core/energy/__init__.py b/scoringrules/core/energy/__init__.py index b9bc987..134bb2b 100644 --- a/scoringrules/core/energy/__init__.py +++ b/scoringrules/core/energy/__init__.py @@ -9,15 +9,40 @@ _owenergy_score_gufunc = None _vrenergy_score_gufunc = None +try: + from ._gufuncs_w import ( + _energy_score_gufunc_w, + _owenergy_score_gufunc_w, + _vrenergy_score_gufunc_w, + ) + +except ImportError: + _energy_score_gufunc = None + _energy_score_gufunc_w = None + _owenergy_score_gufunc = None + _owenergy_score_gufunc_w = None + _vrenergy_score_gufunc = None + _vrenergy_score_gufunc_w = None + from ._score import es_ensemble as nrg from ._score import owes_ensemble as ownrg from ._score import vres_ensemble as vrnrg +from ._score_w import es_ensemble_w as nrg_w +from ._score_w import owes_ensemble_w as ownrg_w +from ._score_w import vres_ensemble_w as vrnrg_w + __all__ = [ "nrg", + "nrg_w", "ownrg", + "ownrg_w", "vrnrg", + "vrnrg_w", "_energy_score_gufunc", + "_energy_score_gufunc_w", "_owenergy_score_gufunc", + "_owenergy_score_gufunc_w", "_vrenergy_score_gufunc", + "_vrenergy_score_gufunc_w", ] diff --git a/scoringrules/core/energy/_gufuncs.py b/scoringrules/core/energy/_gufuncs.py index 13aae9c..630ff5b 100644 --- a/scoringrules/core/energy/_gufuncs.py +++ b/scoringrules/core/energy/_gufuncs.py @@ -19,7 +19,6 @@ def _energy_score_gufunc( ): """Compute the Energy Score for a finite ensemble.""" M = fct.shape[0] - e_1 = 0.0 e_2 = 0.0 for i in range(M): diff --git a/scoringrules/core/energy/_gufuncs_w.py b/scoringrules/core/energy/_gufuncs_w.py new file mode 100644 index 0000000..d856aba --- /dev/null +++ b/scoringrules/core/energy/_gufuncs_w.py @@ -0,0 +1,90 @@ +import numpy as np +from numba import guvectorize + +from scoringrules.core.utils import lazy_gufunc_wrapper_mv + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(m)->()") +def _energy_score_gufunc_w( + obs: np.ndarray, + fct: np.ndarray, + ens_w: np.ndarray, + out: np.ndarray, +): + """Compute the Energy Score for a finite ensemble.""" + M = fct.shape[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(np.linalg.norm(fct[i] - obs)) * ens_w[i] + for j in range(i + 1, M): + e_2 += 2 * float(np.linalg.norm(fct[i] - fct[j])) * ens_w[i] * ens_w[j] + + out[0] = e_1 - 0.5 * e_2 + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(),(m),(m)->()") +def _owenergy_score_gufunc_w( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + ens_w: np.ndarray, + out: np.ndarray, +): + """Compute the Outcome-Weighted Energy Score for a finite ensemble.""" + M = fct.shape[0] + ow = ow[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(np.linalg.norm(fct[i] - obs) * fw[i] * ow) * ens_w[i] + for j in range(i + 1, M): + e_2 += ( + 2 + * float(np.linalg.norm(fct[i] - fct[j]) * fw[i] * fw[j] * ow) + * ens_w[i] + * ens_w[j] + ) + + wbar = np.mean(fw) + + out[0] = e_1 / (wbar) - 0.5 * e_2 / (wbar**2) + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(),(m),(m)->()") +def _vrenergy_score_gufunc_w( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + ens_w: np.ndarray, + out: np.ndarray, +): + """Compute the Vertically Re-scaled Energy Score for a finite ensemble.""" + M = fct.shape[0] + ow = ow[0] + + e_1 = 0.0 + e_2 = 0.0 + wabs_x = 0.0 + for i in range(M): + e_1 += float(np.linalg.norm(fct[i] - obs) * fw[i] * ow) * ens_w[i] + wabs_x += np.linalg.norm(fct[i]) * fw[i] * ens_w[i] + for j in range(i + 1, M): + e_2 += ( + 2 + * float(np.linalg.norm(fct[i] - fct[j]) * fw[i] * fw[j]) + * ens_w[i] + * ens_w[j] + ) + + wbar = np.sum(fw * ens_w) + wabs_y = np.linalg.norm(obs) * ow + + out[0] = e_1 - 0.5 * e_2 + (wabs_x - wabs_y) * (wbar - ow) diff --git a/scoringrules/core/energy/_score.py b/scoringrules/core/energy/_score.py index e264c4b..8bf1059 100644 --- a/scoringrules/core/energy/_score.py +++ b/scoringrules/core/energy/_score.py @@ -13,13 +13,13 @@ def es_ensemble(obs: "Array", fct: "Array", backend=None) -> "Array": The ensemble and variables axes are on the second last and last dimensions respectively. """ B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-2] err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) - E_1 = B.sum(err_norm, -1) / M + E_1 = B.mean(err_norm, axis=-1) spread_norm = B.norm(B.expand_dims(fct, -3) - B.expand_dims(fct, -2), -1) - E_2 = B.sum(spread_norm, (-2, -1)) / (M**2) + E_2 = B.mean(spread_norm, axis=(-2, -1)) + return E_1 - 0.5 * E_2 @@ -32,18 +32,17 @@ def owes_ensemble( ) -> "Array": """Compute the outcome-weighted energy score based on a finite ensemble.""" B = backends.active if backend is None else backends[backend] - M = fct.shape[-2] - wbar = B.sum(fw, -1) / M + wbar = B.mean(fw, axis=-1) err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) # (... M) - E_1 = B.sum(err_norm * fw * B.expand_dims(ow, -1), -1) / (M * wbar) # (...) + E_1 = B.mean(err_norm * fw * B.expand_dims(ow, -1), axis=-1) / wbar # (...) spread_norm = B.norm( B.expand_dims(fct, -2) - B.expand_dims(fct, -3), -1 ) # (... M M) fw_prod = B.expand_dims(fw, -1) * B.expand_dims(fw, -2) # (... M M) spread_norm *= fw_prod * B.expand_dims(ow, (-2, -1)) # (... M M) - E_2 = B.sum(spread_norm, (-2, -1)) / (M**2 * wbar**2) # (...) + E_2 = B.mean(spread_norm, axis=(-2, -1)) / (wbar**2) # (...) return E_1 - 0.5 * E_2 @@ -57,19 +56,18 @@ def vres_ensemble( ) -> "Array": """Compute the vertically re-scaled energy score based on a finite ensemble.""" B = backends.active if backend is None else backends[backend] - M, D = fct.shape[-2:] - wbar = B.sum(fw, -1) / M + wbar = B.mean(fw, axis=-1) err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) # (... M) err_norm *= fw * B.expand_dims(ow, -1) # (... M) - E_1 = B.sum(err_norm, -1) / M # (...) + E_1 = B.mean(err_norm, axis=-1) # (...) spread_norm = B.norm( B.expand_dims(fct, -2) - B.expand_dims(fct, -3), -1 ) # (... M M) fw_prod = B.expand_dims(fw, -2) * B.expand_dims(fw, -1) # (... M M) - E_2 = B.sum(spread_norm * fw_prod, (-2, -1)) / (M**2) # (...) + E_2 = B.mean(spread_norm * fw_prod, axis=(-2, -1)) # (...) - rhobar = B.sum(B.norm(fct, -1) * fw, -1) / M # (...) + rhobar = B.mean(B.norm(fct, -1) * fw, axis=-1) # (...) E_3 = (rhobar - B.norm(obs, -1) * ow) * (wbar - ow) # (...) return E_1 - 0.5 * E_2 + E_3 diff --git a/scoringrules/core/energy/_score_w.py b/scoringrules/core/energy/_score_w.py new file mode 100644 index 0000000..3cd1759 --- /dev/null +++ b/scoringrules/core/energy/_score_w.py @@ -0,0 +1,82 @@ +import typing as tp + +from scoringrules.backend import backends + +if tp.TYPE_CHECKING: + from scoringrules.core.typing import Array, Backend + + +def es_ensemble_w(obs: "Array", fct: "Array", ens_w: "Array", backend=None) -> "Array": + """ + Compute the energy score based on a finite ensemble. + + The ensemble and variables axes are on the second last and last dimensions respectively. + """ + B = backends.active if backend is None else backends[backend] + + err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) + E_1 = B.sum(err_norm * ens_w, -1) + + spread_norm = B.norm(B.expand_dims(fct, -3) - B.expand_dims(fct, -2), -1) + E_2 = B.sum( + spread_norm * B.expand_dims(ens_w, -1) * B.expand_dims(ens_w, -2), (-2, -1) + ) + + return E_1 - 0.5 * E_2 + + +def owes_ensemble_w( + obs: "Array", # (... D) + fct: "Array", # (... M D) + ow: "Array", # (...) + fw: "Array", # (... M) + ens_w: "Array", # (... M) + backend: "Backend" = None, +) -> "Array": + """Compute the outcome-weighted energy score based on a finite ensemble.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(fw * ens_w, -1) + + err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) # (... M) + E_1 = B.sum(err_norm * fw * B.expand_dims(ow, -1) * ens_w, -1) / wbar # (...) + + spread_norm = B.norm( + B.expand_dims(fct, -2) - B.expand_dims(fct, -3), -1 + ) # (... M M) + fw_prod = B.expand_dims(fw, -1) * B.expand_dims(fw, -2) # (... M M) + spread_norm *= fw_prod * B.expand_dims(ow, (-2, -1)) # (... M M) + E_2 = B.sum( + spread_norm * B.expand_dims(ens_w, -1) * B.expand_dims(ens_w, -2), (-2, -1) + ) / (wbar**2) # (...) + + return E_1 - 0.5 * E_2 + + +def vres_ensemble_w( + obs: "Array", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Compute the vertically re-scaled energy score based on a finite ensemble.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(fw * ens_w, -1) + + err_norm = B.norm(fct - B.expand_dims(obs, -2), -1) # (... M) + err_norm *= fw * B.expand_dims(ow, -1) # (... M) + E_1 = B.sum(err_norm * ens_w, -1) # (...) + + spread_norm = B.norm( + B.expand_dims(fct, -2) - B.expand_dims(fct, -3), -1 + ) # (... M M) + fw_prod = B.expand_dims(fw, -2) * B.expand_dims(fw, -1) # (... M M) + E_2 = B.sum( + spread_norm * fw_prod * B.expand_dims(ens_w, -1) * B.expand_dims(ens_w, -2), + (-2, -1), + ) # (...) + + rhobar = B.sum(B.norm(fct, -1) * fw * ens_w, -1) # (...) + E_3 = (rhobar - B.norm(obs, -1) * ow) * (wbar - ow) # (...) + return E_1 - 0.5 * E_2 + E_3 diff --git a/scoringrules/core/kernels/__init__.py b/scoringrules/core/kernels/__init__.py index d0377b7..f7ac807 100644 --- a/scoringrules/core/kernels/__init__.py +++ b/scoringrules/core/kernels/__init__.py @@ -7,15 +7,23 @@ vr_ensemble_mv, ) -try: - from ._gufuncs import estimator_gufuncs -except ImportError: - estimator_gufuncs = None +from ._approx_w import ( + ensemble_uv_w, + ow_ensemble_uv_w, + vr_ensemble_uv_w, + ensemble_mv_w, + ow_ensemble_mv_w, + vr_ensemble_mv_w, +) try: - from ._gufuncs import estimator_gufuncs_mv + from ._gufuncs import estimator_gufuncs_uv, estimator_gufuncs_mv + from ._gufuncs_w import estimator_gufuncs_uv_w, estimator_gufuncs_mv_w except ImportError: + estimator_gufuncs_uv = None estimator_gufuncs_mv = None + estimator_gufuncs_uv_w = None + estimator_gufuncs_mv_w = None __all__ = [ "ensemble_uv", @@ -24,6 +32,14 @@ "ensemble_mv", "ow_ensemble_mv", "vr_ensemble_mv", - "estimator_gufuncs", + "ensemble_uv_w", + "ow_ensemble_uv_w", + "vr_ensemble_uv_w", + "ensemble_mv_w", + "ow_ensemble_mv_w", + "vr_ensemble_mv_w", + "estimator_gufuncs_uv", "estimator_gufuncs_mv", + "estimator_gufuncs_uv_w", + "estimator_gufuncs_mv_w", ] diff --git a/scoringrules/core/kernels/_approx_w.py b/scoringrules/core/kernels/_approx_w.py new file mode 100644 index 0000000..077541d --- /dev/null +++ b/scoringrules/core/kernels/_approx_w.py @@ -0,0 +1,207 @@ +import typing as tp + +from scoringrules.backend import backends + +if tp.TYPE_CHECKING: + from scoringrules.core.typing import Array, ArrayLike, Backend + + +def gauss_kern_uv( + x1: "ArrayLike", x2: "ArrayLike", backend: "Backend" = None +) -> "Array": + """Compute the gaussian kernel evaluated at x1 and x2.""" + B = backends.active if backend is None else backends[backend] + return B.exp(-0.5 * (x1 - x2) ** 2) + + +def gauss_kern_mv( + x1: "ArrayLike", x2: "ArrayLike", backend: "Backend" = None +) -> "Array": + """Compute the gaussian kernel evaluated at vectors x1 and x2.""" + B = backends.active if backend is None else backends[backend] + return B.exp(-0.5 * B.norm(x1 - x2, -1) ** 2) + + +def ensemble_uv_w( + obs: "ArrayLike", + fct: "Array", + ens_w: "Array", + estimator: str = "nrg", + backend: "Backend" = None, +) -> "Array": + """Compute a kernel score for a finite ensemble.""" + B = backends.active if backend is None else backends[backend] + + e_1 = B.sum(gauss_kern_uv(obs[..., None], fct, backend=backend) * ens_w, axis=-1) + e_2 = B.sum( + gauss_kern_uv(fct[..., None], fct[..., None, :], backend=backend) + * ens_w[..., None] + * ens_w[..., None, :], + axis=(-1, -2), + ) + e_3 = gauss_kern_uv(obs, obs) + + if estimator == "fair": + e_2 = e_2 / (1 - B.sum(ens_w * ens_w, axis=-1)) + + out = e_1 - 0.5 * e_2 - 0.5 * e_3 + return -out + + +def ensemble_mv_w( + obs: "ArrayLike", + fct: "Array", + ens_w: "Array", + estimator: str = "nrg", + backend: "Backend" = None, +) -> "Array": + """Compute a kernel score for a finite multivariate ensemble. + + The ensemble and variables axes are on the second last and last dimensions respectively. + """ + B = backends.active if backend is None else backends[backend] + + e_1 = B.sum( + ens_w * gauss_kern_mv(B.expand_dims(obs, -2), fct, backend=backend), axis=-1 + ) + e_2 = B.sum( + gauss_kern_mv(B.expand_dims(fct, -3), B.expand_dims(fct, -2), backend=backend) + * B.expand_dims(ens_w, -1) + * B.expand_dims(ens_w, -2), + axis=(-2, -1), + ) + e_3 = gauss_kern_mv(obs, obs) + + if estimator == "fair": + e_2 = e_2 / (1 - B.sum(ens_w * ens_w, axis=-1)) + + out = e_1 - 0.5 * e_2 - 0.5 * e_3 + return -out + + +def ow_ensemble_uv_w( + obs: "ArrayLike", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Compute an outcome-weighted kernel score for a finite univariate ensemble.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(ens_w * fw, axis=-1) + e_1 = ( + B.sum(ens_w * gauss_kern_uv(obs[..., None], fct, backend=backend) * fw, axis=-1) + * ow + / wbar + ) + e_2 = B.sum( + ens_w[..., None] + * ens_w[..., None, :] + * gauss_kern_uv(fct[..., None], fct[..., None, :], backend=backend) + * fw[..., None] + * fw[..., None, :], + axis=(-1, -2), + ) + e_2 *= ow / (wbar**2) + e_3 = gauss_kern_uv(obs, obs, backend=backend) * ow + + out = e_1 - 0.5 * e_2 - 0.5 * e_3 + return -out + + +def ow_ensemble_mv_w( + obs: "ArrayLike", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Compute an outcome-weighted kernel score for a finite multivariate ensemble. + + The ensemble and variables axes are on the second last and last dimensions respectively. + """ + B = backends.active if backend is None else backends[backend] + wbar = B.sum(fw * ens_w, -1) + + err_kern = gauss_kern_mv(B.expand_dims(obs, -2), fct, backend=backend) + E_1 = B.sum(err_kern * fw * B.expand_dims(ow, -1) * ens_w, axis=-1) / wbar + + spread_kern = gauss_kern_mv( + B.expand_dims(fct, -3), B.expand_dims(fct, -2), backend=backend + ) + fw_prod = B.expand_dims(fw, -1) * B.expand_dims(fw, -2) + spread_kern *= fw_prod * B.expand_dims(ow, (-2, -1)) + E_2 = B.sum( + spread_kern * B.expand_dims(ens_w, -1) * B.expand_dims(ens_w, -2), (-2, -1) + ) / (wbar**2) + + E_3 = gauss_kern_mv(obs, obs, backend=backend) * ow + + out = E_1 - 0.5 * E_2 - 0.5 * E_3 + out = -out + return out + + +def vr_ensemble_uv_w( + obs: "ArrayLike", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Compute a vertically re-scaled kernel score for a finite univariate ensemble.""" + B = backends.active if backend is None else backends[backend] + + e_1 = ( + B.sum(ens_w * gauss_kern_uv(obs[..., None], fct, backend=backend) * fw, axis=-1) + * ow + ) + e_2 = B.sum( + ens_w[..., None] + * ens_w[..., None, :] + * gauss_kern_uv(fct[..., None], fct[..., None, :], backend=backend) + * fw[..., None] + * fw[..., None, :], + axis=(-1, -2), + ) + e_3 = gauss_kern_uv(obs, obs, backend=backend) * ow * ow + + out = e_1 - 0.5 * e_2 - 0.5 * e_3 + out = -out + return out + + +def vr_ensemble_mv_w( + obs: "ArrayLike", + fct: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + backend: "Backend" = None, +) -> "Array": + """Compute a vertically re-scaled kernel score for a finite multivariate ensemble. + + The ensemble and variables axes are on the second last and last dimensions respectively. + """ + B = backends.active if backend is None else backends[backend] + + err_kern = gauss_kern_mv(B.expand_dims(obs, -2), fct, backend=backend) + E_1 = B.sum(err_kern * fw * B.expand_dims(ow, -1) * ens_w, axis=-1) + + spread_kern = gauss_kern_mv( + B.expand_dims(fct, -3), B.expand_dims(fct, -2), backend=backend + ) + fw_prod = B.expand_dims(fw, -1) * B.expand_dims(fw, -2) + spread_kern *= fw_prod + E_2 = B.sum( + spread_kern * B.expand_dims(ens_w, -1) * B.expand_dims(ens_w, -2), (-2, -1) + ) + + E_3 = gauss_kern_mv(obs, obs, backend=backend) * ow * ow + + out = E_1 - 0.5 * E_2 - 0.5 * E_3 + out = -out + return out diff --git a/scoringrules/core/kernels/_gufuncs.py b/scoringrules/core/kernels/_gufuncs.py index 91c8c86..24a8633 100644 --- a/scoringrules/core/kernels/_gufuncs.py +++ b/scoringrules/core/kernels/_gufuncs.py @@ -20,7 +20,6 @@ def _gauss_kern_mv(x1: float, x2: float) -> float: return out -@lazy_gufunc_wrapper_uv @guvectorize("(),(n)->()") def _ks_ensemble_uv_nrg_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray): """Standard version of the kernel score.""" @@ -33,16 +32,15 @@ def _ks_ensemble_uv_nrg_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray e_1 = 0 e_2 = 0 - for x_i in fct: - e_1 += _gauss_kern_uv(x_i, obs) - for x_j in fct: - e_2 += _gauss_kern_uv(x_i, x_j) + for i in range(M): + e_1 += _gauss_kern_uv(fct[i], obs) + for j in range(M): + e_2 += _gauss_kern_uv(fct[i], fct[j]) e_3 = _gauss_kern_uv(obs, obs) - out[0] = -(e_1 / M - 0.5 * e_2 / (M**2) - 0.5 * e_3) + out[0] = -((e_1 / M) - 0.5 * (e_2 / (M**2)) - 0.5 * e_3) -@lazy_gufunc_wrapper_uv @guvectorize("(),(n)->()") def _ks_ensemble_uv_fair_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray): """Fair version of the kernel score.""" @@ -55,16 +53,15 @@ def _ks_ensemble_uv_fair_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarra e_1 = 0 e_2 = 0 - for x_i in fct: - e_1 += _gauss_kern_uv(x_i, obs) - for x_j in fct: - e_2 += _gauss_kern_uv(x_i, x_j) + for i in range(M): + e_1 += _gauss_kern_uv(fct[i], obs) + for j in range(i + 1, M): + e_2 += _gauss_kern_uv(fct[i], fct[j]) e_3 = _gauss_kern_uv(obs, obs) - out[0] = -(e_1 / M - 0.5 * e_2 / (M * (M - 1)) - 0.5 * e_3) + out[0] = -((e_1 / M) - e_2 / (M * (M - 1)) - 0.5 * e_3) -@lazy_gufunc_wrapper_uv @guvectorize("(),(n),(),(n)->()") def _owks_ensemble_uv_gufunc( obs: np.ndarray, @@ -91,10 +88,9 @@ def _owks_ensemble_uv_gufunc( wbar = np.mean(fw) - out[0] = -(e_1 / (M * wbar) - 0.5 * e_2 / ((M * wbar) ** 2) - 0.5 * e_3) + out[0] = -(e_1 / (M * wbar) - 0.5 * e_2 / (M**2 * wbar**2) - 0.5 * e_3) -@lazy_gufunc_wrapper_uv @guvectorize("(),(n),(),(n)->()") def _vrks_ensemble_uv_gufunc( obs: np.ndarray, @@ -119,10 +115,9 @@ def _vrks_ensemble_uv_gufunc( e_2 += _gauss_kern_uv(x_i, x_j) * fw[i] * fw[j] e_3 = _gauss_kern_uv(obs, obs) * ow * ow - out[0] = -(e_1 / M - 0.5 * e_2 / (M**2) - 0.5 * e_3) + out[0] = -((e_1 / M) - 0.5 * (e_2 / (M**2)) - 0.5 * e_3) -@lazy_gufunc_wrapper_mv @guvectorize("(d),(m,d)->()") def _ks_ensemble_mv_nrg_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray): """Standard version of the multivariate kernel score.""" @@ -136,10 +131,9 @@ def _ks_ensemble_mv_nrg_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray e_2 += float(_gauss_kern_mv(fct[i], fct[j])) e_3 = float(_gauss_kern_mv(obs, obs)) - out[0] = -(e_1 / M - 0.5 * e_2 / (M**2) - 0.5 * e_3) + out[0] = -((e_1 / M) - 0.5 * (e_2 / (M**2)) - 0.5 * e_3) -@lazy_gufunc_wrapper_mv @guvectorize("(d),(m,d)->()") def _ks_ensemble_mv_fair_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarray): """Fair version of the multivariate kernel score.""" @@ -149,14 +143,13 @@ def _ks_ensemble_mv_fair_gufunc(obs: np.ndarray, fct: np.ndarray, out: np.ndarra e_2 = 0.0 for i in range(M): e_1 += float(_gauss_kern_mv(fct[i], obs)) - for j in range(M): + for j in range(i + 1, M): e_2 += float(_gauss_kern_mv(fct[i], fct[j])) e_3 = float(_gauss_kern_mv(obs, obs)) - out[0] = -(e_1 / M - 0.5 * e_2 / (M * (M - 1)) - 0.5 * e_3) + out[0] = -((e_1 / M) - (e_2 / (M * (M - 1))) - 0.5 * e_3) -@lazy_gufunc_wrapper_mv @guvectorize("(d),(m,d),(),(m)->()") def _owks_ensemble_mv_gufunc( obs: np.ndarray, @@ -181,7 +174,6 @@ def _owks_ensemble_mv_gufunc( out[0] = -(e_1 / (M * wbar) - 0.5 * e_2 / (M**2 * wbar**2) - 0.5 * e_3) -@lazy_gufunc_wrapper_mv @guvectorize("(d),(m,d),(),(m)->()") def _vrks_ensemble_mv_gufunc( obs: np.ndarray, @@ -201,21 +193,21 @@ def _vrks_ensemble_mv_gufunc( e_2 += float(_gauss_kern_mv(fct[i], fct[j]) * fw[i] * fw[j]) e_3 = float(_gauss_kern_mv(obs, obs)) * ow * ow - out[0] = -(e_1 / M - 0.5 * e_2 / (M**2) - 0.5 * e_3) + out[0] = -((e_1 / M) - 0.5 * (e_2 / (M**2)) - 0.5 * e_3) -estimator_gufuncs = { - "fair": _ks_ensemble_uv_fair_gufunc, - "nrg": _ks_ensemble_uv_nrg_gufunc, - "ow": _owks_ensemble_uv_gufunc, - "vr": _vrks_ensemble_uv_gufunc, +estimator_gufuncs_uv = { + "fair": lazy_gufunc_wrapper_uv(_ks_ensemble_uv_fair_gufunc), + "nrg": lazy_gufunc_wrapper_uv(_ks_ensemble_uv_nrg_gufunc), + "ow": lazy_gufunc_wrapper_uv(_owks_ensemble_uv_gufunc), + "vr": lazy_gufunc_wrapper_uv(_vrks_ensemble_uv_gufunc), } estimator_gufuncs_mv = { - "fair": _ks_ensemble_mv_fair_gufunc, - "nrg": _ks_ensemble_mv_nrg_gufunc, - "ow": _owks_ensemble_mv_gufunc, - "vr": _vrks_ensemble_mv_gufunc, + "fair": lazy_gufunc_wrapper_mv(_ks_ensemble_mv_fair_gufunc), + "nrg": lazy_gufunc_wrapper_mv(_ks_ensemble_mv_nrg_gufunc), + "ow": lazy_gufunc_wrapper_mv(_owks_ensemble_mv_gufunc), + "vr": lazy_gufunc_wrapper_mv(_vrks_ensemble_mv_gufunc), } __all__ = [ diff --git a/scoringrules/core/kernels/_gufuncs_w.py b/scoringrules/core/kernels/_gufuncs_w.py new file mode 100644 index 0000000..8c9a0ee --- /dev/null +++ b/scoringrules/core/kernels/_gufuncs_w.py @@ -0,0 +1,235 @@ +import math + +import numpy as np +from numba import njit, guvectorize + +from scoringrules.core.utils import lazy_gufunc_wrapper_uv, lazy_gufunc_wrapper_mv + + +@njit(["float32(float32, float32)", "float64(float64, float64)"]) +def _gauss_kern_uv(x1: float, x2: float) -> float: + """Gaussian kernel evaluated at x1 and x2.""" + out: float = math.exp(-0.5 * (x1 - x2) ** 2) + return out + + +@njit(["float32(float32[:], float32[:])", "float64(float64[:], float64[:])"]) +def _gauss_kern_mv(x1: float, x2: float) -> float: + """Gaussian kernel evaluated at x1 and x2.""" + out: float = math.exp(-0.5 * np.linalg.norm(x1 - x2) ** 2) + return out + + +@guvectorize("(),(n),(n)->()") +def _ks_ensemble_uv_w_nrg_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """Standard version of the kernel score.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0 + e_2 = 0 + + for i in range(M): + e_1 += _gauss_kern_uv(fct[i], obs) * w[i] + for j in range(M): + e_2 += _gauss_kern_uv(fct[i], fct[j]) * w[i] * w[j] + e_3 = _gauss_kern_uv(obs, obs) + + out[0] = -(e_1 - 0.5 * e_2 - 0.5 * e_3) + + +@guvectorize("(),(n),(n)->()") +def _ks_ensemble_uv_w_fair_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """Fair version of the kernel score.""" + M = fct.shape[-1] + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0 + e_2 = 0 + + for i in range(M): + e_1 += _gauss_kern_uv(fct[i], obs) * w[i] + for j in range(i + 1, M): + e_2 += _gauss_kern_uv(fct[i], fct[j]) * w[i] * w[j] + e_3 = _gauss_kern_uv(obs, obs) + + out[0] = -(e_1 - e_2 / (1 - np.sum(w * w)) - 0.5 * e_3) + + +@guvectorize("(),(n),(),(n),(n)->()") +def _owks_ensemble_uv_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + w: np.ndarray, + out: np.ndarray, +): + """Outcome-weighted kernel score for univariate ensembles.""" + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0.0 + e_2 = 0.0 + + for i, x_i in enumerate(fct): + e_1 += _gauss_kern_uv(x_i, obs) * fw[i] * ow * w[i] + for j, x_j in enumerate(fct): + e_2 += _gauss_kern_uv(x_i, x_j) * fw[i] * fw[j] * ow * w[i] * w[j] + e_3 = _gauss_kern_uv(obs, obs) * ow + + wbar = np.sum(fw * w) + + out[0] = -(e_1 / wbar - 0.5 * e_2 / (wbar**2) - 0.5 * e_3) + + +@guvectorize("(),(n),(),(n),(n)->()") +def _vrks_ensemble_uv_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + w: np.ndarray, + out: np.ndarray, +): + """Vertically re-scaled kernel score for univariate ensembles.""" + + if np.isnan(obs): + out[0] = np.nan + return + + e_1 = 0.0 + e_2 = 0.0 + + for i, x_i in enumerate(fct): + e_1 += _gauss_kern_uv(x_i, obs) * fw[i] * ow * w[i] + for j, x_j in enumerate(fct): + e_2 += _gauss_kern_uv(x_i, x_j) * fw[i] * fw[j] * w[i] * w[j] + e_3 = _gauss_kern_uv(obs, obs) * ow * ow + + out[0] = -(e_1 - 0.5 * e_2 - 0.5 * e_3) + + +@guvectorize("(d),(m,d),(m)->()") +def _ks_ensemble_mv_w_nrg_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """Standard version of the multivariate kernel score.""" + M = fct.shape[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(_gauss_kern_mv(fct[i], obs)) * w[i] + for j in range(M): + e_2 += float(_gauss_kern_mv(fct[i], fct[j])) * w[i] * w[j] + e_3 = float(_gauss_kern_mv(obs, obs)) + + out[0] = -(e_1 - 0.5 * e_2 - 0.5 * e_3) + + +@guvectorize("(d),(m,d),(m)->()") +def _ks_ensemble_mv_w_fair_gufunc( + obs: np.ndarray, fct: np.ndarray, w: np.ndarray, out: np.ndarray +): + """Fair version of the multivariate kernel score.""" + M = fct.shape[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(_gauss_kern_mv(fct[i], obs) * w[i]) + for j in range(i + 1, M): + e_2 += float(_gauss_kern_mv(fct[i], fct[j]) * w[i] * w[j]) + e_3 = float(_gauss_kern_mv(obs, obs)) + + out[0] = -(e_1 - e_2 / (1 - np.sum(w * w)) - 0.5 * e_3) + + +@guvectorize("(d),(m,d),(),(m),(m)->()") +def _owks_ensemble_mv_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + w: np.ndarray, + out: np.ndarray, +): + """Outcome-weighted kernel score for multivariate ensembles.""" + M = fct.shape[0] + ow = ow[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(_gauss_kern_mv(fct[i], obs) * fw[i] * ow * w[i]) + for j in range(M): + e_2 += float( + _gauss_kern_mv(fct[i], fct[j]) * fw[i] * fw[j] * ow * w[i] * w[j] + ) + e_3 = float(_gauss_kern_mv(obs, obs)) * ow + + wbar = np.sum(fw * w) + + out[0] = -(e_1 / wbar - 0.5 * e_2 / (wbar**2) - 0.5 * e_3) + + +def _vrks_ensemble_mv_w_gufunc( + obs: np.ndarray, + fct: np.ndarray, + ow: np.ndarray, + fw: np.ndarray, + w: np.ndarray, + out: np.ndarray, +): + """Vertically re-scaled kernel score for multivariate ensembles.""" + M = fct.shape[0] + ow = ow[0] + + e_1 = 0.0 + e_2 = 0.0 + for i in range(M): + e_1 += float(_gauss_kern_mv(fct[i], obs) * fw[i] * ow * w[i]) + for j in range(M): + e_2 += float(_gauss_kern_mv(fct[i], fct[j]) * fw[i] * fw[j] * w[i] * w[j]) + e_3 = float(_gauss_kern_mv(obs, obs)) * ow * ow + + out[0] = -(e_1 - 0.5 * e_2 - 0.5 * e_3) + + +estimator_gufuncs_uv_w = { + "fair": lazy_gufunc_wrapper_uv(_ks_ensemble_uv_w_fair_gufunc), + "nrg": lazy_gufunc_wrapper_uv(_ks_ensemble_uv_w_nrg_gufunc), + "ow": lazy_gufunc_wrapper_uv(_owks_ensemble_uv_w_gufunc), + "vr": lazy_gufunc_wrapper_uv(_vrks_ensemble_uv_w_gufunc), +} + +estimator_gufuncs_mv_w = { + "fair": lazy_gufunc_wrapper_mv(_ks_ensemble_mv_w_fair_gufunc), + "nrg": lazy_gufunc_wrapper_mv(_ks_ensemble_mv_w_nrg_gufunc), + "ow": lazy_gufunc_wrapper_mv(_owks_ensemble_mv_w_gufunc), + "vr": lazy_gufunc_wrapper_mv(_vrks_ensemble_mv_w_gufunc), +} + +__all__ = [ + "_ks_ensemble_uv_w_fair_gufunc", + "_ks_ensemble_uv_w_nrg_gufunc", + "_owks_ensemble_uv_w_gufunc", + "_vrks_ensemble_uv_w_gufunc", + "_ks_ensemble_mv_w_fair_gufunc", + "_ks_ensemble_mv_w_nrg_gufunc", + "_owks_ensemble_mv_w_gufunc", + "_vrks_ensemble_mv_w_gufunc", +] diff --git a/scoringrules/core/variogram/__init__.py b/scoringrules/core/variogram/__init__.py index 1e74527..8b9749b 100644 --- a/scoringrules/core/variogram/__init__.py +++ b/scoringrules/core/variogram/__init__.py @@ -9,15 +9,36 @@ _variogram_score_gufunc = None _vrvariogram_score_gufunc = None +try: + from ._gufuncs_w import ( + _owvariogram_score_gufunc_w, + _variogram_score_gufunc_w, + _vrvariogram_score_gufunc_w, + ) +except ImportError: + _owvariogram_score_gufunc_w = None + _variogram_score_gufunc_w = None + _vrvariogram_score_gufunc_w = None + from ._score import owvs_ensemble as owvs from ._score import vs_ensemble as vs from ._score import vrvs_ensemble as vrvs +from ._score_w import owvs_ensemble_w as owvs_w +from ._score_w import vs_ensemble_w as vs_w +from ._score_w import vrvs_ensemble_w as vrvs_w + __all__ = [ "vs", + "vs_w", "owvs", + "owvs_w", "vrvs", + "vrvs_w", "_variogram_score_gufunc", + "_variogram_score_gufunc_w", "_owvariogram_score_gufunc", + "_owvariogram_score_gufunc_w", "_vrvariogram_score_gufunc", + "_vrvariogram_score_gufunc_w", ] diff --git a/scoringrules/core/variogram/_gufuncs.py b/scoringrules/core/variogram/_gufuncs.py index e3eb1f4..25465db 100644 --- a/scoringrules/core/variogram/_gufuncs.py +++ b/scoringrules/core/variogram/_gufuncs.py @@ -6,12 +6,12 @@ @guvectorize( [ - "void(float32[:], float32[:,:], float32, float32[:])", - "void(float64[:], float64[:,:], float64, float64[:])", + "void(float32[:], float32[:,:], float32[:,:], float32, float32[:])", + "void(float64[:], float64[:,:], float64[:,:], float64, float64[:])", ], - "(d),(m,d),()->()", + "(d),(m,d),(d,d),()->()", ) -def _variogram_score_gufunc(obs, fct, p, out): +def _variogram_score_gufunc(obs, fct, w, p, out): M = fct.shape[-2] D = fct.shape[-1] out[0] = 0.0 @@ -22,12 +22,12 @@ def _variogram_score_gufunc(obs, fct, p, out): vfct += abs(fct[m, i] - fct[m, j]) ** p vfct = vfct / M vobs = abs(obs[i] - obs[j]) ** p - out[0] += (vobs - vfct) ** 2 + out[0] += w[i, j] * (vobs - vfct) ** 2 @lazy_gufunc_wrapper_mv -@guvectorize("(d),(m,d),(),(),(m)->()") -def _owvariogram_score_gufunc(obs, fct, p, ow, fw, out): +@guvectorize("(d),(m,d),(d,d),(),(m),()->()") +def _owvariogram_score_gufunc(obs, fct, w, ow, fw, p, out): M = fct.shape[-2] D = fct.shape[-1] @@ -38,13 +38,13 @@ def _owvariogram_score_gufunc(obs, fct, p, ow, fw, out): for j in range(D): rho1 = abs(fct[k, i] - fct[k, j]) ** p rho2 = abs(obs[i] - obs[j]) ** p - e_1 += (rho1 - rho2) ** 2 * fw[k] * ow + e_1 += w[i, j] * (rho1 - rho2) ** 2 * fw[k] * ow for m in range(k + 1, M): for i in range(D): for j in range(D): rho1 = abs(fct[k, i] - fct[k, j]) ** p rho2 = abs(fct[m, i] - fct[m, j]) ** p - e_2 += 2 * ((rho1 - rho2) ** 2) * fw[k] * fw[m] * ow + e_2 += 2 * w[i, j] * ((rho1 - rho2) ** 2) * fw[k] * fw[m] * ow wbar = np.mean(fw) @@ -52,8 +52,8 @@ def _owvariogram_score_gufunc(obs, fct, p, ow, fw, out): @lazy_gufunc_wrapper_mv -@guvectorize("(d),(m,d),(),(),(m)->()") -def _vrvariogram_score_gufunc(obs, fct, p, ow, fw, out): +@guvectorize("(d),(m,d),(d,d),(),(m),()->()") +def _vrvariogram_score_gufunc(obs, fct, w, ow, fw, p, out): M = fct.shape[-2] D = fct.shape[-1] @@ -65,14 +65,14 @@ def _vrvariogram_score_gufunc(obs, fct, p, ow, fw, out): for j in range(D): rho1 = abs(fct[k, i] - fct[k, j]) ** p rho2 = abs(obs[i] - obs[j]) ** p - e_1 += (rho1 - rho2) ** 2 * fw[k] * ow - e_3_x += (rho1) ** 2 * fw[k] + e_1 += w[i, j] * (rho1 - rho2) ** 2 * fw[k] * ow + e_3_x += w[i, j] * (rho1) ** 2 * fw[k] for m in range(k + 1, M): for i in range(D): for j in range(D): rho1 = abs(fct[k, i] - fct[k, j]) ** p rho2 = abs(fct[m, i] - fct[m, j]) ** p - e_2 += 2 * ((rho1 - rho2) ** 2) * fw[k] * fw[m] + e_2 += 2 * w[i, j] * ((rho1 - rho2) ** 2) * fw[k] * fw[m] e_3_x *= 1 / M wbar = np.mean(fw) @@ -80,6 +80,6 @@ def _vrvariogram_score_gufunc(obs, fct, p, ow, fw, out): for i in range(D): for j in range(D): rho1 = abs(obs[i] - obs[j]) ** p - e_3_y += (rho1) ** 2 * ow + e_3_y += w[i, j] * (rho1) ** 2 * ow out[0] = e_1 / M - 0.5 * e_2 / (M**2) + (e_3_x - e_3_y) * (wbar - ow) diff --git a/scoringrules/core/variogram/_gufuncs_w.py b/scoringrules/core/variogram/_gufuncs_w.py new file mode 100644 index 0000000..71ef326 --- /dev/null +++ b/scoringrules/core/variogram/_gufuncs_w.py @@ -0,0 +1,88 @@ +import numpy as np +from numba import guvectorize + +from scoringrules.core.utils import lazy_gufunc_wrapper_mv + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(d,d),(m),()->()") +def _variogram_score_gufunc_w(obs, fct, w, ens_w, p, out): + M = fct.shape[-2] + D = fct.shape[-1] + out[0] = 0.0 + for i in range(D): + for j in range(D): + vfct = 0.0 + for m in range(M): + vfct += ens_w[m] * abs(fct[m, i] - fct[m, j]) ** p + vobs = abs(obs[i] - obs[j]) ** p + out[0] += w[i, j] * (vobs - vfct) ** 2 + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(d,d),(),(m),(m),()->()") +def _owvariogram_score_gufunc_w(obs, fct, w, ow, fw, ens_w, p, out): + M = fct.shape[-2] + D = fct.shape[-1] + + e_1 = 0.0 + e_2 = 0.0 + for k in range(M): + for i in range(D): + for j in range(D): + rho1 = abs(fct[k, i] - fct[k, j]) ** p + rho2 = abs(obs[i] - obs[j]) ** p + e_1 += w[i, j] * (rho1 - rho2) ** 2 * fw[k] * ow * ens_w[k] + for m in range(k + 1, M): + for i in range(D): + for j in range(D): + rho1 = abs(fct[k, i] - fct[k, j]) ** p + rho2 = abs(fct[m, i] - fct[m, j]) ** p + e_2 += w[i, j] * ( + 2 + * ((rho1 - rho2) ** 2) + * fw[k] + * fw[m] + * ow + * ens_w[k] + * ens_w[m] + ) + + wbar = np.sum(fw * ens_w) + + out[0] = e_1 / wbar - 0.5 * e_2 / (wbar**2) + + +@lazy_gufunc_wrapper_mv +@guvectorize("(d),(m,d),(d,d),(),(m),(m),()->()") +def _vrvariogram_score_gufunc_w(obs, fct, w, ow, fw, ens_w, p, out): + M = fct.shape[-2] + D = fct.shape[-1] + + e_1 = 0.0 + e_2 = 0.0 + e_3_x = 0.0 + for k in range(M): + for i in range(D): + for j in range(D): + rho1 = abs(fct[k, i] - fct[k, j]) ** p + rho2 = abs(obs[i] - obs[j]) ** p + e_1 += w[i, j] * (rho1 - rho2) ** 2 * fw[k] * ow * ens_w[k] + e_3_x += w[i, j] * (rho1) ** 2 * fw[k] * ens_w[k] + for m in range(k + 1, M): + for i in range(D): + for j in range(D): + rho1 = abs(fct[k, i] - fct[k, j]) ** p + rho2 = abs(fct[m, i] - fct[m, j]) ** p + e_2 += w[i, j] * ( + 2 * ((rho1 - rho2) ** 2) * fw[k] * fw[m] * ens_w[k] * ens_w[m] + ) + + wbar = np.sum(fw * ens_w) + e_3_y = 0.0 + for i in range(D): + for j in range(D): + rho1 = abs(obs[i] - obs[j]) ** p + e_3_y += w[i, j] * (rho1) ** 2 * ow + + out[0] = e_1 - 0.5 * e_2 + (e_3_x - e_3_y) * (wbar - ow) diff --git a/scoringrules/core/variogram/_score.py b/scoringrules/core/variogram/_score.py index 8850cb5..e232f94 100644 --- a/scoringrules/core/variogram/_score.py +++ b/scoringrules/core/variogram/_score.py @@ -9,22 +9,23 @@ def vs_ensemble( obs: "Array", # (... D) fct: "Array", # (... M D) + w: "Array", # (..., D D) p: float = 1, backend: "Backend" = None, ) -> "Array": """Compute the Variogram Score for a multivariate finite ensemble.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-2] fct_diff = B.expand_dims(fct, -2) - B.expand_dims(fct, -1) # (... M D D) - vfct = B.sum(B.abs(fct_diff) ** p, axis=-3) / M # (... D D) + vfct = B.mean(B.abs(fct_diff) ** p, axis=-3) # (... D D) obs_diff = B.expand_dims(obs, -2) - B.expand_dims(obs, -1) # (... D D) vobs = B.abs(obs_diff) ** p # (... D D) - return B.sum((vobs - vfct) ** 2, axis=(-2, -1)) # (...) + return B.sum(w * (vobs - vfct) ** 2, axis=(-2, -1)) # (...) def owvs_ensemble( obs: "Array", fct: "Array", + w: "Array", ow: "Array", fw: "Array", p: float = 1, @@ -32,7 +33,6 @@ def owvs_ensemble( ) -> "Array": """Compute the Outcome-Weighted Variogram Score for a multivariate finite ensemble.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-2] wbar = B.mean(fw, axis=-1) fct_diff = B.expand_dims(fct, -2) - B.expand_dims(fct, -1) # (... M D D) @@ -43,16 +43,18 @@ def owvs_ensemble( del obs, fct E_1 = (fct_diff - B.expand_dims(obs_diff, -3)) ** 2 # (... M D D) - E_1 = B.sum(E_1, axis=(-2, -1)) # (... M) - E_1 = B.sum(E_1 * fw * B.expand_dims(ow, -1), axis=-1) / (M * wbar) # (...) + E_1 = B.sum(B.expand_dims(w, -3) * E_1, axis=(-2, -1)) # (... M) + E_1 = B.mean(E_1 * fw * B.expand_dims(ow, -1), axis=-1) / wbar # (...) fct_diff_spread = B.expand_dims(fct_diff, -3) - B.expand_dims( fct_diff, -4 ) # (... M M D D) fw_prod = B.expand_dims(fw, -2) * B.expand_dims(fw, -1) # (... M M) - E_2 = B.sum(fct_diff_spread**2, axis=(-2, -1)) # (... M M) + E_2 = B.sum( + B.expand_dims(w, (-3, -4)) * fct_diff_spread**2, axis=(-2, -1) + ) # (... M M) E_2 *= fw_prod * B.expand_dims(ow, (-2, -1)) # (... M M) - E_2 = B.sum(E_2, axis=(-2, -1)) / (M**2 * wbar**2) # (...) + E_2 = B.mean(E_2, axis=(-2, -1)) / (wbar**2) # (...) return E_1 - 0.5 * E_2 @@ -60,6 +62,7 @@ def owvs_ensemble( def vrvs_ensemble( obs: "Array", fct: "Array", + w: "Array", ow: "Array", fw: "Array", p: float = 1, @@ -67,7 +70,6 @@ def vrvs_ensemble( ) -> "Array": """Compute the Vertically Re-scaled Variogram Score for a multivariate finite ensemble.""" B = backends.active if backend is None else backends[backend] - M: int = fct.shape[-2] wbar = B.mean(fw, axis=-1) fct_diff = ( @@ -76,20 +78,20 @@ def vrvs_ensemble( obs_diff = B.abs(B.expand_dims(obs, -2) - B.expand_dims(obs, -1)) ** p # (... D D) E_1 = (fct_diff - B.expand_dims(obs_diff, axis=-3)) ** 2 # (... M D D) - E_1 = B.sum(E_1, axis=(-2, -1)) # (... M) - E_1 = B.sum(E_1 * fw * B.expand_dims(ow, axis=-1), axis=-1) / M # (...) + E_1 = B.sum(B.expand_dims(w, -3) * E_1, axis=(-2, -1)) # (... M) + E_1 = B.mean(E_1 * fw * B.expand_dims(ow, axis=-1), axis=-1) # (...) E_2 = ( B.expand_dims(fct_diff, -3) - B.expand_dims(fct_diff, -4) ) ** 2 # (... M M D D) - E_2 = B.sum(E_2, axis=(-2, -1)) # (... M M) + E_2 = B.sum(B.expand_dims(w, (-3, -4)) * E_2, axis=(-2, -1)) # (... M M) fw_prod = B.expand_dims(fw, axis=-2) * B.expand_dims(fw, axis=-1) # (... M M) - E_2 = B.sum(E_2 * fw_prod, axis=(-2, -1)) / (M**2) # (...) + E_2 = B.mean(E_2 * fw_prod, axis=(-2, -1)) # (...) - E_3 = B.sum(fct_diff**2, axis=(-2, -1)) # (... M) - E_3 = B.sum(E_3 * fw, axis=-1) / M # (...) - E_3 -= B.sum(obs_diff**2, axis=(-2, -1)) * ow # (...) + E_3 = B.sum(B.expand_dims(w, -3) * fct_diff**2, axis=(-2, -1)) # (... M) + E_3 = B.mean(E_3 * fw, axis=-1) # (...) + E_3 -= B.sum(w * obs_diff**2, axis=(-2, -1)) * ow # (...) E_3 *= wbar - ow # (...) return E_1 - 0.5 * E_2 + E_3 diff --git a/scoringrules/core/variogram/_score_w.py b/scoringrules/core/variogram/_score_w.py new file mode 100644 index 0000000..b8dad99 --- /dev/null +++ b/scoringrules/core/variogram/_score_w.py @@ -0,0 +1,104 @@ +import typing as tp + +from scoringrules.backend import backends + +if tp.TYPE_CHECKING: + from scoringrules.core.typing import Array, Backend + + +def vs_ensemble_w( + obs: "Array", # (... D) + fct: "Array", # (... M D) + w: "Array", # (..., D D) + ens_w: "Array", # (... M) + p: float = 1, + backend: "Backend" = None, +) -> "Array": + """Compute the Variogram Score for a multivariate finite ensemble.""" + B = backends.active if backend is None else backends[backend] + fct_diff = B.expand_dims(fct, -2) - B.expand_dims(fct, -1) # (... M D D) + vfct = B.sum( + B.expand_dims(ens_w, (-1, -2)) * B.abs(fct_diff) ** p, axis=-3 + ) # (... D D) + obs_diff = B.expand_dims(obs, -2) - B.expand_dims(obs, -1) # (... D D) + vobs = B.abs(obs_diff) ** p # (... D D) + return B.sum(w * (vobs - vfct) ** 2, axis=(-2, -1)) # (...) + + +def owvs_ensemble_w( + obs: "Array", + fct: "Array", + w: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + p: float = 1, + backend: "Backend" = None, +) -> "Array": + """Compute the Outcome-Weighted Variogram Score for a multivariate finite ensemble.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(fw * ens_w, axis=-1) + + fct_diff = B.expand_dims(fct, -2) - B.expand_dims(fct, -1) # (... M D D) + fct_diff = B.abs(fct_diff) ** p # (... M D D) + + obs_diff = B.expand_dims(obs, -2) - B.expand_dims(obs, -1) # (... D D) + obs_diff = B.abs(obs_diff) ** p # (... D D) + del obs, fct + + E_1 = (fct_diff - B.expand_dims(obs_diff, -3)) ** 2 # (... M D D) + E_1 = B.sum(B.expand_dims(w, -3) * E_1, axis=(-2, -1)) # (... M) + E_1 = B.sum(E_1 * fw * B.expand_dims(ow, -1) * ens_w, axis=-1) / wbar # (...) + + fct_diff_spread = B.expand_dims(fct_diff, -3) - B.expand_dims( + fct_diff, -4 + ) # (... M M D D) + fw_prod = B.expand_dims(fw, -2) * B.expand_dims(fw, -1) # (... M M) + ew_prod = B.expand_dims(ens_w, -2) * B.expand_dims(ens_w, -1) # (... M M) + E_2 = B.sum( + B.expand_dims(w, (-3, -4)) * fct_diff_spread**2, axis=(-2, -1) + ) # (... M M) + E_2 *= fw_prod * B.expand_dims(ow, (-2, -1)) * ew_prod # (... M M) + E_2 = B.sum(E_2, axis=(-2, -1)) / (wbar**2) # (...) + + return E_1 - 0.5 * E_2 + + +def vrvs_ensemble_w( + obs: "Array", + fct: "Array", + w: "Array", + ow: "Array", + fw: "Array", + ens_w: "Array", + p: float = 1, + backend: "Backend" = None, +) -> "Array": + """Compute the Vertically Re-scaled Variogram Score for a multivariate finite ensemble.""" + B = backends.active if backend is None else backends[backend] + wbar = B.sum(fw * ens_w, axis=-1) + + fct_diff = ( + B.abs(B.expand_dims(fct, -2) - B.expand_dims(fct, -1)) ** p + ) # (... M D D) + obs_diff = B.abs(B.expand_dims(obs, -2) - B.expand_dims(obs, -1)) ** p # (... D D) + + E_1 = (fct_diff - B.expand_dims(obs_diff, axis=-3)) ** 2 # (... M D D) + E_1 = B.sum(B.expand_dims(w, -3) * E_1, axis=(-2, -1)) # (... M) + E_1 = B.sum(E_1 * fw * B.expand_dims(ow, axis=-1) * ens_w, axis=-1) # (...) + + E_2 = ( + B.expand_dims(fct_diff, -3) - B.expand_dims(fct_diff, -4) + ) ** 2 # (... M M D D) + E_2 = B.sum(B.expand_dims(w, (-3, -4)) * E_2, axis=(-2, -1)) # (... M M) + + fw_prod = B.expand_dims(fw, axis=-2) * B.expand_dims(fw, axis=-1) # (... M M) + ew_prod = B.expand_dims(ens_w, -2) * B.expand_dims(ens_w, -1) # (... M M) + E_2 = B.sum(E_2 * fw_prod * ew_prod, axis=(-2, -1)) # (...) + + E_3 = B.sum(B.expand_dims(w, -3) * fct_diff**2, axis=(-2, -1)) # (... M) + E_3 = B.sum(E_3 * fw * ens_w, axis=-1) # (...) + E_3 -= B.sum(w * obs_diff**2, axis=(-2, -1)) * ow # (...) + E_3 *= wbar - ow # (...) + + return E_1 - 0.5 * E_2 + E_3 diff --git a/tests/test_crps.py b/tests/test_crps.py index d61672b..9abee88 100644 --- a/tests/test_crps.py +++ b/tests/test_crps.py @@ -21,7 +21,7 @@ def test_crps_ensemble(estimator, backend): # test exceptions if backend in ["numpy", "jax", "torch", "tensorflow"]: - if estimator not in ["nrg", "fair", "pwm"]: + if estimator == "int": with pytest.raises(ValueError): sr.crps_ensemble(obs, fct, estimator=estimator, backend=backend) return @@ -43,7 +43,7 @@ def test_crps_ensemble(estimator, backend): assert res.shape == (N,) # non-negative values - if estimator not in ["akr", "akr_circperm"]: + if estimator not in ["akr", "akr_circperm", "int"]: res = sr.crps_ensemble(obs, fct, estimator=estimator, backend=backend) res = np.asarray(res) assert not np.any(res < 0.0) @@ -55,6 +55,57 @@ def test_crps_ensemble(estimator, backend): assert not np.any(res - 0.0 > 0.0001) +@pytest.mark.parametrize("backend", BACKENDS) +def test_crps_ensemble_corr(backend): + obs = np.random.randn(N) + mu = obs + np.random.randn(N) * 0.3 + sigma = abs(np.random.randn(N)) * 0.5 + fct = np.random.randn(N, ENSEMBLE_SIZE) * sigma[..., None] + mu[..., None] + + # test equivalence of different estimators + res_nrg = sr.crps_ensemble(obs, fct, estimator="nrg", backend=backend) + res_qd = sr.crps_ensemble(obs, fct, estimator="qd", backend=backend) + res_fair = sr.crps_ensemble(obs, fct, estimator="fair", backend=backend) + res_pwm = sr.crps_ensemble(obs, fct, estimator="pwm", backend=backend) + if backend in ["torch", "jax"]: + assert np.allclose(res_nrg, res_qd, rtol=1e-03) + assert np.allclose(res_fair, res_pwm, rtol=1e-03) + else: + assert np.allclose(res_nrg, res_qd) + assert np.allclose(res_fair, res_pwm) + + w = np.abs(np.random.randn(N, ENSEMBLE_SIZE) * sigma[..., None]) + res_nrg = sr.crps_ensemble(obs, fct, ens_w=w, estimator="nrg", backend=backend) + res_qd = sr.crps_ensemble(obs, fct, ens_w=w, estimator="qd", backend=backend) + if backend in ["torch", "jax"]: + assert np.allclose(res_nrg, res_qd, rtol=1e-03) + else: + assert np.allclose(res_nrg, res_qd) + + # test correctness + obs = -0.6042506 + fct = np.array( + [ + 1.7812118, + 0.5863797, + 0.7038174, + -0.7743998, + -0.2751647, + 1.1863249, + 1.2990966, + -0.3242982, + -0.5968781, + 0.9064937, + ] + ) + res = sr.crps_ensemble(obs, fct, estimator="qd") + assert np.isclose(res, 0.6126602) + + w = np.arange(10) + res = sr.crps_ensemble(obs, fct, ens_w=w, estimator="qd") + assert np.isclose(res, 0.4923673) + + @pytest.mark.parametrize("backend", BACKENDS) def test_crps_quantile(backend): # test shapes @@ -111,7 +162,17 @@ def test_crps_beta(backend): # test exceptions with pytest.raises(ValueError): - sr.crps_beta(0.3, 0.7, 1.1, lower=1.0, upper=0.0, backend=backend) + sr.crps_beta( + 0.3, 0.7, 1.1, lower=1.0, upper=0.0, backend=backend, check_pars=True + ) + return + + with pytest.raises(ValueError): + sr.crps_beta(0.3, -0.7, 1.1, backend=backend, check_pars=True) + return + + with pytest.raises(ValueError): + sr.crps_beta(0.3, 0.7, -1.1, backend=backend, check_pars=True) return # correctness tests @@ -182,17 +243,17 @@ def test_crps_exponential(backend): @pytest.mark.parametrize("backend", BACKENDS) def test_crps_exponentialM(backend): obs, mass, location, scale = 0.3, 0.1, 0.0, 1.0 - res = sr.crps_exponentialM(obs, mass, location, scale, backend=backend) + res = sr.crps_exponentialM(obs, location, scale, mass, backend=backend) expected = 0.2384728 assert np.isclose(res, expected) obs, mass, location, scale = 0.3, 0.1, -2.0, 3.0 - res = sr.crps_exponentialM(obs, mass, location, scale, backend=backend) + res = sr.crps_exponentialM(obs, location, scale, mass, backend=backend) expected = 0.6236187 assert np.isclose(res, expected) obs, mass, location, scale = -1.2, 0.1, -2.0, 3.0 - res = sr.crps_exponentialM(obs, mass, location, scale, backend=backend) + res = sr.crps_exponentialM(obs, location, scale, mass, backend=backend) expected = 0.751013 assert np.isclose(res, expected) diff --git a/tests/test_kernels.py b/tests/test_kernels.py index b02f511..3c63564 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -20,23 +20,23 @@ def test_gksuv(estimator, backend): sigma = abs(np.random.randn(N)) * 0.3 fct = np.random.randn(N, ENSEMBLE_SIZE) * sigma[..., None] + mu[..., None] - # non-negative values - res = sr.gksuv_ensemble(obs, fct, estimator=estimator, backend=backend) - res = np.asarray(res) - assert not np.any(res < 0.0) - - # m_axis keyword - res = sr.gksuv_ensemble( - obs, - np.random.randn(ENSEMBLE_SIZE, N), - m_axis=0, - estimator=estimator, - backend=backend, - ) - res = np.asarray(res) - assert not np.any(res < 0.0) - if estimator == "nrg": + # non-negative values + res = sr.gksuv_ensemble(obs, fct, estimator=estimator, backend=backend) + res = np.asarray(res) + assert not np.any(res < 0.0) + + # m_axis keyword + res = sr.gksuv_ensemble( + obs, + np.random.randn(ENSEMBLE_SIZE, N), + m_axis=0, + estimator=estimator, + backend=backend, + ) + res = np.asarray(res) + assert not np.any(res < 0.0) + # approx zero when perfect forecast perfect_fct = obs[..., None] + np.random.randn(N, ENSEMBLE_SIZE) * 0.00001 res = sr.gksuv_ensemble(obs, perfect_fct, estimator=estimator, backend=backend) @@ -49,13 +49,6 @@ def test_gksuv(estimator, backend): expected = 0.2490516 assert np.isclose(res, expected) - elif estimator == "fair": - # test correctness - obs, fct = 11.6, np.array([9.8, 8.7, 11.9, 12.1, 13.4]) - res = sr.gksuv_ensemble(obs, fct, estimator=estimator, backend=backend) - expected = 0.2987752 - assert np.isclose(res, expected) - # test exceptions with pytest.raises(ValueError): est = "undefined_estimator" @@ -94,16 +87,6 @@ def test_gksmv(estimator, backend): expected = 0.5868737 assert np.isclose(res, expected) - elif estimator == "fair": - # test correctness - obs = np.array([11.6, -23.1]) - fct = np.array( - [[9.8, 8.7, 11.9, 12.1, 13.4], [-24.8, -18.5, -29.9, -18.3, -21.0]] - ).transpose() - res = sr.gksmv_ensemble(obs, fct, estimator=estimator, backend=backend) - expected = 0.6120162 - assert np.isclose(res, expected) - # test exceptions with pytest.raises(ValueError): est = "undefined_estimator" @@ -212,43 +195,6 @@ def v_func2(x): ) np.testing.assert_allclose(res, 0.0089314, rtol=1e-6) - elif estimator == "fair": - res = np.mean( - np.float64( - sr.twgksuv_ensemble( - obs, fct, v_func=v_func1, estimator=estimator, backend=backend - ) - ) - ) - np.testing.assert_allclose(res, 0.130842, rtol=1e-6) - - res = np.mean( - np.float64( - sr.twgksuv_ensemble( - obs, fct, a=-1.0, estimator=estimator, backend=backend - ) - ) - ) - np.testing.assert_allclose(res, 0.130842, rtol=1e-6) - - res = np.mean( - np.float64( - sr.twgksuv_ensemble( - obs, fct, v_func=v_func2, estimator=estimator, backend=backend - ) - ) - ) - np.testing.assert_allclose(res, 0.1283745, rtol=1e-6) - - res = np.mean( - np.float64( - sr.twgksuv_ensemble( - obs, fct, b=1.85, estimator=estimator, backend=backend - ) - ) - ) - np.testing.assert_allclose(res, 0.1283745, rtol=1e-6) - @pytest.mark.parametrize("backend", BACKENDS) def test_twgksmv(backend): @@ -404,8 +350,8 @@ def w_func(x): @pytest.mark.parametrize("backend", BACKENDS) def test_vrgksuv(backend): - if backend == "jax": - pytest.skip("Not implemented in jax backend") + if backend in ["jax", "torch"]: + pytest.skip("Not implemented in torch and jax backends") obs = np.random.randn(N) mu = obs + np.random.randn(N) * 0.1 sigma = abs(np.random.randn(N)) * 0.3 @@ -426,19 +372,19 @@ def test_vrgksuv(backend): # no argument given resw = sr.vrgksuv_ensemble(obs, fct, backend=backend) - np.testing.assert_allclose(res, resw, rtol=1e-5) + np.testing.assert_allclose(res, resw, rtol=1e-6) # a and b resw = sr.vrgksuv_ensemble( obs, fct, a=float("-inf"), b=float("inf"), backend=backend ) - np.testing.assert_allclose(res, resw, rtol=1e-5) + np.testing.assert_allclose(res, resw, rtol=1e-6) # w_func as identity function resw = sr.vrgksuv_ensemble( obs, fct, w_func=lambda x: x * 0.0 + 1.0, backend=backend ) - np.testing.assert_allclose(res, resw, rtol=1e-5) + np.testing.assert_allclose(res, resw, rtol=1e-6) # test correctness fct = np.array( diff --git a/tests/test_wcrps.py b/tests/test_wcrps.py index e39d02d..10d951f 100644 --- a/tests/test_wcrps.py +++ b/tests/test_wcrps.py @@ -10,34 +10,36 @@ @pytest.mark.parametrize("backend", BACKENDS) def test_owcrps_ensemble(backend): - # test exceptions - with pytest.raises(ValueError): - est = "not_nrg" - sr.owcrps_ensemble(1, 1.1, w_func=lambda x: x, estimator=est, backend=backend) - # test shapes obs = np.random.randn(N) - res = sr.owcrps_ensemble(obs, np.random.randn(N, M), w_func=lambda x: x * 0.0 + 1.0) + res = sr.owcrps_ensemble( + obs, np.random.randn(N, M), w_func=lambda x: x * 0.0 + 1.0, backend=backend + ) assert res.shape == (N,) res = sr.owcrps_ensemble( - obs, np.random.randn(M, N), w_func=lambda x: x * 0.0 + 1.0, m_axis=0 + obs, + np.random.randn(M, N), + w_func=lambda x: x * 0.0 + 1.0, + m_axis=0, + backend=backend, ) assert res.shape == (N,) @pytest.mark.parametrize("backend", BACKENDS) def test_vrcrps_ensemble(backend): - # test exceptions - with pytest.raises(ValueError): - est = "not_nrg" - sr.vrcrps_ensemble(1, 1.1, w_func=lambda x: x, estimator=est, backend=backend) - # test shapes obs = np.random.randn(N) - res = sr.vrcrps_ensemble(obs, np.random.randn(N, M), w_func=lambda x: x * 0.0 + 1.0) + res = sr.vrcrps_ensemble( + obs, np.random.randn(N, M), w_func=lambda x: x * 0.0 + 1.0, backend=backend + ) assert res.shape == (N,) res = sr.vrcrps_ensemble( - obs, np.random.randn(M, N), w_func=lambda x: x * 0.0 + 1.0, m_axis=0 + obs, + np.random.randn(M, N), + w_func=lambda x: x * 0.0 + 1.0, + m_axis=0, + backend=backend, ) assert res.shape == (N,) @@ -75,23 +77,21 @@ def test_owcrps_vs_crps(backend): sigma = abs(np.random.randn(N)) * 0.5 fct = np.random.randn(N, M) * sigma[..., None] + mu[..., None] - res = sr.crps_ensemble(obs, fct, backend=backend, estimator="nrg") + res = sr.crps_ensemble(obs, fct, estimator="qd", backend=backend) # no argument given - resw = sr.owcrps_ensemble(obs, fct, estimator="nrg", backend=backend) - np.testing.assert_allclose(res, resw, rtol=1e-5) + resw = sr.owcrps_ensemble(obs, fct, backend=backend) + np.testing.assert_allclose(res, resw, rtol=1e-4) # a and b resw = sr.owcrps_ensemble( - obs, fct, a=float("-inf"), b=float("inf"), estimator="nrg", backend=backend + obs, fct, a=float("-inf"), b=float("inf"), backend=backend ) - np.testing.assert_allclose(res, resw, rtol=1e-5) + np.testing.assert_allclose(res, resw, rtol=1e-4) # w_func as identity function - resw = sr.owcrps_ensemble( - obs, fct, w_func=lambda x: x * 0.0 + 1.0, estimator="nrg", backend=backend - ) - np.testing.assert_allclose(res, resw, rtol=1e-5) + resw = sr.owcrps_ensemble(obs, fct, w_func=lambda x: x * 0.0 + 1.0, backend=backend) + np.testing.assert_allclose(res, resw, rtol=1e-4) @pytest.mark.parametrize("backend", BACKENDS) @@ -104,19 +104,17 @@ def test_vrcrps_vs_crps(backend): res = sr.crps_ensemble(obs, fct, backend=backend, estimator="nrg") # no argument given - resw = sr.vrcrps_ensemble(obs, fct, estimator="nrg", backend=backend) + resw = sr.vrcrps_ensemble(obs, fct, backend=backend) np.testing.assert_allclose(res, resw, rtol=1e-5) # a and b resw = sr.vrcrps_ensemble( - obs, fct, a=float("-inf"), b=float("inf"), estimator="nrg", backend=backend + obs, fct, a=float("-inf"), b=float("inf"), backend=backend ) np.testing.assert_allclose(res, resw, rtol=1e-5) # w_func as identity function - resw = sr.vrcrps_ensemble( - obs, fct, w_func=lambda x: x * 0.0 + 1.0, estimator="nrg", backend=backend - ) + resw = sr.vrcrps_ensemble(obs, fct, w_func=lambda x: x * 0.0 + 1.0, backend=backend) np.testing.assert_allclose(res, resw, rtol=1e-5) diff --git a/tests/test_wenergy.py b/tests/test_wenergy.py index 5dca3fa..4931f9d 100644 --- a/tests/test_wenergy.py +++ b/tests/test_wenergy.py @@ -23,7 +23,7 @@ def test_owes_vs_es(backend): lambda x: backends[backend].mean(x) * 0.0 + 1.0, backend=backend, ) - np.testing.assert_allclose(res, resw, atol=1e-7) + np.testing.assert_allclose(res, resw, atol=1e-6) @pytest.mark.parametrize("backend", BACKENDS) @@ -48,7 +48,7 @@ def test_vres_vs_es(backend): lambda x: backends[backend].mean(x) * 0.0 + 1.0, backend=backend, ) - np.testing.assert_allclose(res, resw, rtol=1e-10) + np.testing.assert_allclose(res, resw, rtol=1e-6) @pytest.mark.parametrize("backend", BACKENDS) @@ -59,13 +59,13 @@ def test_owenergy_score_correctness(backend): obs = np.array([0.2743836, 0.8146400]) def w_func(x): - return backends[backend].all(x > 0.2) + return backends[backend].all(x > 0.2) * 1.0 res = sr.owes_ensemble(obs, fct, w_func, backend=backend) np.testing.assert_allclose(res, 0.2274243, rtol=1e-6) def w_func(x): - return backends[backend].all(x < 1.0) + return backends[backend].all(x < 1.0) * 1.0 res = sr.owes_ensemble(obs, fct, w_func, backend=backend) np.testing.assert_allclose(res, 0.3345418, rtol=1e-6) diff --git a/tests/test_wvariogram.py b/tests/test_wvariogram.py index 688db83..5592a1d 100644 --- a/tests/test_wvariogram.py +++ b/tests/test_wvariogram.py @@ -23,9 +23,7 @@ def test_owvs_vs_vs(backend): lambda x: backends[backend].mean(x) * 0.0 + 1.0, backend=backend, ) - np.testing.assert_allclose( - res, resw, rtol=1e-3 - ) # TODO: not sure why tolerance must be so high + np.testing.assert_allclose(res, resw, rtol=1e-3) @pytest.mark.parametrize("backend", BACKENDS)