Skip to content
GitLab

implement quantile scores and CRPS

Merged François Bourgin requested to merge 1-add-probabilistic-scores into main
Compare
and
11 files
+ 290
− 54
include/evalhyd/determinist/evaluator.hpp
+ 4
− 4
@@ -30,11 +30,11 @@ namespace evalhyd
}
};
// members for evaluation metrics
A nse;
A NSE;
// methods to compute elements
// TODO
// methods to compute metrics
void calc_nse();
void calc_NSE();
};
// Compute the Nash-Sutcliffe Efficiency (NSE).
@@ -54,7 +54,7 @@ namespace evalhyd
// Array of computed Nash-Sutcliffe efficiencies.
// shape: (...)
template <class A>
void Evaluator<A>::calc_nse()
void Evaluator<A>::calc_NSE()
{
// compute average observed flow
A q_avg = xt::mean(q_obs, -1, xt::keep_dims);
@@ -64,7 +64,7 @@ namespace evalhyd
A f_den = xt::sum(xt::square(q_obs - q_avg), -1, xt::keep_dims);
// return computed NSE
nse = 1 - (f_num / f_den);
NSE = 1 - (f_num / f_den);
}
}
}
include/evalhyd/probabilist/evaluator.h
+ 27
− 6
@@ -19,6 +19,7 @@ namespace evalhyd
xt::xtensor<double, 2> o_k;
xt::xtensor<double, 1> bar_o;
xt::xtensor<double, 2> y_k;
xt::xtensor<double, 2> q_qnt;
public:
// constructor method
@@ -27,22 +28,42 @@ namespace evalhyd
const xt::xtensor<double, 1>& thr) :
q_obs{obs}, q_frc{frc}, q_thr{thr} {};
// members for evaluation metrics
// members for intermediate evaluation metrics
// (i.e. before the reduction along the temporal axis)
xt::xtensor<double, 2> bs;
xt::xtensor<double, 2> bs_crd;
xt::xtensor<double, 2> bs_lbd;
xt::xtensor<double, 2> bss;
xt::xtensor<double, 2> qs;
xt::xtensor<double, 2> crps;
// members for evaluation metrics
xt::xtensor<double, 2> BS;
xt::xtensor<double, 2> BS_CRD;
xt::xtensor<double, 2> BS_LBD;
xt::xtensor<double, 2> BSS;
xt::xtensor<double, 2> QS;
xt::xtensor<double, 2> CRPS;
// methods to compute derived data
void calc_q_qnt();
// methods to compute elements
void calc_o_k();
void calc_bar_o();
void calc_y_k();
// methods to compute metrics
// methods to compute intermediate metrics
void calc_bs();
void calc_bs_crd();
void calc_bs_lbd();
void calc_bss();
void calc_qs();
void calc_crps();
// methods to compute metrics
void calc_BS();
void calc_BS_CRD();
void calc_BS_LBD();
void calc_BSS();
void calc_QS();
void calc_CRPS();
};
xt::xtensor<double, 3> is_above_threshold(
include/evalhyd/probabilist/evaluator_brier.cpp
+ 53
− 23
@@ -17,7 +17,7 @@ namespace evalhyd
{
namespace probabilist
{
// Compute the Brier score (BS).
// Compute the Brier score for each time step.
//
// \require o_k:
// Observed event outcome.
@@ -26,13 +26,28 @@ namespace evalhyd
// Event probability forecast.
// shape: (thresholds, time)
// \assign bs:
// 2D array of Brier Score for each threshold.
// shape: (thresholds, 1)
// 2D array of Brier score for each threshold for each time step.
// shape: (thresholds, time)
void Evaluator::calc_bs()
{
// return computed Brier score(s)
// $bs = (o_k - y_k)^2$
bs = xt::square(o_k - y_k);
}
// Compute the Brier score (BS).
//
// \require bs:
// 2D array of Brier score for each threshold for each time step.
// shape: (thresholds, time)
// \assign BS:
// 2D array of Brier score for each threshold.
// shape: (thresholds, 1)
void Evaluator::calc_BS()
{
// calculate the mean over the time steps
// $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
bs = xt::mean(xt::square(o_k - y_k), -1, xt::keep_dims);
BS = xt::mean(bs, -1, xt::keep_dims);
}
// Compute the calibration-refinement decomposition of the Brier score
@@ -49,11 +64,11 @@ namespace evalhyd
// \require y_k:
// Event probability forecast.
// shape: (thresholds, time)
// \assign bs_crd:
// \assign BS_CRD:
// 2D array of Brier score components (reliability, resolution,
// uncertainty) for each threshold.
// shape: (thresholds, components)
void Evaluator::calc_bs_crd()
void Evaluator::calc_BS_CRD()
{
// declare internal variables
// shape: (bins, thresholds, time)
@@ -71,7 +86,7 @@ namespace evalhyd
// initialise output variable
// shape: (components, thresholds)
bs_crd = xt::zeros<double>({n_thr, n_cmp});
BS_CRD = xt::zeros<double>({n_thr, n_cmp});
// compute range of forecast values $y_i$
y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr;
@@ -103,7 +118,7 @@ namespace evalhyd
// calculate reliability =
// $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
xt::col(bs_crd, 0) =
xt::col(BS_CRD, 0) =
xt::sum(
xt::square(
xt::view(y_i, xt::all(), xt::newaxis())
@@ -114,7 +129,7 @@ namespace evalhyd
// calculate resolution =
// $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
xt::col(bs_crd, 1) =
xt::col(BS_CRD, 1) =
xt::sum(
xt::square(
bar_o_i - bar_o
@@ -123,7 +138,7 @@ namespace evalhyd
) / n;
// calculate uncertainty = $\bar{o} (1 - \bar{o})$
xt::col(bs_crd, 2) = bar_o * (1 - bar_o);
xt::col(BS_CRD, 2) = bar_o * (1 - bar_o);
}
// Compute the likelihood-base rate decomposition of the Brier score
@@ -137,11 +152,11 @@ namespace evalhyd
// \require y_k:
// Event probability forecast.
// shape: (thresholds, time)
// \return bs_lbd:
// \return BS_LBD:
// 2D array of Brier score components (type 2 bias, discrimination,
// sharpness) for each threshold.
// shape: (thresholds, components)
void Evaluator::calc_bs_lbd()
void Evaluator::calc_BS_LBD()
{
// declare internal variables
// shape: (bins, thresholds, time)
@@ -160,7 +175,7 @@ namespace evalhyd
// declare and initialise output variable
// shape: (components, thresholds)
bs_lbd = xt::zeros<double>({n_thr, n_cmp});
BS_LBD = xt::zeros<double>({n_thr, n_cmp});
// set the range of observed values $o_j$
o_j = {0., 1.};
@@ -196,7 +211,7 @@ namespace evalhyd
// calculate type 2 bias =
// $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
xt::col(bs_lbd, 0) =
xt::col(BS_LBD, 0) =
xt::sum(
xt::square(
xt::view(o_j, xt::all(), xt::newaxis())
@@ -207,7 +222,7 @@ namespace evalhyd
// calculate discrimination =
// $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
xt::col(bs_lbd, 1) =
xt::col(BS_LBD, 1) =
xt::sum(
xt::square(
bar_y_j - bar_y
@@ -217,7 +232,7 @@ namespace evalhyd
// calculate sharpness =
// $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
xt::col(bs_lbd, 2) =
xt::col(BS_LBD, 2) =
xt::sum(
xt::square(
y_k -
@@ -227,7 +242,7 @@ namespace evalhyd
) / n;
}
// Compute the Brier skill score (BSS).
// Compute the Brier skill score for each time step.
//
// \require o_k:
// Observed event outcome.
@@ -236,15 +251,15 @@ namespace evalhyd
// Mean event observed outcome.
// shape: (thresholds,)
// \require bs:
// Brier score(s).
// shape: (thresholds,)
// Brier scores for each time step.
// shape: (thresholds, time)
// \assign bss:
// 2D array of Brier skill score for each threshold.
// shape: (thresholds, 1)
// 2D array of Brier skill score for each threshold for each time step.
// shape: (thresholds, time)
void Evaluator::calc_bss()
{
// calculate reference Brier score(s)
// $BS_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o})^2$
// $bs_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o})^2$
xt::xtensor<double, 2> bs_ref =
xt::mean(
xt::square(
@@ -255,8 +270,23 @@ namespace evalhyd
);
// return computed Brier skill score(s)
// $BSS = 1 - \frac{BS}{BS_{ref}}
// $bss = 1 - \frac{bs}{bs_{ref}}
bss = 1 - (bs / bs_ref);
}
// Compute the Brier skill score (BSS).
//
// \require bss:
// 2D array of Brier skill score for each threshold for each time step.
// shape: (thresholds, time)
// \assign BSS:
// 2D array of Brier skill score for each threshold.
// shape: (thresholds, 1)
void Evaluator::calc_BSS()
{
// calculate the mean over the time steps
// $BSS = \frac{1}{n} \sum_{k=1}^{n} bss$
BSS = xt::mean(bss, -1, xt::keep_dims);
}
}
}
include/evalhyd/probabilist/evaluator_elements.cpp
+ 14
− 0
#include <xtensor/xmath.hpp>
#include <xtensor/xsort.hpp>
#include "evaluator.h"
@@ -67,5 +68,18 @@ namespace evalhyd
// determine probability of threshold(s) exceedance
y_k = n_frc / n_mbr;
}
// Compute the forecast quantiles from the ensemble members.
//
// \require q_frc:
// 2D array of streamflow forecasts.
// shape: (members, time)
// \assign q_qnt:
// 2D array of streamflow forecast quantiles.
// shape: (quantiles, time)
void Evaluator::calc_q_qnt()
{
q_qnt = xt::sort(q_frc, 0);
}
}
}
include/evalhyd/probabilist/evaluator_quantiles.cpp 0 → 100644
+ 104
− 0
#include <unordered_map>
#include <xtensor/xmath.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xoperation.hpp>
#include "evaluator.h"
namespace eh = evalhyd;
namespace evalhyd
{
namespace probabilist
{
// Compute the quantile scores for each time step.
//
// \require q_obs:
// 2D array of streamflow observations.
// shape: (1, time)
// \require q_qnt:
// 2D array of streamflow forecasts.
// shape: (quantiles, time)
// \assign qs:
// 2D array of the quantile scores for each time step.
// shape: (quantiles, time)
void Evaluator::calc_qs()
{
// define some dimensions
std::size_t n_qnt = q_qnt.shape(0);
// compute the quantile order $alpha$
xt::xtensor<double, 1> alpha =
xt::arange<double>(1., double(n_qnt + 1))
/ double(n_qnt + 1);
// calculate the difference
xt::xtensor<double, 2> diff = q_qnt - q_obs;
// calculate the quantile scores
qs = xt::where(
diff > 0,
2 * (1 - xt::view(alpha, xt::all(), xt::newaxis())) * diff,
- 2 * xt::view(alpha, xt::all(), xt::newaxis()) * diff
);
}
// Compute the quantile score (QS).
//
// \require qs:
// 2D array of quantile scores for each time step.
// shape: (quantiles, time)
// \assign QS:
// 2D array of quantile scores.
// shape: (quantiles, 1)
void Evaluator::calc_QS()
{
// calculate the mean over the time steps
// $QS = \frac{1}{n} \sum_{k=1}^{n} qs$
QS = xt::mean(qs, -1, xt::keep_dims);
}
// Compute the continuous rank probability score(s) based
// on quantile scores for each time step, and integrating using the
// trapezoidal rule.
//
// /!\ The number of quantiles must be sufficiently large so that the
// cumulative distribution is smooth enough for the numerical
// integration to be accurate.
//
// \require qs:
// 2D array of quantile scores for each time step.
// shape: (quantiles, time)
// \assign crps:
// 2D array of CRPS for each time step.
// shape: (1, time)
void Evaluator::calc_crps()
{
// define some dimensions
std::size_t n_qnt = q_qnt.shape(0);
// integrate with trapezoidal rule
crps = xt::view(
// xt::trapz(y, dx=1/(n+1), axis=0)
xt::trapz(qs, 1./(double(n_qnt) + .1), 0),
xt::newaxis(), xt::all()
);
}
// Compute the continuous rank probability score (CRPS) based
// on quantile scores.
//
// \require crps:
// 2D array of CRPS for each time step.
// shape: (quantiles, time)
// \assign CRPS:
// 2D array containing the CRPS.
// shape: (1, 1)
void Evaluator::calc_CRPS()
{
// calculate the mean over the time steps
// $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$
CRPS = xt::mean(crps, -1, xt::keep_dims);
}
}
}