Feature #2058 bias_ratio

data/table_files/met_header_columns_V11.0.txt

@@ -19,7 +19,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 CRPS_EMP_FAIR
+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 BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS
 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]*

docs/Users_Guide/appendixC.rst

@@ -470,6 +470,7 @@ Mean Error (ME)
 ---------------
 
 Called "ME" in CNT output :numref:`table_PS_format_info_CNT`
+Called "ME_OERR", "ME_GE_OBS", and "ME_LT_OBS" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
 
 The Mean Error, ME, is a measure of overall bias for continuous variables; in particular ME = Bias. It is defined as 
 
@@ -1003,6 +1004,19 @@ The continuous ranked probability skill score (CRPSS) is similar to the MSESS an
 
 For the normal distribution fit (CRPSS), the reference CRPS is computed using the climatological mean and standard deviation. For the empirical distribution (CRPSS_EMP), the reference CRPS is computed by sampling from the assumed normal climatological distribution defined by the mean and standard deviation.  
 
+Bias Ratio
+----------
+
+Called "BIAS_RATIO" in ECNT output :numref:`table_ES_header_info_es_out_ECNT`
+
+The bias ratio (BIAS_RATIO) is computed when verifying an ensemble against gridded analyses or point observations. It is defined as the mean error (ME) of ensemble member values greater than or equal to the observation value to which they are matched divided by the absolute value of the mean error (ME) of ensemble member values less than the observation values.
+
+.. math:: \text{BIAS_RATIO} = \frac{ \text{ME}_{f >= o} }{ |\text{ME}_{f < o}| }
+
+A perfect forecast has ME = 0. Since BIAS_RATIO is computed as the high bias (ME_GE_OBS) divide by the absolute value of the low bias (ME_LT_OBS), a perfect forecast has BIAS_RATIO = 0/0, which is undefined. In practice, the high and low bias values are unlikely to be 0.
+
+The range for BIAS_RATIO is 0 to infinity. A score of 1 indicates that the high and low biases are equal. A score greater than 1 indicates that the high bias is larger than the magnitude of the low bias. A score less than 1 indicates the opposite behavior.
+
 IGN
 ---
 

docs/Users_Guide/ensemble-stat.rst

@@ -623,6 +623,21 @@ The format of the STAT and ASCII output of the Ensemble-Stat tool are described
   * - 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)
+  * - 42
+    - BIAS_RATIO
+    - The Bias Ratio
+  * - 43
+    - N_GE_OBS
+    - The number of ensemble values greater than or equal to their observations
+  * - 44
+    - ME_GE_OBS
+    - The Mean Error of the ensemble values greater than or equal to their observations
+  * - 45
+    - N_LT_OBS
+    - The number of ensemble values less than their observations
+  * - 46
+    - ME_LT_OBS
+    - The Mean Error of the ensemble values less than or equal to their observations
 
 .. _table_ES_header_info_es_out_RPS:
 

internal/test_unit/hdr/met_11_0.hdr

@@ -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 CRPS_EMP_FAIR
+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 BIAS_RATIO N_GE_OBS ME_GE_OBS N_LT_OBS ME_LT_OBS
 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_

src/basic/vx_util/stat_column_defs.h

@@ -262,12 +262,14 @@ static const char * isc_columns [] = {
 };
 
 static const char * ecnt_columns [] = {
-   "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"
+   "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", "BIAS_RATIO",
+   "N_GE_OBS",    "ME_GE_OBS",     "N_LT_OBS",
+   "ME_LT_OBS"
 };
 
 static const char * rps_columns [] = {

src/libcode/vx_stat_out/stat_columns.cc

@@ -4220,12 +4220,14 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
    //
    // Ensemble Continuous Statistics
    // Dump out the ECNT line:
-   //    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
+   //    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, BIAS RATIO,
+   //    N_GE_OBS,         ME_GE_OBS,     N_LT_OBS,
+   //    ME_LT_OBS
    //
    at.set_entry(r, c+0,  // Total Number of Pairs
       ecnt_info.n_pair);
@@ -4278,6 +4280,21 @@ void write_ecnt_cols(const ECNTInfo &ecnt_info,
    at.set_entry(r, c+16,  // Empirical ensemble CRPS FAIR
       ecnt_info.crps_emp_fair);
 
+   at.set_entry(r, c+17,  // Bias Ratio
+      ecnt_info.bias_ratio);
+
+   at.set_entry(r, c+18,  // Number of ensemble values >= observations
+      ecnt_info.n_ge_obs);
+
+   at.set_entry(r, c+19,  // ME of ensemble values >= observations
+      ecnt_info.me_ge_obs);
+
+   at.set_entry(r, c+20,  // Number of ensemble values < observations
+      ecnt_info.n_lt_obs);
+
+   at.set_entry(r, c+21,  // ME of ensemble values < observations
+      ecnt_info.me_lt_obs);
+
    return;
 }
 

src/libcode/vx_statistics/compute_stats.cc

@@ -1372,6 +1372,21 @@ void compute_ecnt_mean(const ECNTInfo *ecnt_info, int n,
    for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].spread_plus_oerr);
    ecnt_mean.spread_plus_oerr = na.mean();
 
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].bias_ratio);
+   ecnt_mean.bias_ratio = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].n_ge_obs);
+   ecnt_mean.n_ge_obs = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].me_ge_obs);
+   ecnt_mean.me_ge_obs = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].n_lt_obs);
+   ecnt_mean.n_lt_obs = na.mean();
+
+   for(i=0,na.erase(); i<n; i++) na.add(ecnt_info[i].me_lt_obs);
+   ecnt_mean.me_lt_obs = na.mean();
+
    return;
 }
 

src/libcode/vx_statistics/ens_stats.cc

@@ -183,6 +183,9 @@ void ECNTInfo::clear() {
    me         = rmse       = spread      = bad_data_double;
    me_oerr    = rmse_oerr  = spread_oerr = bad_data_double;
    spread_plus_oerr        = bad_data_double;
+
+   n_ge_obs   = n_lt_obs   = 0;
+   me_ge_obs  = me_lt_obs  = bias_ratio  = bad_data_double;
 
    return;
 }
@@ -215,6 +218,12 @@ void ECNTInfo::assign(const ECNTInfo &c) {
    spread_oerr      = c.spread_oerr;
    spread_plus_oerr = c.spread_plus_oerr;
 
+   n_ge_obs         = c.n_ge_obs;
+   n_lt_obs         = c.n_lt_obs;
+   me_ge_obs        = c.me_ge_obs;
+   me_lt_obs        = c.me_lt_obs;
+   bias_ratio       = c.bias_ratio;
+
    return;
 }
 
@@ -230,7 +239,7 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
    n_ens = pd.n_ens;
 
    // Compute empirical CRPS scores
-   crps_emp   = pd.crps_emp_na.wmean(pd.wgt_na);
+   crps_emp  = pd.crps_emp_na.wmean(pd.wgt_na);
    if(is_eq(crps_emp, 0.0)) crps_emp = 0.0;
 
    crps_emp_fair = pd.crps_emp_fair_na.wmean(pd.wgt_na);
@@ -337,6 +346,15 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
    // Compute the square root of the average variance plus oerr
    spread_plus_oerr = square_root(pd.var_plus_oerr_na.wmean(pd.wgt_na));
 
+   // Compute bias ratio terms 
+   n_ge_obs  = nint(pd.n_ge_obs_na.sum());
+   me_ge_obs = pd.me_ge_obs_na.wmean(pd.n_ge_obs_na);
+   n_lt_obs  = nint(pd.n_lt_obs_na.sum());
+   me_lt_obs = pd.me_lt_obs_na.wmean(pd.n_lt_obs_na);
+
+   // Compute bias ratio
+   bias_ratio = compute_bias_ratio(me_ge_obs, me_lt_obs);
+
    return;
 }
 

src/libcode/vx_statistics/ens_stats.h

@@ -82,6 +82,11 @@ class ECNTInfo {
       double me_oerr, rmse_oerr, spread_oerr;
       double spread_plus_oerr;
 
+      // Bias ratio information
+      int n_ge_obs, n_lt_obs;
+      double me_ge_obs, me_lt_obs;
+      double bias_ratio;
+
       // Compute statistics
       void set(const PairDataEnsemble &);
 

src/libcode/vx_statistics/pair_data_ensemble.cc

@@ -107,6 +107,11 @@ void PairDataEnsemble::clear() {
    ign_na.clear();
    pit_na.clear();
 
+   n_ge_obs_na.clear();
+   me_ge_obs_na.clear();
+   n_lt_obs_na.clear();
+   me_lt_obs_na.clear();
+
    n_ens = 0;
    n_pair = 0;
    ctrl_index = bad_data_int;
@@ -141,6 +146,8 @@ void PairDataEnsemble::clear() {
    me_oerr        = bad_data_double;
    rmse_oerr      = bad_data_double;
 
+   bias_ratio     = bad_data_double;
+
    return;
 }
 
@@ -168,6 +175,10 @@ void PairDataEnsemble::extend(int n) {
    crpscl_gaus_na.extend     (n);
    ign_na.extend             (n);
    pit_na.extend             (n);
+   n_ge_obs_na.extend        (n);
+   me_ge_obs_na.extend       (n);
+   n_lt_obs_na.extend        (n);
+   me_lt_obs_na.extend       (n);
    skip_ba.extend            (n);
    var_na.extend             (n);
    var_oerr_na.extend        (n);
@@ -227,6 +238,12 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) {
    crpscl_gaus_na = pd.crpscl_gaus_na;
    ign_na         = pd.ign_na;
    pit_na         = pd.pit_na;
+
+   n_ge_obs_na    = pd.n_ge_obs_na;
+   me_ge_obs_na   = pd.me_ge_obs_na;
+   n_lt_obs_na    = pd.n_lt_obs_na;
+   me_lt_obs_na   = pd.me_lt_obs_na;
+
    n_pair         = pd.n_pair;
    ctrl_index     = pd.ctrl_index;
    skip_const     = pd.skip_const;
@@ -265,6 +282,8 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) {
    me_oerr        = pd.me_oerr;
    rmse_oerr      = pd.rmse_oerr;
 
+   bias_ratio     = pd.bias_ratio;
+
    set_ens_size(pd.n_ens);
 
    for(i=0; i<n_ens; i++) e_na[i] = pd.e_na[i];
@@ -370,8 +389,8 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
       // Initialize
       cur_ens.erase();
 
-      // Compute the number of ensemble values less than the observation
-      for(j=0, n_vld=0, n_bel=0, n_tie=0; j<n_ens; j++) {
+      // Compute the number of ensemble values above and below the observation
+      for(j=0, n_vld = n_bel = n_tie = 0; j<n_ens; j++) {
 
          // Skip bad data
          if(e_na[j].n() > i && !is_bad_data(e_na[j][i])) {
@@ -382,7 +401,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
             // Store the current ensemble value
             cur_ens.add(e_na[j][i]);
 
-            // Keep track of the number of ties and the number below
+            // Track running counts and sums
             if(is_eq(e_na[j][i], o_na[i])) n_tie++;
             else if(e_na[j][i] < o_na[i])  n_bel++;
          }
@@ -422,6 +441,10 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
          crpscl_gaus_na.add(bad_data_double);
          ign_na.add(bad_data_double);
          pit_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);
+         me_lt_obs_na.add(bad_data_double);
       }
       // Otherwise, compute scores
       else {
@@ -492,6 +515,19 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
          crpscl_gaus_na.add(compute_crps_gaus(o_na[i], cmn_na[i], csd_na[i]));
          ign_na.add(compute_ens_ign(o_na[i], mean, stdev));
          pit_na.add(compute_ens_pit(o_na[i], mean, stdev));
+
+         // Compute the Bias Ratio terms 
+         int n_ge_obs, n_lt_obs;
+         double me_ge_obs, me_lt_obs;
+         compute_bias_ratio_terms(o_na[i], cur_ens,
+                                  n_ge_obs, me_ge_obs,
+                                  n_lt_obs, me_lt_obs);
+
+         // Store the Bias Ratio terms 
+         n_ge_obs_na.add(n_ge_obs);
+         me_ge_obs_na.add(me_ge_obs);
+         n_lt_obs_na.add(n_lt_obs);
+         me_lt_obs_na.add(me_lt_obs);
       }
    } // end for i
 
@@ -820,7 +856,8 @@ 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, 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,
+      //   ign_na, pit_na, 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
       //
       // Exclude from subset:
@@ -841,6 +878,10 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o
       pd.crpscl_gaus_na.add(crpscl_gaus_na[i]);
       pd.ign_na.add(ign_na[i]);
       pd.pit_na.add(pit_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]);
+      pd.me_lt_obs_na.add(me_lt_obs_na[i]);
       pd.skip_ba.add(false);
       pd.var_na.add(var_na[i]);
       pd.var_oerr_na.add(var_oerr_na[i]);
@@ -2006,5 +2047,59 @@ double compute_ens_pit(double obs, double m, double s) {
 
    return(v);
 }
+
+////////////////////////////////////////////////////////////////////////
+//
+// Compute the bias ratio terms
+//
+////////////////////////////////////////////////////////////////////////
+
+void compute_bias_ratio_terms(double obs, const NumArray &e_na,
+                              int &n_ge_obs, double &me_ge_obs,
+                              int &n_lt_obs, double &me_lt_obs) {
+
+   // Initialize
+   n_ge_obs  = n_lt_obs  = 0;
+   me_ge_obs = me_lt_obs = 0.0;
+
+   // Loop over ensemble member values
+   for(int i=0; i<e_na.n(); i++) {
+
+      // Track counts and sums
+      if(e_na[i] >= obs) {
+         n_ge_obs  += 1;
+         me_ge_obs += (e_na[i] - obs);
+      }
+      else {
+         n_lt_obs  += 1;
+         me_lt_obs += (e_na[i] - obs);
+      }
+   }
+
+   // Convert sums to means
+   if(n_ge_obs > 0) me_ge_obs /= n_ge_obs;
+   else             me_ge_obs  = bad_data_double;
+   if(n_lt_obs > 0) me_lt_obs /= n_lt_obs;
+   else             me_lt_obs  = bad_data_double;
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
+double compute_bias_ratio(double me_ge_obs, double me_lt_obs) {
+   double v;
+
+   // Compute bias ratio
+   if(is_bad_data(me_ge_obs) ||
+      is_bad_data(me_lt_obs) || is_eq(me_lt_obs, 0.0)) {
+      v = bad_data_double;
+   }
+   else {
+      v = me_ge_obs / abs(me_lt_obs);
+   }
+
+   return(v);
+}
 
 ////////////////////////////////////////////////////////////////////////