ft: finalized brier score implementation as well as brier consistency test. It is implemented only using a Binary calculation for spatial-magnitude bins.

tests: Implemented brier score tests.

Author: pabloitu

csep/core/brier_evaluations.py

@@ -7,8 +7,22 @@
 
 
 def _brier_score_ndarray(forecast, observations):
-    """ Computes the brier score
+    """ Computes the brier (binary) score for spatial-magnitude cells
+    using the formula:
+
+    Q(Lambda, Sigma) = 1/N sum_{i=1}^N (Lambda_i - Ind(Sigma_i > 0 ))^2
+
+    where Lambda is the forecast array, Sigma is the observed catalog, N the
+    number of spatial-magnitude cells and Ind is the indicator function, which
+    is 1 if Sigma_i > 0 and 0 otherwise.
+
+    Args:
+        forecast: 2d array of forecasted rates
+        observations: 2d array of observed counts
+    Returns
+        brier: float, brier score
     """
+
     prob_success = 1 - poisson.cdf(0, forecast)
     brier_cell = np.square(prob_success.ravel() - (observations.ravel() > 0))
     brier = -2 * brier_cell.sum()
@@ -19,6 +33,18 @@ def _brier_score_ndarray(forecast, observations):
 
 
 def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None):
+    """
+    Simulates a catalog by sampling from the sampling_weights array.
+    Identical to binomial_evaluations._simulate_catalog
+
+    Args:
+        sim_cells:
+        sampling_weights:
+        random_numbers:
+
+    Returns:
+
+    """
     sim_fore = numpy.zeros(sampling_weights.shape)
 
     if random_numbers is None:
@@ -46,10 +72,19 @@ def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None):
 
 def _brier_score_test(forecast_data, observed_data, num_simulations=1000,
                       random_numbers=None, seed=None, verbose=True):
-    """  Computes binary conditional-likelihood test from CSEP using an
-    efficient simulation based approach.
+    """  Computes the Brier consistency test conditional on the total observed
+    number of events
 
     Args:
+        forecast_data: 2d array of forecasted rates for spatial_magnitude cells
+        observed_data: 2d array of a catalog resampled to spatial_magnitude
+            cells
+        num_simulations: number of synthetic catalog simulations
+        random_numbers: numpy array of random numbers to use for simulation
+        seed: seed for random number generator
+        verbose: print status updates
+
+
 
     """
     # Array-masking that avoids log singularities:
@@ -58,8 +93,6 @@ def _brier_score_test(forecast_data, observed_data, num_simulations=1000,
     # set seed for the likelihood test
     if seed is not None:
         numpy.random.seed(seed)
-    import time
-    start = time.process_time()
 
     # used to determine where simulated earthquake should
     # be placed, by definition of cumsum these are sorted
@@ -93,7 +126,6 @@ def _brier_score_test(forecast_data, observed_data, num_simulations=1000,
                 print(f'... {idx + 1} catalogs simulated.')
 
     obs_brier = _brier_score_ndarray(forecast_data.data, observed_data)
-
     # quantile score
     qs = numpy.sum(simulated_brier <= obs_brier) / num_simulations
 
@@ -107,13 +139,11 @@ def brier_score_test(gridded_forecast,
                      seed=None,
                      random_numbers=None,
                      verbose=False):
-    """ Performs the Brier conditional test on Gridded Forecast using an
-    Observed Catalog.
-    Normalizes the forecast so the forecasted rate are consistent with the
-     observations. This modification
-    eliminates the strong impact differences in the number distribution
-     have on the forecasted rates.
-
+    """
+    Performs the Brier conditional test on a Gridded Forecast using an
+    Observed Catalog. Normalizes the forecast so the forecasted rate are
+    consistent with the observations. This modification eliminates the strong
+    impact differences in the number distribution have on the forecasted rates.
     """
 
     # grid catalog onto spatial grid

tests/test_brier.py

@@ -1,10 +1,12 @@
 import os
 import numpy
+import numpy as np
+import scipy
 import unittest
 
 import csep.core.poisson_evaluations as poisson
 import csep.core.binomial_evaluations as binary
-
+import csep.core.brier_evaluations as brier
 def get_datadir():
     root_dir = os.path.dirname(os.path.abspath(__file__))
     data_dir = os.path.join(root_dir, 'artifacts', 'Comcat')
@@ -115,5 +117,62 @@ def test_binomial_likelihood(self):
         numpy.testing.assert_allclose(simulated_ll[0], -7.921741654647629)
 
 
+class TestPoissonBrier(unittest.TestCase):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.seed = 1
+        self.forecast_data = numpy.array([[0.3, 0.2, 0.1], [0.2, 0.1, 0.1]])
+        self.observed_data = numpy.array([[0, 1, 2], [1, 1, 0]])
+
+    def test_brier_score_calculation(self):
+
+        # 1 bin
+        rate = 1
+        prob = 1 - scipy.stats.poisson.pmf(0, rate)
+        brier_score_hit = -2 * (prob - 1)**2
+        brier_score = brier._brier_score_ndarray(numpy.array([[rate]]),
+                                                 numpy.array([[1]]))
+        numpy.testing.assert_allclose(brier_score_hit, brier_score)
+
+        brier_score_nohit = -2 * prob**2
+        brier_score = brier._brier_score_ndarray(numpy.array([[rate]]),
+                                                 numpy.array([[0]]))
+        numpy.testing.assert_allclose(brier_score_nohit, brier_score)
+        # 2 bins
+        rate = np.array([1, 0.5])
+        hits = np.array([0, 1])
+        prob = 1 - scipy.stats.poisson.pmf(0, rate)
+        brier_score_2bins = (-2 * (prob - hits)**2).sum() / len(rate)
+        brier_score = brier._brier_score_ndarray(np.array(rate),
+                                                 np.array([hits]))
+        numpy.testing.assert_allclose(brier_score_2bins, brier_score)
+
+        hits = np.array([0, 0])
+        brier_score_2bins = (-2 * (prob - hits)**2).sum() / len(rate)
+        brier_score = brier._brier_score_ndarray(np.array(rate),
+                                                 np.array([hits]))
+        numpy.testing.assert_allclose(brier_score_2bins, brier_score)
+
+    def test_brier_test(self):
+
+        expected_sim = numpy.array([1, 1, 1, 1, 0, 0])
+        qs, obs_brier, simulated_brier = brier._brier_score_test(
+            self.forecast_data,
+            self.observed_data,
+            num_simulations=1,
+            seed=self.seed,
+            verbose=True)
+
+        probs = 1 - scipy.stats.poisson.pmf(0, self.forecast_data.ravel())
+        sim_brier_analytic = ((-2 * (probs - expected_sim)**2).sum() /
+                              len(probs))
+        numpy.testing.assert_allclose(simulated_brier[0], sim_brier_analytic)
+        obs_brier_analytic = (
+            (-2 * (probs - self.observed_data.ravel().astype(bool))**2).sum() /
+            len(probs))
+        numpy.testing.assert_allclose(obs_brier, obs_brier_analytic)
+
+
 if __name__ == '__main__':
     unittest.main()