docs: 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)