Skip to content

Commit

Permalink
ft: finalized brier score implementation as well as brier consistency…
Browse files Browse the repository at this point in the history
… test. It is implemented only using a Binary calculation for spatial-magnitude bins.

tests: Implemented brier score tests.
  • Loading branch information
pabloitu committed Jan 22, 2024
1 parent 5fb55ca commit a50d1e3
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 30 deletions.
56 changes: 43 additions & 13 deletions csep/core/brier_evaluations.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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:
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down
16 changes: 0 additions & 16 deletions tests/test_brier.py

This file was deleted.

61 changes: 60 additions & 1 deletion tests/test_evaluations.py
Original file line number Diff line number Diff line change
@@ -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')
Expand Down Expand Up @@ -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()

0 comments on commit a50d1e3

Please sign in to comment.