Merge pull request #526 from nci/release_9_0

Release 0.9.0

README.md

@@ -1,8 +1,5 @@
 # Scores: Verification and Evaluation for Forecasts and Models
 
-> **Notice:**
-> **`scores` is undergoing final testing and review. When this is completed, this notice will be removed.**
-
 > 
 > **A list of over 50 metrics, statistical techniques and data processing tools contained in `scores` is [available here](https://scores.readthedocs.io/en/stable/included.html).**
 

docs/conf.py

@@ -9,7 +9,7 @@
 
 project = "scores"
 copyright = "Licensed under Apache 2.0 - https://www.apache.org/licenses/LICENSE-2.0"
-release = "0.8.6"
+release = "0.9.0"
 
 version = __version__
 

docs/included.md

@@ -391,14 +391,6 @@
     [Tutorial](project:./tutorials/Binary_Contingency_Scores.md)
   - 
     [Wikipedia](https://en.wikipedia.org/wiki/Precision_and_recall); [Probability of detection (WWRP/WGNE Joint Working Group on Forecast Verification Research)](https://www.cawcr.gov.au/projects/verification/#POD)
-* -  
-    - Symmetric Extremal Dependence Index (SEDI)
-  -
-    [API](api.md#scores.categorical.BasicContingencyManager.symmetric_extremal_dependence_index)
-  - 
-    [Tutorial](project:./tutorials/Binary_Contingency_Scores.md)
-  - 
-    [Ferro and Stephenson (2011)](https://doi.org/10.1175/WAF-D-10-05030.1)
 * -  
     - Sensitivity (Hit Rate, Probability of Detection (POD), True Positive Rate, Recall)
   -
@@ -423,6 +415,14 @@
     [Tutorial](project:./tutorials/Binary_Contingency_Scores.md)
   - 
     [Success ratio (WWRP/WGNE Joint Working Group on Forecast Verification Research)](https://www.cawcr.gov.au/projects/verification/#SR)    
+* -  
+    - Symmetric Extremal Dependence Index (SEDI)
+  -
+    [API](api.md#scores.categorical.BasicContingencyManager.symmetric_extremal_dependence_index)
+  - 
+    [Tutorial](project:./tutorials/Binary_Contingency_Scores.md)
+  - 
+    [Ferro and Stephenson (2011)](https://doi.org/10.1175/WAF-D-10-05030.1)
 * -  
     - Threat Score (Critical Success Index)
   -

pyproject.toml

@@ -48,7 +48,8 @@ tutorial = [
     "h5netcdf",
     "rasterio",
     "rioxarray",
-    "plotly"
+    "plotly",
+    "dask"
 ]
 maintainer = ["build",
               "hatch",

src/scores/__init__.py

@@ -13,7 +13,7 @@
 import scores.sample_data
 import scores.stats.statistical_tests  # noqa: F401
 
-__version__ = "0.8.6"
+__version__ = "0.9.0"
 
 __all__ = [
     "scores.categorical",

src/scores/continuous/standard_impl.py

@@ -23,13 +23,11 @@ def mse(
 ) -> XarrayLike:
     """Calculates the mean squared error from forecast and observed data.
 
-    Dimensional reduction is not supported for pandas and the user should
-    convert their data to xarray to formulate the call to the metric. At
-    most one of reduce_dims and preserve_dims may be specified.
-    Specifying both will result in an exception.
-
     See "Mean squared error" section at https://www.cawcr.gov.au/projects/verification/#MSE for more information
 
+    .. math ::
+        \\frac{1}{n} \\sum_{i=1}^n (\\text{forecast}_i - \\text{observed}_i)^2
+
     Args:
         fcst (Union[xr.Dataset, xr.DataArray, pd.Dataframe, pd.Series]):
             Forecast or predicted variables in xarray or pandas.
@@ -211,6 +209,16 @@ def correlation(
     """
     Calculates the Pearson's correlation coefficient between two xarray DataArrays
 
+    .. math::
+        \\rho = \\frac{\\sum_{i=1}^{n}{(x_i - \\bar{x})(y_i - \\bar{y})}}{\\sqrt{\\sum_{i=1}^{n}{(x_i-\\bar{x})^2}\\sum_{i=1}^{n}{(y_i - \\bar{y})^2}}}
+
+    where:
+        - :math:`\\rho` = Pearson's correlation coefficient
+        - :math:`x_i` = the values of x in a sample (i.e. forecast values)
+        - :math:`\\bar{x}` = the mean value of the forecast sample
+        - :math:`y_i` = the values of y in a sample (i.e. observed values)
+        - :math:`\\bar{y}` = the mean value of the observed sample value
+
     Args:
         fcst: Forecast or predicted variables
         obs: Observed variables.
@@ -225,7 +233,7 @@ def correlation(
             point (i.e. single-value comparison against observed), and the
             forecast and observed dimensions must match precisely.
     Returns:
-        An xarray object with Pearson's correlation coefficient values
+        xr.DataArray: An xarray object with Pearson's correlation coefficient values
     """
     reduce_dims = scores.utils.gather_dimensions(
         fcst.dims, obs.dims, reduce_dims=reduce_dims, preserve_dims=preserve_dims

src/scores/pandas/continuous.py

@@ -16,6 +16,9 @@ def mse(
 
     A detailed explanation is on https://en.wikipedia.org/wiki/Mean_squared_error
 
+    .. math ::
+        \\frac{1}{n} \\sum_{i=1}^n (\\text{forecast}_i - \\text{observed}_i)^2
+
 
     Notes:
         Dimensional reduction is not supported for pandas and the user should

src/scores/probability/crps_impl.py

@@ -786,20 +786,20 @@ def crps_for_ensemble(
 
 
     .. math::
-        CRPS(x_i, y) = (1 / M) * sum(|x_i - y|) - (1 / 2 * K) * sum(|x_i - x_j|)
+        CRPS(x_i, y) = \\frac{\\sum_{i=1}^{M}(|x_i - y|)}{M} - \\frac{\\sum_{i=1}^{M}(|x_i - x_j|)}{2K}
 
     where the first sum is iterated over 1 <= i <= M and the second sum is iterated over
     1 <= i <= M and 1 <= j <= M.
 
     The value of the constant K in this formula depends on the method:
-        - If `method="ecdf"` then :math:`K = M ** 2`. In this case the CRPS value returned is \
+        - If `method="ecdf"` then :math:`K = M ^ 2`. In this case the CRPS value returned is \
             the exact CRPS value for the emprical cumulation distribution function \
             constructed using the ensemble values.
-        - If `method="fair"` then :math:`K = M * (M - 1)`. In this case the CRPS value returned \
+        - If `method="fair"` then :math:`K = M(M - 1)`. In this case the CRPS value returned \
             is the approximated CRPS where the ensemble values can be interpreted as a \
             random sample from the underlying predictive distribution. This interpretation \
-            stems from the formula :math:`\\text{CRPS}(F, Y) = E|X - Y| - E|X - X'|/2`, where X and X' \
-            are independent samples of the predictive distribution F, Y is the observation \
+            stems from the formula :math:`\\text{CRPS}(F, y) = \\mathbb{E}(|X - y|) - \\frac{1}{2}\\mathbb{E}(|X - X'|)`, where X and X' \
+            are independent samples of the predictive distribution F, y is the observation \
             (possibly unknown) and E denotes the expectation. This choice of K gives an \
             unbiased estimate for the second expectation.
 
@@ -820,7 +820,7 @@ def crps_for_ensemble(
         ValueError: when method is not one of "ecdf" or "fair".
 
     See also:
-        `scores.probability.crps_cdf`
+        :py:func:`scores.probability.crps_cdf`
 
     References:
         - C. Ferro (2014), "Fair scores for ensemble forecasts", Quarterly Journal of the \

tests/categorical/test_contingency.py

@@ -71,6 +71,15 @@
     dims=["height", "lat", "lon"],
 )
 
+finleys_table = {
+    "tp_count": xr.DataArray(28),
+    "fp_count": xr.DataArray(72),
+    "fn_count": xr.DataArray(23),
+    "tn_count": xr.DataArray(2680),
+    "total_count": xr.DataArray(2803),
+}
+
+
 # This truth table shows where the forecast and obs match within 0.2
 somewhat_near_matches = xr.DataArray(
     [
@@ -368,3 +377,26 @@ def test_dask_if_available_categorical():
 
     # And that transformed things get the same numbers
     assert table.false_alarm_rate() == table.transform().false_alarm_rate()
+
+
+def test_examples_with_finley():
+    """
+    Test some of the complex scores with the Finley tornado example
+    """
+
+    table = finleys_table
+    cm = scores.categorical.BasicContingencyManager(table)
+    heidke = cm.heidke_skill_score()
+    gilbert = cm.gilberts_skill_score()
+
+    # Note - the reference in the verification site has 0.36 for the expected
+    # result, but presumably this is rounded to two decimal places
+    # See https://www.cawcr.gov.au/projects/verification/Finley/Finley_Tornados.html
+    heidke_expected = xr.DataArray(0.355325)
+    xr.testing.assert_allclose(heidke, heidke_expected)
+
+    # Note - the reference in the verification site has 0.22 for the expected
+    # result, but presumably this is rounded to two decimal places
+    # See https://www.cawcr.gov.au/projects/verification/Finley/Finley_Tornados.html
+    gilbert_expected = xr.DataArray(0.216046)
+    xr.testing.assert_allclose(gilbert_expected, gilbert)