Skip to content

Commit

Permalink
Per #2583, remove the DS_ADD_OERR and DS_MULT_OERR ECNT columns and r…
Browse files Browse the repository at this point in the history
…ename DS_OERR as DSS, since observation error is not actually involved in its computation.
  • Loading branch information
JohnHalleyGotway committed Feb 22, 2024
1 parent a217030 commit cc9d981
Show file tree
Hide file tree
Showing 13 changed files with 70 additions and 159 deletions.
2 changes: 1 addition & 1 deletion data/table_files/met_header_columns_V12.0.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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]*
Expand Down
8 changes: 4 additions & 4 deletions docs/Users_Guide/appendixC.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
---
Expand Down
8 changes: 1 addition & 7 deletions docs/Users_Guide/ensemble-stat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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 <Ferro-2017>`
* - 52
- DS_OERR
- DSS
- Scoring rule from Equation 16 of :ref:`Dawid and Sebastiani, 1999 <Dawid-Sebastiani-1999>`
* - 53
- DS_ADD_OERR
- Additive observation error scoring rule from Equation 17 of :ref:`Dawid and Sebastiani, 1999 <Dawid-Sebastiani-1999>`
* - 54
- DS_MULT_OERR
- Multiplicative observation error scoring rule from Equation 18 of :ref:`Dawid and Sebastiani, 1999 <Dawid-Sebastiani-1999>`

.. _table_ES_header_info_es_out_RPS:

Expand Down
2 changes: 1 addition & 1 deletion internal/test_unit/hdr/met_12_0.hdr
Original file line number Diff line number Diff line change
Expand Up @@ -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_
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 @@ -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 [] = {
Expand Down
12 changes: 3 additions & 9 deletions src/libcode/vx_stat_out/stat_columns.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
}
Expand Down
10 changes: 3 additions & 7 deletions src/libcode/vx_statistics/ens_stats.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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());
Expand Down
5 changes: 2 additions & 3 deletions src/libcode/vx_statistics/ens_stats.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
128 changes: 34 additions & 94 deletions src/libcode/vx_statistics/pair_data_ensemble.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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
Expand All @@ -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]);
Expand Down Expand Up @@ -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);
}

////////////////////////////////////////////////////////////////////////
9 changes: 3 additions & 6 deletions src/libcode/vx_statistics/pair_data_ensemble.h
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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);

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

Expand Down
Loading

0 comments on commit cc9d981

Please sign in to comment.