code: def threshold_brier_score(observations, forecasts, threshold, issorted=False,
                            axis=-1):
      observations = np.asarray(observations)
      threshold = np.asarray(threshold)
      forecasts = np.asarray(forecasts)
      if axis != -1:
          forecasts = move_axis_to_end(forecasts, axis)
      if forecasts.shape == observations.shape:
          forecasts = forecasts[..., np.newaxis]
      if observations.shape != forecasts.shape[:-1]:
          raise ValueError('observations and forecasts must have matching '
                           'shapes or matching shapes except along `axis=%s`'
                           % axis)
      scalar_threshold = threshold.ndim == 0
      if threshold.ndim > 1:
          raise ValueError('threshold must be scalar or 1-dimensional')
      if threshold.ndim == 1 and not (np.sort(threshold) == threshold).all():
          raise ValueError('1D thresholds must be sorted')
      threshold = threshold.reshape((1,) * observations.ndim + (-1,))
      if not issorted:
          forecasts = np.sort(forecasts, axis=-1)
      result = _threshold_brier_score_core(observations, forecasts, threshold)
      if scalar_threshold:
          result = result.squeeze(axis=-1)
      return result