crps_kernel

hydrostats.ens_metrics.crps_kernel(obs, fcst_ens, remove_neg=False, remove_zero=False)

Compute the kernel representation of the continuous ranked probability score (CRPS).

Calculates the kernel representation of the continuous ranked probability score (CRPS) as per equation 3 in Leutbecher et al. (2018) and the adjusted (or fair) crps as per equation 6 in the same paper. Note that it was Gneiting and Raftery (2007) who show the kernel representation as calculated here is equivalent to the standard definition based on the integral over the squared error of the cumulative distribution.

It is strongly recommended to use the hydrostats.ens_metric.ens_crps() function for the improved performance, instead of using this particular implementation, this is meant more for a proof of concept and an algorithm implementation.

Parameters:
  • obs (1D ndarray) – array of observations for each start date.
  • fcst_ens (2D ndarray) – Array of ensemble forecast of dimension n x M, where n = number of start dates and M = number of ensemble members.
  • remove_neg (bool) – If True, when a negative value is found at the i-th position in the observed OR ensemble array, the i-th value of the observed AND ensemble array are removed before the computation.
  • remove_zero (bool) – If true, when a zero value is found at the i-th position in the observed OR ensemble array, the i-th value of the observed AND ensemble array are removed before the computation.
Returns:

Dictionary of outputs includes:
  • [“crps”] 1D ndarray with crps values of length n.
  • [“crpsAdjusted”] 1D ndarray with adjusted crps values of length n.
  • [“crpsMean”] Arithmetic mean of crps values as a float.
  • [“crpsAdjustedMean”] Arithmetic mean of adjusted crps values as a float.
Return type:

dict

Notes

NaN treatment: If any start date in obs is NaN, then the corresponding row in fcst (for all ensemble members) will also be deleted.

Examples

>>> import numpy as np
>>> import hydrostats.ens_metrics as em
>>> np.random.seed(3849590438)

Creating an observed 1D array and an ensemble 2D array with all random numbers

>>> ens_array_random = (np.random.rand(15, 52) + 1) * 100
>>> obs_array_random = (np.random.rand(15) + 1) * 100

Creating an observed 1D array and an ensemble 2D array with noise.

>>> noise = np.random.normal(scale=1, size=(15, 52))
>>> x = np.linspace(1, 10, 15)
>>> observed_array = np.sin(x) + 10
>>> ensemble_array_noise = (np.ones((15, 52)).T * observed_array).T + noise

Computing the Hersbach CRPS values between the ensemble mean and the observed data with the random data.

>>> crps_dictionary_rand = em.crps_kernel(obs_array_random, ens_array_random)
>>> print(crps_dictionary_rand['crps'])
[ 7.73360237  9.59248626 34.46719655 30.10271075  7.451665   16.07882352
 14.59543529  8.55181637 15.4833089   8.32422363 16.55108154 19.20821296
  8.39452279 12.59949378 27.82543302]
>>> print(crps_dictionary_rand['crpsAdjusted'])
[ 7.43000827  9.29100065 34.14067524 29.76359191  7.14776152 15.75147589
 14.25192856  8.23647876 15.19419171  8.05998301 16.26113448 18.90686679
  8.09725139 12.24021268 27.45673444]
>>> crps_dictionary_rand['crpsMean']
15.797334183277723
>>> crps_dictionary_rand['crpsAdjustedMean']
15.481953018707593

Computing the Hersbach CRPS values between the ensemble mean and the observed data with noise in the ensemble data.

>>> crps_dictionary_noise = em.crps_kernel(obs=observed_array, fcst_ens=ensemble_array_noise)
>>> print(crps_dictionary_noise['crps'])
[0.26921152 0.21388687 0.24927151 0.26047667 0.30234843 0.1996493
 0.2779844  0.29478927 0.275383   0.25682693 0.21485236 0.22824711
 0.2813889  0.21264652 0.18141063]
>>> print(crps_dictionary_noise['crpsAdjusted'])
[0.25850726 0.20482797 0.23908004 0.25032814 0.28996894 0.18961905
 0.2670867  0.2821429  0.2634554  0.24573507 0.20457832 0.21730326
 0.26951946 0.2034818  0.17198615]
>>> crps_dictionary_noise['crpsMean']
0.2478915604121471
>>> crps_dictionary_noise['crpsAdjustedMean']
0.23717469718824744

References

  • Gneiting, T. and Raftery, A. E. (2007) Strictly proper scoring rules, prediction and estimation, J. Amer. Stat. Asoc., 102, 359-378.
  • Leutbecher, M. (2018) Ensemble size: How suboptimal is less than infinity?, Q. J. R. Meteorol., Accepted.