Skip to content
GitLab
Select Git revision
  • dev
  • 11-check-the-computation-of-the-ranks-for-ties
  • main default
  • add-determinist-to-evalp
  • v0.1.2
  • v0.1.1
  • v0.1.0
  • v0.0.1
Find file BlameHistoryPermalink
evaluator_brier.cpp 9.63 KiB
1 
#include <unordered_map>
2 
#include <xtensor/xmath.hpp>
3 
#include <xtensor/xview.hpp>
4 
#include <xtensor/xoperation.hpp>
5
6 
#include "evaluator.h"
7
8 
namespace eh = evalhyd;
9
10 
// NOTE ------------------------------------------------------------------------
11 
// All equations in metrics below are following notations from
12 
// "Wilks, D. S. (2011). Statistical methods in the atmospheric sciences.
13 
// Amsterdam; Boston: Elsevier Academic Press. ISBN: 9780123850225".
14 
// In particular, pp. 302-303, 332-333.
15 
// -----------------------------------------------------------------------------
16
17 
namespace evalhyd
18 
{
19 
    namespace probabilist
20 
    {
21 
        // Compute the Brier score (BS).
22 
        //
23 
        // \require o_k:
24 
        //     Observed event outcome.
25 
        //     shape: (thresholds, time)
26 
        // \require y_k:
27 
        //     Event probability forecast.
28 
        //     shape: (thresholds, time)
29 
        // \assign bs:
30 
        //     2D array of Brier Score for each threshold.
31 
        //     shape: (thresholds, 1)
32 
        void Evaluator::calc_bs()
33 
        {
34 
            // return computed Brier score(s)
35 
            // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
36 
            bs = xt::mean(xt::square(o_k - y_k), -1, xt::keep_dims);
37 
        }
38
39 
        // Compute the calibration-refinement decomposition of the Brier score
40 
        // into reliability, resolution, and uncertainty.
41 
        //
42 
        // BS = reliability - resolution + uncertainty
43 
        //
44 
        // \require q_thr:
45 
        //     1D array of streamflow exceedance threshold(s).
46 
        //     shape: (thresholds,)
47 
        // \require o_k:
48 
        //     Observed event outcome.
49 
        //     shape: (thresholds, time)
50 
        // \require y_k:
51 
        //     Event probability forecast.
52 
        //     shape: (thresholds, time)
53 
        // \assign bs_crd:
54 
        //     2D array of Brier score components (reliability, resolution,
55 
        //     uncertainty) for each threshold.
56 
        //     shape: (thresholds, components)
57 
        void Evaluator::calc_bs_crd()
58 
        {
59 
            // declare internal variables
60 
            // shape: (bins, thresholds, time)
61 
            xt::xtensor<double, 3> msk_bins;
62 
            // shape: (thresholds,)
63 
            xt::xtensor<double, 1> bar_o;
64 
            // shape: (bins, thresholds)
65 
            xt::xtensor<double, 2> N_i, bar_o_i;
66 
            // shape: (bins,)
67 
            xt::xtensor<double, 1> y_i;
68
69 
            // define some dimensions
70 
            std::size_t n = o_k.shape(1);
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
std::size_t n_thr = o_k.shape(0); std::size_t n_mbr = q_frc.shape(0); std::size_t n_cmp = 3; // initialise output variable // shape: (components, thresholds) 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; // compute mask to subsample time steps belonging to same forecast bin // (where bins are defined as the range of forecast values) msk_bins = xt::equal( y_k, xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis()) ); // compute number of forecasts in each forecast bin $N_i$ N_i = xt::sum(msk_bins, -1); // compute subsample relative frequency // $\bar{o_i} = \frac{1}{N_i} \sum_{k \in N_i} o_k$ bar_o_i = xt::where( N_i > 0, xt::sum( xt::where( msk_bins, xt::view(o_k, xt::newaxis(), xt::all(), xt::all()), 0. ), -1 ) / N_i, 0. ); // compute mean "climatology" relative frequency of the event // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$ bar_o = xt::mean(o_k, -1); // calculate reliability = // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$ xt::col(bs_crd, 0) = xt::sum( xt::square( xt::view(y_i, xt::all(), xt::newaxis()) - bar_o_i ) * N_i, 0 ) / n; // calculate resolution = // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$ xt::col(bs_crd, 1) = xt::sum( xt::square( bar_o_i - bar_o ) * N_i, 0 ) / n; // calculate uncertainty = $\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 // into type 2 bias, discrimination, and sharpness (a.k.a. refinement). // // BS = type 2 bias - discrimination + sharpness //
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
// \require o_k: // Observed event outcome. // shape: (thresholds, time) // \require y_k: // Event probability forecast. // shape: (thresholds, time) // \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() { // declare internal variables // shape: (bins, thresholds, time) xt::xtensor<double, 3> msk_bins; // shape: (thresholds,) xt::xtensor<double, 1> bar_y; // shape: (bins, thresholds) xt::xtensor<double, 2> M_j, bar_y_j; // shape: (bins,) xt::xtensor<double, 1> o_j; // define some dimensions std::size_t n = o_k.shape(1); std::size_t n_thr = o_k.shape(0); std::size_t n_cmp = 3; // declare and initialise output variable // shape: (components, thresholds) bs_lbd = xt::zeros<double>({n_thr, n_cmp}); // set the range of observed values $o_j$ o_j = {0., 1.}; // compute mask to subsample time steps belonging to same observation bin // (where bins are defined as the range of forecast values) msk_bins = xt::equal( o_k, xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis()) ); // compute number of observations in each observation bin $M_j$ M_j = xt::sum(msk_bins, -1); // compute subsample relative frequency // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$ bar_y_j = xt::where( M_j > 0, xt::sum( xt::where( msk_bins, xt::view(y_k, xt::newaxis(), xt::all(), xt::all()), 0. ), -1 ) / M_j, 0. ); // compute mean "climatology" forecast probability // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$ bar_y = xt::mean(y_k, -1); // 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::sum( xt::square( xt::view(o_j, xt::all(), xt::newaxis()) - bar_y_j
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
) * M_j, 0 ) / n; // calculate discrimination = // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$ xt::col(bs_lbd, 1) = xt::sum( xt::square( bar_y_j - bar_y ) * M_j, 0 ) / n; // calculate sharpness = // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$ xt::col(bs_lbd, 2) = xt::sum( xt::square( y_k - xt::view(bar_y, xt::all(), xt::newaxis()) ), -1 ) / n; } // Compute the Brier skill score (BSS). // // \require o_k: // Observed event outcome. // shape: (thresholds, time) // \require bs: // Brier score(s). // shape: (thresholds,) // \assign bss: // 2D array of Brier skill score for each threshold. // shape: (thresholds, 1) void Evaluator::calc_bss() { // calculate mean observed event outcome // $\bar{o_k} = \frac{1}{n} \sum_{k=1}^{n} o_k$ xt::xtensor<double, 2> bar_o_k = xt::mean(o_k, -1, xt::keep_dims); // calculate reference Brier score(s) // $BS_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o_k})^2$ xt::xtensor<double, 2> bs_ref = xt::mean( xt::square(o_k - bar_o_k), -1, xt::keep_dims ); // return computed Brier skill score(s) // $BSS = 1 - \frac{BS}{BS_{ref}} bss = 1 - (bs / bs_ref); } } }