From faa6a1326e8b7720b4b7a525b2de8ece93073620 Mon Sep 17 00:00:00 2001 From: Seth Linden Date: Wed, 24 Aug 2022 12:07:40 -0600 Subject: [PATCH] Per issue #2206, started stubbing code pieced to calculate crps_emp_fair. SL --- src/libcode/vx_statistics/pair_data_ensemble.cc | 16 +++++++++++++++- src/libcode/vx_statistics/pair_data_ensemble.h | 5 +++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/libcode/vx_statistics/pair_data_ensemble.cc b/src/libcode/vx_statistics/pair_data_ensemble.cc index 3c16d13033..a9df65d59e 100644 --- a/src/libcode/vx_statistics/pair_data_ensemble.cc +++ b/src/libcode/vx_statistics/pair_data_ensemble.cc @@ -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(); @@ -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); @@ -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; @@ -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); @@ -472,6 +476,15 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) { // Store empirical CRPS stats crps_emp_na.add(compute_crps_emp(o_na[i], cur_ens)); + + // Stub in for now? + // Like this? + //crps_emp_fair_na.add(compute_crps_emp(o_na[i], cur_ens) - compute_mean_abs_diff(cur_ens)); + + // Or like this? + // Note mean_abs_diff function should be added to NumArray class + crps_emp_fair_na.add( crps_emp_na[i] - ens_na[i].mean_abs_diff() ); + crpscl_emp_na.add(compute_crps_emp(o_na[i], cur_clm)); // Ensemble mean and standard deviation @@ -811,7 +824,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 // @@ -827,6 +840,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]); diff --git a/src/libcode/vx_statistics/pair_data_ensemble.h b/src/libcode/vx_statistics/pair_data_ensemble.h index 86ad381395..ce4210c0a6 100644 --- a/src/libcode/vx_statistics/pair_data_ensemble.h +++ b/src/libcode/vx_statistics/pair_data_ensemble.h @@ -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]