diff --git a/data/table_files/met_header_columns_V12.0.txt b/data/table_files/met_header_columns_V12.0.txt index 74b1b1f723..b9ad70022c 100644 --- a/data/table_files/met_header_columns_V12.0.txt +++ b/data/table_files/met_header_columns_V12.0.txt @@ -19,7 +19,7 @@ V12.0 : STAT : PJC : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID V12.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]* V12.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]* V12.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]* -V12.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 SPREAD_MD MAE MAE_OERR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS IGN_OERR_CONV IGN_OERR_CORR DS_OERR DS_ADD_OERR DS_MULT_OERR +V12.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 SPREAD_MD MAE MAE_OERR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS IGN_OERR_CONV IGN_OERR_CORR DSS V12.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 V12.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]* V12.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]* diff --git a/docs/Users_Guide/appendixC.rst b/docs/Users_Guide/appendixC.rst index 1b02e201c6..759ab22653 100644 --- a/docs/Users_Guide/appendixC.rst +++ b/docs/Users_Guide/appendixC.rst @@ -1094,12 +1094,12 @@ Called "IGN_CONV_OERR" and "IGN_CORR_OERR" in ECNT output :numref:`table_ES_head TODO: Eric add 2 equations and interpretation advice here -Dawid-Sebastiani Scoring Rules ------------------------------- +Dawid-Sebastiani Scoring Rule +----------------------------- -Called "DS_OERR", "DS_ADD_OERR", and "DS_MULT_OERR" in ECNT output :numref:`table_ES_header_info_es_out_ECNT` +Called "DSS" in ECNT output :numref:`table_ES_header_info_es_out_ECNT` -TODO: Eric add 3 equations and interpretation advice here +TODO: Eric add equation 16 and interpretation advice here IGN --- diff --git a/docs/Users_Guide/ensemble-stat.rst b/docs/Users_Guide/ensemble-stat.rst index 8d39b1fa9e..0036285a2e 100644 --- a/docs/Users_Guide/ensemble-stat.rst +++ b/docs/Users_Guide/ensemble-stat.rst @@ -656,14 +656,8 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described - IGN_CORR_OERR - Error-corrected logarithmic scoring rule (i.e. ignornance score) from Equation 7 of :ref:`Ferro, 2017 ` * - 52 - - DS_OERR + - DSS - Scoring rule from Equation 16 of :ref:`Dawid and Sebastiani, 1999 ` - * - 53 - - DS_ADD_OERR - - Additive observation error scoring rule from Equation 17 of :ref:`Dawid and Sebastiani, 1999 ` - * - 54 - - DS_MULT_OERR - - Multiplicative observation error scoring rule from Equation 18 of :ref:`Dawid and Sebastiani, 1999 ` .. _table_ES_header_info_es_out_RPS: diff --git a/internal/test_unit/hdr/met_12_0.hdr b/internal/test_unit/hdr/met_12_0.hdr index 23d3717cb4..813ecb116d 100644 --- a/internal/test_unit/hdr/met_12_0.hdr +++ b/internal/test_unit/hdr/met_12_0.hdr @@ -19,7 +19,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 CRPS_EMP_FAIR SPREAD_MD MAE MAE_OERR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS IGN_CONV_OERR IGN_CORR_OERR DS_OERR DS_ADD_OERR DS_MULT_OERR +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 SPREAD_MD MAE MAE_OERR BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS IGN_CONV_OERR IGN_CORR_OERR DSS 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_ diff --git a/src/basic/vx_util/stat_column_defs.h b/src/basic/vx_util/stat_column_defs.h index 273c1e5d3d..8e14327ca7 100644 --- a/src/basic/vx_util/stat_column_defs.h +++ b/src/basic/vx_util/stat_column_defs.h @@ -277,7 +277,7 @@ static const char * ecnt_columns [] = { "MAE", "MAE_OERR", "BIAS_RATIO", "N_GE_OBS", "ME_GE_OBS", "N_LT_OBS", "ME_LT_OBS", "IGN_CONV_OERR", "IGN_CORR_OERR", - "DS_OERR", "DS_ADD_OERR", "DS_MULT_OERR" + "DSS" }; static const char * rps_columns [] = { diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index b697b5d49c..8a6bad221f 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -4276,7 +4276,7 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info, // MAE, MAE_OERR, BIAS_RATIO, // N_GE_OBS, ME_GE_OBS, N_LT_OBS, // ME_LT_OBS, IGN_CONV_OERR, IGN_CORR_OERR, - // DS_OERR, DS_ADD_OERR, DS_MULT_OERR + // DSS // at.set_entry(r, c+0, // Total Number of Pairs ecnt_info.n_pair); @@ -4359,14 +4359,8 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info, at.set_entry(r, c+26, // Ignorance Score, observation error corrected ecnt_info.ign_corr_oerr); - at.set_entry(r, c+27, // Dawid-Sebastiani - ecnt_info.ds_oerr); - - at.set_entry(r, c+28, // Dawid-Sebastiani additive error - ecnt_info.ds_add_oerr); - - at.set_entry(r, c+29, // Dawid-Sebastiani multiplicative error - ecnt_info.ds_mult_oerr); + at.set_entry(r, c+27, // Dawid-Sebastiani Score + ecnt_info.dss); return; } diff --git a/src/libcode/vx_statistics/ens_stats.cc b/src/libcode/vx_statistics/ens_stats.cc index 74fc1e1427..2508adda97 100644 --- a/src/libcode/vx_statistics/ens_stats.cc +++ b/src/libcode/vx_statistics/ens_stats.cc @@ -186,7 +186,7 @@ void ECNTInfo::clear() { spread_plus_oerr = bad_data_double; ign_conv_oerr = ign_corr_oerr = bad_data_double; - ds_oerr = ds_add_oerr = ds_mult_oerr = bad_data_double; + dss = bad_data_double; n_ge_obs = n_lt_obs = 0; me_ge_obs = me_lt_obs = bias_ratio = bad_data_double; @@ -227,9 +227,7 @@ void ECNTInfo::assign(const ECNTInfo &c) { ign_conv_oerr = c.ign_conv_oerr; ign_corr_oerr = c.ign_corr_oerr; - ds_oerr = c.ds_oerr; - ds_add_oerr = c.ds_add_oerr; - ds_mult_oerr = c.ds_mult_oerr; + dss = c.dss; n_ge_obs = c.n_ge_obs; n_lt_obs = c.n_lt_obs; @@ -373,9 +371,7 @@ void ECNTInfo::set(const PairDataEnsemble &pd) { // Compute log scores with observational uncertainty ign_conv_oerr = pd.ign_conv_oerr_na.wmean(pd.wgt_na); ign_corr_oerr = pd.ign_corr_oerr_na.wmean(pd.wgt_na); - ds_oerr = pd.ds_oerr_na.wmean(pd.wgt_na); - ds_add_oerr = pd.ds_add_oerr_na.wmean(pd.wgt_na); - ds_mult_oerr = pd.ds_mult_oerr_na.wmean(pd.wgt_na); + dss = pd.dss_na.wmean(pd.wgt_na); // Compute bias ratio terms n_ge_obs = nint(pd.n_ge_obs_na.sum()); diff --git a/src/libcode/vx_statistics/ens_stats.h b/src/libcode/vx_statistics/ens_stats.h index db6e26f346..e7ea1d3146 100644 --- a/src/libcode/vx_statistics/ens_stats.h +++ b/src/libcode/vx_statistics/ens_stats.h @@ -83,9 +83,8 @@ class ECNTInfo { double spread_plus_oerr; // Log scores that incorporate observational uncertainty - // as advised in Ferro (2017) - double ign_conv_oerr, ign_corr_oerr; - double ds_oerr, ds_add_oerr, ds_mult_oerr; + // and Dawid-Sebastiani score, as advised in Ferro (2017) + double ign_conv_oerr, ign_corr_oerr, dss; // Bias ratio information int n_ge_obs, n_lt_obs; diff --git a/src/libcode/vx_statistics/pair_data_ensemble.cc b/src/libcode/vx_statistics/pair_data_ensemble.cc index 6ec1646d2d..6e1603f64c 100644 --- a/src/libcode/vx_statistics/pair_data_ensemble.cc +++ b/src/libcode/vx_statistics/pair_data_ensemble.cc @@ -110,9 +110,7 @@ void PairDataEnsemble::clear() { ign_conv_oerr_na.clear(); ign_corr_oerr_na.clear(); - ds_oerr_na.clear(); - ds_add_oerr_na.clear(); - ds_mult_oerr_na.clear(); + dss_na.clear(); n_ge_obs_na.clear(); me_ge_obs_na.clear(); @@ -187,9 +185,7 @@ void PairDataEnsemble::extend(int n) { pit_na.extend (n); ign_conv_oerr_na.extend (n); ign_corr_oerr_na.extend (n); - ds_oerr_na.extend (n); - ds_add_oerr_na.extend (n); - ds_mult_oerr_na.extend (n); + dss_na.extend (n); n_ge_obs_na.extend (n); me_ge_obs_na.extend (n); n_lt_obs_na.extend (n); @@ -259,9 +255,7 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) { ign_conv_oerr_na = pd.ign_conv_oerr_na; ign_corr_oerr_na = pd.ign_corr_oerr_na; - ds_oerr_na = pd.ds_oerr_na; - ds_add_oerr_na = pd.ds_add_oerr_na; - ds_mult_oerr_na = pd.ds_mult_oerr_na; + dss_na = pd.dss_na; n_ge_obs_na = pd.n_ge_obs_na; me_ge_obs_na = pd.me_ge_obs_na; @@ -470,9 +464,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { pit_na.add(bad_data_double); ign_conv_oerr_na.add(bad_data_double); ign_corr_oerr_na.add(bad_data_double); - ds_oerr_na.add(bad_data_double); - ds_add_oerr_na.add(bad_data_double); - ds_mult_oerr_na.add(bad_data_double); + dss_na.add(bad_data_double); n_ge_obs_na.add(bad_data_double); me_ge_obs_na.add(bad_data_double); n_lt_obs_na.add(bad_data_double); @@ -485,48 +477,42 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { var_unperturbed = compute_variance(esum_na[i], esumsq_na[i], esumn_na[i]); var_na.add(var_unperturbed); + double emn_unperturbed = compute_mean(esum_na[i], esumn_na[i]); + double esd_unperturbed = compute_stdev(esum_na[i], esumsq_na[i], esumn_na[i]); + + // Compute the Dawid Sebastiani scores + dss_na.add( + compute_dawid_sebastiani( + emn_unperturbed, esd_unperturbed, o_na[i])); + // Process the observation error information ObsErrorEntry * e = (has_obs_error() ? obs_error_entry[i] : 0); if(e) { - double emn_unperturbed = compute_mean(esum_na[i], esumn_na[i]); - double esd_unperturbed = compute_stdev(esum_na[i], esumsq_na[i], esumn_na[i]); - double v_conv, v_corr; - // Compute the observation error log scores + double v_conv, v_corr; compute_obs_error_log_scores( - emn_unperturbed, esd_unperturbed, o_na[i], e->variance(), + emn_unperturbed, esd_unperturbed, + o_na[i], e->variance(), v_conv, v_corr); ign_conv_oerr_na.add(v_conv); ign_corr_oerr_na.add(v_corr); - // Compute the Dawid Sebastiani scores - double v_ds, v_ds_add, v_ds_mult; - compute_dawid_sebastiani( - emn_unperturbed, esd_unperturbed, o_na[i], e->variance(), - e->bias_scale, e->bias_offset, - v_ds, v_ds_add, v_ds_mult); - ds_oerr_na.add(v_ds); - ds_add_oerr_na.add(v_ds_add); - ds_mult_oerr_na.add(v_ds_mult); - // Compute perturbed ensemble mean and variance // Exclude the control member from the variance mn_oerr_na.add(cur_ens.mean()); var_oerr_na.add(cur_ens.variance(ctrl_index)); // Compute the variance plus observation error variance. - var_plus_oerr_na.add(var_unperturbed + - dist_var(e->dist_type, - e->dist_parm[0], e->dist_parm[1])); + var_plus_oerr_na.add( + var_unperturbed + + dist_var(e->dist_type, + e->dist_parm[0], e->dist_parm[1])); } // If no observation error specified, store bad data values else { ign_conv_oerr_na.add(bad_data_double); ign_corr_oerr_na.add(bad_data_double); - ds_oerr_na.add(bad_data_double); - ds_add_oerr_na.add(bad_data_double); - ds_mult_oerr_na.add(bad_data_double); mn_oerr_na.add(bad_data_double); var_oerr_na.add(bad_data_double); var_plus_oerr_na.add(bad_data_double); @@ -918,8 +904,7 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o // crps_emp_na, crps_emp_fair_na, spread_md_na, // crpscl_emp_na, crps_gaus_na, crpscl_gaus_na, // ign_na, pit_na, - // ign_conv_oerr, ign_corr_oerr, - // ds_oerr, ds_add_oerr, ds_mult_oerr, + // ign_conv_oerr, ign_corr_oerr, dss, // n_gt_obs_na, me_gt_obs_na, n_lt_obs_na, me_lt_obs_na, // var_na, var_oerr_na, var_plus_oerr_na, // mn_na, mn_oerr_na, e_na @@ -945,9 +930,7 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o pd.pit_na.add(pit_na[i]); pd.ign_conv_oerr_na.add(ign_conv_oerr_na[i]); pd.ign_corr_oerr_na.add(ign_corr_oerr_na[i]); - pd.ds_oerr_na.add(ds_oerr_na[i]); - pd.ds_add_oerr_na.add(ds_add_oerr_na[i]); - pd.ds_mult_oerr_na.add(ds_mult_oerr_na[i]); + pd.dss_na.add(dss_na[i]); pd.n_ge_obs_na.add(n_ge_obs_na[i]); pd.me_ge_obs_na.add(me_ge_obs_na[i]); pd.n_lt_obs_na.add(n_lt_obs_na[i]); @@ -2220,79 +2203,36 @@ void compute_obs_error_log_scores(double emn, double esd, //////////////////////////////////////////////////////////////////////// -void compute_dawid_sebastiani(double emn, double esd, - double obs, double oerr_var, - double oerr_mbias, double oerr_abias, - double &ds_oerr, double &ds_add_oerr, - double &ds_mult_oerr) { +double compute_dawid_sebastiani(double emn, double esd, double obs) { const char *method_name = "compute_dawid_sebastiani() -> "; + double v = bad_data_double; + // Compute the Dawid-Sebastiani scoring rules in // Ferro (2017, Eqs 17 and 18) doi:10.1002/qj.3115 // Dawid and Sebastiani (1999) doi:10.1214/aos/1018031101 // These are the recommended scores in Ferro (2017). - // Equation 16 (no obs uncertainty) - if(is_bad_data(emn) || - is_bad_data(esd) || - is_bad_data(obs) || - is_eq(esd, 0.0)) { - ds_oerr = bad_data_double; - } - else { - ds_oerr = log(esd) + - (obs - emn) * (obs - emn) / - (2.0 * esd * esd); - } + // Check for bad data + if(!is_bad_data(emn) && !is_bad_data(esd) && + !is_bad_data(obs) && !is_eq(esd, 0.0)) { - // Equations 17 and 18 - if(is_bad_data(emn) || - is_bad_data(esd) || - is_bad_data(obs) || - is_bad_data(oerr_var)) { - ds_add_oerr = ds_mult_oerr = bad_data_double; - } - else { - - // Default additive and multiplicative biases to 0 and 1 - double a = (is_bad_data(oerr_abias) ? 0.0 : oerr_abias); - double b = (is_bad_data(oerr_mbias) ? 1.0 : oerr_mbias); - - double b2s2 = 2.0 * b * b * esd * esd; - double ov2 = oerr_var * oerr_var; - - if(is_eq(b2s2, 0.0)) { - ds_add_oerr = ds_mult_oerr = bad_data_double; - } - else { - ds_add_oerr = log(esd) + - ((obs - a - b * emn) * - (obs - a - b * emn) - ov2) / - b2s2; - - ds_mult_oerr = log(esd) + - ((obs - b * emn) * - (obs - b * emn) - obs * obs * ov2 / - (b * b + ov2)) / - b2s2; - } + // Equation 16 + v = log(esd) + + (obs - emn) * (obs - emn) / + (2.0 * esd * esd); } if(mlog.verbosity_level() >= 10) { mlog << Debug(10) << method_name - << "inputs (emn = " << emn + << "for input emn = " << emn << ", esd = " << esd << ", obs = " << obs - << ", oerr_var = " << oerr_var - << ", oerr_mbias = " << oerr_mbias - << ", oerr_abias = " << oerr_abias - << ") and outputs (ds_oerr = " << ds_oerr - << ", ds_add_oerr = " << ds_add_oerr - << ", ds_mult_oerr = " << ds_mult_oerr << ")\n"; + << ", output dss = " << v << "\n"; } - return; + return(v); } //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_statistics/pair_data_ensemble.h b/src/libcode/vx_statistics/pair_data_ensemble.h index b9eaea5837..0b5b2a4c79 100644 --- a/src/libcode/vx_statistics/pair_data_ensemble.h +++ b/src/libcode/vx_statistics/pair_data_ensemble.h @@ -93,9 +93,7 @@ class PairDataEnsemble : public PairBase { NumArray ign_conv_oerr_na; // Error convolved log score [n_obs] NumArray ign_corr_oerr_na; // Error corrected log score [n_obs] - NumArray ds_oerr_na; // Dawid Sebastiani score [n_obs] - NumArray ds_add_oerr_na; // Dawid Sebastiani additive error [n_obs] - NumArray ds_mult_oerr_na; // Dawid Sebastiani multiplicative error [n_obs] + NumArray dss_na; // Dawid Sebastiani score [n_obs] NumArray n_ge_obs_na; // Number of ensemble memebers >= obs [n_obs] NumArray me_ge_obs_na; // Mean error of ensemble members >= obs [n_obs] @@ -334,9 +332,8 @@ extern double compute_bias_ratio(double, double); extern void compute_obs_error_log_scores( double, double, double, double, double &, double &); -extern void compute_dawid_sebastiani( - double, double, double, double, double, double, - double &, double &, double &); +extern double compute_dawid_sebastiani( + double, double, double); //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/stat_analysis/aggr_stat_line.cc b/src/tools/core/stat_analysis/aggr_stat_line.cc index 64b2969258..c687e85a13 100644 --- a/src/tools/core/stat_analysis/aggr_stat_line.cc +++ b/src/tools/core/stat_analysis/aggr_stat_line.cc @@ -2658,9 +2658,7 @@ void aggr_ecnt_lines(LineDataFile &f, STATAnalysisJob &job, m[key].ens_pd.ign_na.add(cur.ign); m[key].ens_pd.ign_conv_oerr_na.add(cur.ign_conv_oerr); m[key].ens_pd.ign_corr_oerr_na.add(cur.ign_corr_oerr); - m[key].ens_pd.ds_oerr_na.add(cur.ds_oerr); - m[key].ens_pd.ds_add_oerr_na.add(cur.ds_add_oerr); - m[key].ens_pd.ds_mult_oerr_na.add(cur.ds_mult_oerr); + m[key].ens_pd.dss_na.add(cur.dss); m[key].ens_pd.n_ge_obs_na.add(cur.n_ge_obs); m[key].ens_pd.me_ge_obs_na.add(cur.me_ge_obs); m[key].ens_pd.n_lt_obs_na.add(cur.n_lt_obs); @@ -3242,24 +3240,6 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, square(cur.spread); } - // Compute observation error log scores - double v_conv, v_corr; - compute_obs_error_log_scores( - cur.ens_mean, cur.spread, cur.obs, oerr_var, - v_conv, v_corr); - m[key].ens_pd.ign_conv_oerr_na.add(v_conv); - m[key].ens_pd.ign_corr_oerr_na.add(v_corr); - - // Compute the Dawid Sebastiani scores - double v_ds, v_ds_add, v_ds_mult; - compute_dawid_sebastiani( - cur.ens_mean, cur.spread, cur.obs, oerr_var, - bad_data_double, bad_data_double, - v_ds, v_ds_add, v_ds_mult); - m[key].ens_pd.ds_oerr_na.add(v_ds); - m[key].ens_pd.ds_add_oerr_na.add(v_ds_add); - m[key].ens_pd.ds_mult_oerr_na.add(v_ds_mult); - // Store BIAS_RATIO terms int n_ge_obs, n_lt_obs; double me_ge_obs, me_lt_obs; @@ -3272,6 +3252,19 @@ void aggr_orank_lines(LineDataFile &f, STATAnalysisJob &job, m[key].ens_pd.n_lt_obs_na.add(n_lt_obs); m[key].ens_pd.me_lt_obs_na.add(me_lt_obs); + // Compute observation error log scores + double v_conv, v_corr; + compute_obs_error_log_scores( + cur.ens_mean, cur.spread, cur.obs, oerr_var, + v_conv, v_corr); + m[key].ens_pd.ign_conv_oerr_na.add(v_conv); + m[key].ens_pd.ign_corr_oerr_na.add(v_corr); + + // Compute the Dawid Sebastiani score + m[key].ens_pd.dss_na.add( + compute_dawid_sebastiani( + cur.ens_mean, cur.spread, cur.obs)); + // // Increment the RHIST counts // diff --git a/src/tools/core/stat_analysis/parse_stat_line.cc b/src/tools/core/stat_analysis/parse_stat_line.cc index dcabc811c9..2eb6d8644d 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/src/tools/core/stat_analysis/parse_stat_line.cc @@ -401,9 +401,7 @@ void parse_ecnt_line(STATLine &l, ECNTData &e_data) { e_data.ign_conv_oerr = atof(l.get_item("IGN_CONV_OERR")); e_data.ign_corr_oerr = atof(l.get_item("IGN_CONV_OERR")); - e_data.ds_oerr = atof(l.get_item("DS_OERR")); - e_data.ds_add_oerr = atof(l.get_item("DS_ADD_OERR")); - e_data.ds_mult_oerr = atof(l.get_item("DS_MULT_OERR")); + e_data.dss = atof(l.get_item("DSS")); return; } diff --git a/src/tools/core/stat_analysis/parse_stat_line.h b/src/tools/core/stat_analysis/parse_stat_line.h index c2ca0d4b77..b8eff792d6 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.h +++ b/src/tools/core/stat_analysis/parse_stat_line.h @@ -72,10 +72,10 @@ struct ECNTData { double me_oerr, mae_oerr, rmse_oerr, spread_oerr; double spread_plus_oerr; double bias_ratio; - double ign_conv_oerr, ign_corr_oerr; - double ds_oerr, ds_add_oerr, ds_mult_oerr; int n_ge_obs, n_lt_obs; double me_ge_obs, me_lt_obs; + double ign_conv_oerr, ign_corr_oerr; + double dss; }; // Ranked Histogram (RHIST) data structure