Commit feeb9a95 authored by Pablo Iturrieta Rebolledo's avatar Pablo Iturrieta Rebolledo
Browse files

brier

parent 4701971a
import matplotlib.ticker as tick
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import csep
from csep.core import poisson_evaluations as poisson
from csep.utils import plots
import numpy as np
import copy
import scipy.stats
import scipy.special
import pickle
import os
from codes import get_catalogs, get_models, paths
def mod_bessel0(lambd, iter):
I_0 = 0
for i in range(iter):
I_0 += (lambd**2 / 4) ** (i) / \
(scipy.special.factorial(i))**2
return I_0
#
def p_series(lambd, iter = 10):
p2 = np.exp(-2*lambd) * mod_bessel0(lambd, iter)
return p2
def p_series(lambd, iter = 100):
p2 = 0
for i in range(iter):
p2 += scipy.stats.poisson.pmf(i, lambd) ** 2
return p2
def brier_score(Forecast, Catalog, spatial_only=True, binary=True):
if spatial_only:
data = Forecast.spatial_counts()
obs = Catalog.spatial_counts()
else:
data = Forecast.data
obs = Catalog.spatial_magnitude_counts()
brier = 0
if binary:
obs = (obs >= 1).astype(int)
prob_success = 1 - scipy.stats.poisson.cdf(0, data)
for p, o in zip(prob_success.ravel(), obs.ravel()):
if o == 0:
brier += -2 * p**2
else:
brier += -2 * (p - 1)**2
else:
prob_success = scipy.stats.poisson.pmf(obs, data)
brier = 2 * (prob_success) - p_series(data) - 1
brier = np.sum(brier)
for n_dim in obs.shape:
brier /= n_dim
return brier
if __name__ == '__main__':
# use_saved = True
# forecast1 = get_models.ten_years()[1]
# forecast2 = get_models.ten_years()[5]
# forecast_ref = get_models.ten_years()[11]
#
catalog = get_catalogs.filter_10yr_01_2010(get_catalogs.bollettino())
# with open('trial.obj', 'wb') as file_:
# pickle.dump((forecast1, forecast2, forecast_ref, catalog), file_)
# with open('trial.obj', 'rb') as file_:
# forecast1, forecast2, forecast_ref, catalog = pickle.load(file_)
brier = []
name = []
for model in get_models.ten_years():
br = brier_score(model, catalog, binary=False, spatial_only=False)
brier.append(br)
name.append(model.name)
ind = np.argsort(brier)
brier = np.sort(brier)
name = [name[i] for i in ind]
for i, j in zip(brier, name):
print(i, j)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment