Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature 2206 fair crps to ecnt #2247

Merged
merged 25 commits into from
Sep 13, 2022
Merged
Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
b01447c
Per issue #2206, for the ECNT line type, added CRPS_EMP_FAIR. SL
Aug 22, 2022
7b473f4
Per issue #2206, for ECNT line type, added CRPS_EMP_FAIR. SL
Aug 22, 2022
3a18c1d
Per issue #2206, for ecnt_columns, added CRPS_EMP_FAIR. SL
Aug 22, 2022
017448c
Per issue #2206, in write_ecnt_cols(), added crps_emp_fair. SL
Aug 22, 2022
3552644
Per issue #2206, in ECNTInfo, added crps_emp_fair. SL
Aug 22, 2022
faa6a13
Per issue #2206, started stubbing code pieced to calculate crps_emp_f…
Aug 24, 2022
53c72f5
Per issue #2206, added new function mean_abs_diff(). SL
Aug 24, 2022
eb77572
Per issue #2206, in compute_pair_vals() updated the code to calculate…
Aug 25, 2022
7ffe9b4
Per issue #2206, added crps_emp_fair (CRPS_EMP_FAIR) to relevant code…
Aug 29, 2022
1eeed47
Per issue #2206, added wording for CRPS_EMP_FAIR, also added it to la…
Aug 29, 2022
f4039c5
Per issue #2206, added info for CRPS_EMP_FAIR and also added the math…
Aug 29, 2022
6648bcd
Per issue #2206, update the math equation for CRPS_EMP_FAIR. SL
Aug 30, 2022
f16f217
Merge branch 'develop' into feature_2206_fair_crps_to_ecnt
Aug 30, 2022
2c13f4a
Per issue #2206, added new function weighted_mean_absolute_diff(). SL
Sep 7, 2022
51d5f67
Per issue #2206, updated the crps_emp_fair calculation to subtract th…
Sep 7, 2022
44309c6
Per issue #2206, in aggr_orank_lines(), updated calculation for crps_…
Sep 7, 2022
e326222
Merge branch 'develop' into feature_2206_fair_crps_to_ecnt
Sep 7, 2022
3442897
Per issue #2206, renamed weighted_mean_abs_diff() to wmean_abs_diff()…
Sep 8, 2022
48648cf
Per issue #2206, for the crps_emp_fair calculation, changed weighted_…
Sep 8, 2022
c93090d
Per issue #2206, in aggr_orank_lines(), for the crps_emp_fair calcula…
Sep 8, 2022
a451368
Per issue #2206, updated the math equation for the CRPS_EMP_FAIR calc…
Sep 8, 2022
ea93131
Per issue #2206, fixed a bug in the wmean_abs_diff() function. SL
Sep 9, 2022
32c060a
Per issue #2206, added code for crps_emp_fair. SL
Sep 12, 2022
8db8ce8
Merge branch 'develop' into feature_2206_fair_crps_to_ecnt
Sep 12, 2022
f9e17db
Per issue #2206, in set() function, added checks for crps_emp and crp…
Sep 12, 2022
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion data/table_files/met_header_columns_V11.0.txt
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ V11.0 : STAT : PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID
V11.0 : STAT : PRC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) THRESH_[0-9]* PODY_[0-9]* POFD_[0-9]*
V11.0 : STAT : PSTD : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_THRESH) BASER BASER_NCL BASER_NCU RELIABILITY RESOLUTION UNCERTAINTY ROC_AUC BRIER BRIER_NCL BRIER_NCU BRIERCL BRIERCL_NCL BRIERCL_NCU BSS BSS_SMPL THRESH_[0-9]*
V11.0 : STAT : ECLV : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BASER VALUE_BASER (N_PTS) CL_[0-9]* VALUE_[0-9]*
V11.0 : STAT : ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP
V11.0 : STAT : ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR
V11.0 : STAT : RPS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_PROB RPS_REL RPS_RES RPS_UNC RPS RPSS RPSS_SMPL RPS_COMP
V11.0 : STAT : RHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL (N_RANK) RANK_[0-9]*
V11.0 : STAT : PHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BIN_SIZE (N_BIN) BIN_[0-9]*
Expand Down
8 changes: 6 additions & 2 deletions docs/Users_Guide/appendixC.rst
Original file line number Diff line number Diff line change
Expand Up @@ -972,9 +972,9 @@ where :math:`BS_m` is the Brier score for the m-th category (:ref:`Tödter and A
CRPS
----

Called "CRPS", "CRPSCL", "CRPS_EMP", and "CRPSCL_EMP" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
Called "CRPS", "CRPSCL", "CRPS_EMP", "CRPS_EMP_FAIR" and "CRPSCL_EMP" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`

The continuous ranked probability score (CRPS) is the integral, over all possible thresholds, of the Brier scores (:ref:`Gneiting et al., 2004 <Gneiting-2004>`). In MET, the CRPS is calculated two ways: using a normal distribution fit to the ensemble forecasts (CRPS and CRPSCL), and using the empirical ensemble distribution (CRPS_EMP and CRPSCL_EMP). In some cases, use of other distributions would be better.
The continuous ranked probability score (CRPS) is the integral, over all possible thresholds, of the Brier scores (:ref:`Gneiting et al., 2004 <Gneiting-2004>`). In MET, the CRPS is calculated two ways: using a normal distribution fit to the ensemble forecasts (CRPS and CRPSCL), and using the empirical ensemble distribution (CRPS_EMP and CRPSCL_EMP). The empirical ensemble CRPS is also adjusted by subtracting 1/2(m) times the mean absolute difference of the ensemble members (where m is the ensemble size), this is saved as CRPS_EMP_FAIR. In some cases, use of other distributions would be better.

WARNING: The normal distribution is probably a good fit for temperature and pressure, and possibly a not horrible fit for winds. However, the normal approximation will not work on most precipitation forecasts and may fail for many other atmospheric variables.

Expand All @@ -988,6 +988,10 @@ The overall CRPS is calculated as the average of the individual measures. In equ

The score can be interpreted as a continuous version of the mean absolute error (MAE). Thus, the score is negatively oriented, so smaller is better. Further, similar to MAE, bias will inflate the CRPS. Thus, bias should also be calculated and considered when judging forecast quality using CRPS.

To calculate CRPS_EMP_FAIR (bias adjusted, empirical ensemble CRPS)

.. math:: \text{CRPS_EMP_FAIR} = \text{CRPS_EMP} - \frac{1}{n*(n-1)} \sum|f_{i} - f_{j}|

CRPS Skill Score
----------------

Expand Down
5 changes: 4 additions & 1 deletion docs/Users_Guide/ensemble-stat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ Often, the goal of ensemble forecasting is to reproduce the distribution of obse

The relative position (RELP) is a count of the number of times each ensemble member is closest to the observation. For stochastic or randomly derived ensembles, this statistic is meaningless. For specified ensemble members, however, it can assist users in determining if any ensemble member is performing consistently better or worse than the others.

The ranked probability score (RPS) is included in the Ranked Probability Score (RPS) line type. It is the mean of the Brier scores computed from ensemble probabilities derived for each probability category threshold (prob_cat_thresh) specified in the configuration file. The continuous ranked probability score (CRPS) is the average the distance between the forecast (ensemble) cumulative distribution function and the observation cumulative distribution function. It is an analog of the Brier score, but for continuous forecast and observation fields. The CRPS statistic is computed using two methods: assuming a normal distribution defined by the ensemble mean and spread (:ref:`Gneiting et al., 2004 <Gneiting-2004>`) and using the empirical ensemble distribution (:ref:`Hersbach, 2000 <Hersbach-2000>`). The CRPS statistic is included in the Ensemble Continuous Statistics (ECNT) line type, along with other statistics quantifying the ensemble spread and ensemble mean skill.
The ranked probability score (RPS) is included in the Ranked Probability Score (RPS) line type. It is the mean of the Brier scores computed from ensemble probabilities derived for each probability category threshold (prob_cat_thresh) specified in the configuration file. The continuous ranked probability score (CRPS) is the average the distance between the forecast (ensemble) cumulative distribution function and the observation cumulative distribution function. It is an analog of the Brier score, but for continuous forecast and observation fields. The CRPS statistic is computed using two methods: assuming a normal distribution defined by the ensemble mean and spread (:ref:`Gneiting et al., 2004 <Gneiting-2004>`) and using the empirical ensemble distribution (:ref:`Hersbach, 2000 <Hersbach-2000>`). The CRPS statistic using the empirical ensemble distribution can be adjusted (bias corrected) by subtracting 1/2(m) times the mean absolute difference of the ensemble members (where m is the ensemble size), this is saved as a separate stastics called CRPS_EMP_FAIR. The CRPS (and CRPS FAIR) statistic is included in the Ensemble Continuous Statistics (ECNT) line type, along with other statistics quantifying the ensemble spread and ensemble mean skill.

The Ensemble-Stat tool can derive ensemble relative frequencies and verify them as probability forecasts all in the same run. Note however that these simple ensemble relative frequencies are not actually calibrated probability forecasts. If probabilistic line types are requested (output_flag), this logic is applied to each pair of fields listed in the forecast (fcst) and observation (obs) dictionaries of the configuration file. Each probability category threshold (prob_cat_thresh) listed for the forecast field is applied to the input ensemble members to derive a relative frequency forecast. The probability category threshold (prob_cat_thresh) parsed from the corresponding observation entry is applied to the (gridded or point) observations to determine whether or not the event actually occurred. The paired ensemble relative freqencies and observation events are used to populate an Nx2 probabilistic contingency table. The dimension of that table is determined by the probability PCT threshold (prob_pct_thresh) configuration file option parsed from the forecast dictionary. All probabilistic output types requested are derived from the this Nx2 table and written to the ascii output files. Note that the FCST_VAR name header column is automatically reset as "PROB({FCST_VAR}{THRESH})" where {FCST_VAR} is the current field being evaluated and {THRESH} is the threshold that was applied.

Expand Down Expand Up @@ -718,6 +718,9 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described
* - 40
- CRPSS_EMP
- The Continuous Ranked Probability Skill Score (empirical distribution)
* - 41
- CRPS_EMP_FAIR
- The Continuous Ranked Probability Skill Score (empirical distribution) adjusted by subtracting 1/2(m) times the mean absolute difference of the ensemble members (m is the ensemble size)

.. _table_ES_header_info_es_out_RPS:

Expand Down
2 changes: 1 addition & 1 deletion internal/test_unit/hdr/met_11_0.hdr
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L
PRC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_THRESH _VAR_
PSTD : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_THRESH BASER BASER_NCL BASER_NCU RELIABILITY RESOLUTION UNCERTAINTY ROC_AUC BRIER BRIER_NCL BRIER_NCU BRIERCL BRIERCL_NCL BRIERCL_NCU BSS BSS_SMPL _VAR_
ECLV : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BASE N_PTS _VAR_
ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP
ECNT : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_ENS CRPS CRPSS IGN ME RMSE SPREAD ME_OERR RMSE_OERR SPREAD_OERR SPREAD_PLUS_OERR CRPSCL CRPS_EMP CRPSCL_EMP CRPSS_EMP CRPS_EMP_FAIR
RPS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL N_PROB RPS_REL RPS_RES RPS_UNC RPS RPSS RPSS_SMPL RPS_COMP
RHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL CRPS IGN N_RANK CRPSS SPREAD _VAR_
PHIST : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE TOTAL BIN_SIZE N_BIN _VAR_
Expand Down
29 changes: 29 additions & 0 deletions src/basic/vx_util/num_array.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1273,4 +1273,33 @@ ConcatString write_css_hhmmss(const NumArray &na)
}


////////////////////////////////////////////////////////////////////////


double NumArray::mean_abs_diff() const

{

int i, j, count;
double sum, mad;

int n = n_elements();

for(i=0, count=0, sum=0.0; i<n; i++) {
for(j=i+1; j<n; j++) {

if( is_bad_data(e[i]) || is_bad_data(e[j]) ) continue;
sum += abs(e[i]-e[j]);
count++;
}
}

if(count == 0) mad = bad_data_double;
else mad = sum / (n*(n-1));

return(mad);

}


////////////////////////////////////////////////////////////////////////
3 changes: 2 additions & 1 deletion src/basic/vx_util/num_array.h
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,8 @@ class NumArray {
double mean() const;
double mean_sqrt() const;
double mean_fisher() const;

double mean_abs_diff() const;

double variance(int skip_index = bad_data_int) const;
double stdev(int skip_index = bad_data_int) const;

Expand Down
2 changes: 1 addition & 1 deletion src/basic/vx_util/stat_column_defs.h
Original file line number Diff line number Diff line change
Expand Up @@ -266,7 +266,7 @@ static const char * ecnt_columns [] = {
"RMSE", "SPREAD", "ME_OERR",
"RMSE_OERR", "SPREAD_OERR", "SPREAD_PLUS_OERR",
"CRPSCL", "CRPS_EMP", "CRPSCL_EMP",
"CRPSS_EMP"
"CRPSS_EMP", "CRPS_EMP_FAIR"
};

static const char * rps_columns [] = {
Expand Down
5 changes: 4 additions & 1 deletion src/libcode/vx_stat_out/stat_columns.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4044,7 +4044,7 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
// RMSE, SPREAD, ME_OERR,
// RMSE_OERR, SPREAD_OERR, SPREAD_PLUS_OERR,
// CRPSCL, CRPS_EMP, CRPSCL_EMP,
// CRPSS_EMP
// CRPSS_EMP CRPS_EMP_FAIR
//
at.set_entry(r, c+0, // Total Number of Pairs
ecnt_info.n_pair);
Expand Down Expand Up @@ -4094,6 +4094,9 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
at.set_entry(r, c+15, // Empirical CRPSS
ecnt_info.crpss_emp);

at.set_entry(r, c+16, // Empirical ensemble CRPS FAIR
ecnt_info.crps_emp_fair);

return;
}

Expand Down
2 changes: 1 addition & 1 deletion src/libcode/vx_statistics/ens_stats.h
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ class ECNTInfo {
// Number of ensemble members and pairs
int n_ens, n_pair;

double crps_emp, crpscl_emp, crpss_emp;
double crps_emp, crpscl_emp, crpss_emp, crps_emp_fair;
double crps_gaus, crpscl_gaus, crpss_gaus;
double ign, me, rmse, spread;
double me_oerr, rmse_oerr, spread_oerr;
Expand Down
13 changes: 11 additions & 2 deletions src/libcode/vx_statistics/pair_data_ensemble.cc
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,7 @@ void PairDataEnsemble::clear() {
r_na.clear();

crps_emp_na.clear();
crps_emp_fair_na.clear();
crpscl_emp_na.clear();
crps_gaus_na.clear();
crpscl_gaus_na.clear();
Expand Down Expand Up @@ -161,6 +162,7 @@ void PairDataEnsemble::extend(int n) {
v_na.extend (n);
r_na.extend (n);
crps_emp_na.extend (n);
crps_emp_fair_na.extend (n);
crpscl_emp_na.extend (n);
crps_gaus_na.extend (n);
crpscl_gaus_na.extend (n);
Expand Down Expand Up @@ -219,6 +221,7 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) {
v_na = pd.v_na;
r_na = pd.r_na;
crps_emp_na = pd.crps_emp_na;
crps_emp_fair_na = pd.crps_emp_fair_na;
crpscl_emp_na = pd.crpscl_emp_na;
crps_gaus_na = pd.crps_gaus_na;
crpscl_gaus_na = pd.crpscl_gaus_na;
Expand Down Expand Up @@ -413,6 +416,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
var_plus_oerr_na.add(bad_data_double);
r_na.add(bad_data_int);
crps_emp_na.add(bad_data_double);
crps_emp_fair_na.add(bad_data_double);
crpscl_emp_na.add(bad_data_double);
crps_gaus_na.add(bad_data_double);
crpscl_gaus_na.add(bad_data_double);
Expand Down Expand Up @@ -471,7 +475,11 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
derive_climo_vals(cdf_info_ptr, cmn_na[i], csd_na[i], cur_clm);

// Store empirical CRPS stats
crps_emp_na.add(compute_crps_emp(o_na[i], cur_ens));
//
// For crps_emp use temporary, local variable so we can use it for the crps_emp_fair calculation
double crps_emp = compute_crps_emp(o_na[i], cur_ens);
crps_emp_na.add(crps_emp);
crps_emp_fair_na.add(crps_emp - cur_ens.mean_abs_diff());
crpscl_emp_na.add(compute_crps_emp(o_na[i], cur_clm));

// Ensemble mean and standard deviation
Expand Down Expand Up @@ -811,7 +819,7 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o
//
// Include in subset:
// wgt_na, o_na, cmn_na, csd_na, v_na, r_na,
// crps_emp_na, crpscl_emp_na, crps_gaus_na, crpscl_gaus_na,
// crps_emp_na, crps_emp_fair_na, crpscl_emp_na, crps_gaus_na, crpscl_gaus_na,
// ign_na, pit_na, var_na, var_oerr_na, var_plus_oerr_na,
// mn_na, mn_oerr_na, e_na
//
Expand All @@ -827,6 +835,7 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o
pd.v_na.add(v_na[i]);
pd.r_na.add(r_na[i]);
pd.crps_emp_na.add(crps_emp_na[i]);
pd.crps_emp_fair_na.add(crps_emp_fair_na[i]);
pd.crpscl_emp_na.add(crpscl_emp_na[i]);
pd.crps_gaus_na.add(crps_gaus_na[i]);
pd.crpscl_gaus_na.add(crpscl_gaus_na[i]);
Expand Down
5 changes: 3 additions & 2 deletions src/libcode/vx_statistics/pair_data_ensemble.h
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,9 @@ class PairDataEnsemble : public PairBase {
NumArray v_na; // Number of valid ensemble values [n_obs]
NumArray r_na; // Observation ranks [n_obs]

NumArray crps_emp_na; // Empirical Continuous Ranked Probability Score [n_obs]
NumArray crpscl_emp_na; // Empirical climatological CRPS [n_obs]
NumArray crps_emp_na; // Empirical Continuous Ranked Probability Score [n_obs]
NumArray crps_emp_fair_na; // Empirical Continuous Ranked Probability Score [n_obs], adjusted with the mean absolute difference of the ensmble members
NumArray crpscl_emp_na; // Empirical climatological CRPS [n_obs]

NumArray crps_gaus_na; // Gaussian CRPS [n_obs]
NumArray crpscl_gaus_na; // Gaussian climatological CRPS [n_obs]
Expand Down
Loading