def _crps_ensemble_vectorized(observations, forecasts, weights=1):
    """
    An alternative but simpler implementation of CRPS for testing purposes

    This implementation is based on the identity:

    .. math::
        CRPS(F, x) = E_F|X - x| - 1/2 * E_F|X - X'|

    where X and X' denote independent random variables drawn from the forecast
    distribution F, and E_F denotes the expectation value under F.

    Hence it has runtime O(n^2) instead of O(n log(n)) where n is the number of
    ensemble members.

    Reference
    ---------
    Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules,
        prediction, and estimation, 2005. University of Washington Department of
        Statistics Technical Report no. 463R.
        https://www.stat.washington.edu/research/reports/2004/tr463R.pdf
    """
    observations = np.asarray(observations)
    forecasts = np.asarray(forecasts)
    weights = np.asarray(weights)
    if weights.ndim > 0:
        weights = np.where(~np.isnan(forecasts), weights, np.nan)
        weights = weights / np.nanmean(weights, axis=-1, keepdims=True)

    if observations.ndim == forecasts.ndim - 1:
        # sum over the last axis
        assert observations.shape == forecasts.shape[:-1]
        observations = observations[..., np.newaxis]
        with suppress_warnings('Mean of empty slice'):
            score = np.nanmean(weights * abs(forecasts - observations), -1)
        # insert new axes along last and second to last forecast dimensions so
        # forecasts_diff expands with the array broadcasting
        forecasts_diff = (np.expand_dims(forecasts, -1) -
                          np.expand_dims(forecasts, -2))
        weights_matrix = (np.expand_dims(weights, -1) *
                          np.expand_dims(weights, -2))
        with suppress_warnings('Mean of empty slice'):
            score += -0.5 * np.nanmean(weights_matrix * abs(forecasts_diff),
                                       axis=(-2, -1))
        return score
    elif observations.ndim == forecasts.ndim:
        # there is no 'realization' axis to sum over (this is a deterministic
        # forecast)
        return abs(observations - forecasts)
def _crps_ensemble_vectorized ( observations , forecasts , weights = 1 ) : observations = np . asarray ( observations ) forecasts = np . asarray ( forecasts ) weights = np . asarray ( weights ) if weights . ndim > 0 : weights = np . where ( ~ np . isnan ( forecasts ) , weights , np . nan ) weights = weights / np . nanmean ( weights , axis = - 1 , keepdims = True ) if observations . ndim == forecasts . ndim - 1 : assert observations . shape == forecasts . shape [ : - 1 ] observations = observations [ ... , np . newaxis ] with suppress_warnings ( 'Mean of empty slice' ) : score = np . nanmean ( weights * abs ( forecasts - observations ) , - 1 ) forecasts_diff = ( np . expand_dims ( forecasts , - 1 ) - np . expand_dims ( forecasts , - 2 ) ) weights_matrix = ( np . expand_dims ( weights , - 1 ) * np . expand_dims ( weights , - 2 ) ) with suppress_warnings ( 'Mean of empty slice' ) : score += - 0.5 * np . nanmean ( weights_matrix * abs ( forecasts_diff ) , axis = ( - 2 , - 1 ) ) return score elif observations . ndim == forecasts . ndim : return abs ( observations - forecasts )
An alternative but simpler implementation of CRPS for testing purposes

    This implementation is based on the identity:

    .. math::
        CRPS(F, x) = E_F|X - x| - 1/2 * E_F|X - X'|

    where X and X' denote independent random variables drawn from the forecast
    distribution F, and E_F denotes the expectation value under F.

    Hence it has runtime O(n^2) instead of O(n log(n)) where n is the number of
    ensemble members.

    Reference
    ---------
    Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules,
        prediction, and estimation, 2005. University of Washington Department of
        Statistics Technical Report no. 463R.
        https://www.stat.washington.edu/research/reports/2004/tr463R.pdf