// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
// ** Copyright UCAR (c) 1992 - 2007
// ** University Corporation for Atmospheric Research (UCAR)
// ** National Center for Atmospheric Research (NCAR)
// ** Research Applications Lab (RAL)
// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*


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

using namespace std;

#include <math.h>

#include "vx_met_util/compute_ci.h"
#include "vx_math/constants.h"

////////////////////////////////////////////////////////////////////////
//
// Compute a confidence interval for a proportion.
//
////////////////////////////////////////////////////////////////////////

void compute_proportion_ci(double p, int n, double alpha, double &p_cl, double &p_cu) {

   //
   // Compute the confidence interval using the Wilson method for all
   // sizes of n, since it provides a better approximation
   //
   compute_wilson_ci(p, n, alpha, p_cl, p_cu);
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Compute a confidence interval for a proportion using the Wald method.
// Taken from Statistical Methods in the Atmospheric Sciences 2nd Ed.
// Daniel S. Wilks, page 327
//
////////////////////////////////////////////////////////////////////////

void compute_wald_ci(double p, int n, double alpha, double &p_cl, double &p_cu) {
   double cv_normal_l, cv_normal_u;

   if(p == bad_data_double) {
      p_cl = p_cu = bad_data_double;
      return;
   }

   //
   // Compute the upper and lower critical values from the
   // normal distribution.
   //
   cv_normal_l = normal_cdf_inv(alpha/2.0, 0.0, 1.0);
   cv_normal_u = normal_cdf_inv(1.0 - (alpha/2.0), 0.0, 1.0);

   //
   // Compute the upper and lower bounds of the confidence interval
   //
   p_cl = p + cv_normal_l * sqrt(p*(1.0-p)/n);
   p_cu = p + cv_normal_u * sqrt(p*(1.0-p)/n);

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Compute a confidence interval for a proportion using the Wilson method.
// Taken from Statistical Methods in the Atmospheric Sciences 2nd Ed.
// Daniel S. Wilks, page 327
//
////////////////////////////////////////////////////////////////////////

void compute_wilson_ci(double p, int n, double alpha, double &p_cl, double &p_cu) {
   double cv_normal_l, cv_normal_u;
   
   if(p == bad_data_double) {
      p_cl = p_cu = bad_data_double;
      return;
   }

   //
   // Compute the upper and lower critical values from the 
   // normal distribution.
   //
   cv_normal_l = normal_cdf_inv(alpha/2.0, 0.0, 1.0);
   cv_normal_u = normal_cdf_inv(1.0 - (alpha/2.0), 0.0, 1.0);
   
   //
   // Compute the upper and lower bounds of the confidence interval
   //
   p_cl = ( p 
            + (cv_normal_l*cv_normal_l)/(2.0*n) 
            + cv_normal_l * sqrt(p*(1.0-p)/n + cv_normal_l*cv_normal_l/(4*n*n)) )
          / (1.0 + cv_normal_l*cv_normal_l/n);
   p_cu = ( p 
            + (cv_normal_u*cv_normal_u)/(2.0*n) 
            + cv_normal_u * sqrt(p*(1.0-p)/n + cv_normal_u*cv_normal_u/(4*n*n)) )
          / (1.0 + cv_normal_u*cv_normal_u/n);

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Compute a confidence interval for the odds ratio using Woolf's 
// method.
//
////////////////////////////////////////////////////////////////////////

void compute_woolf_ci(double odds, double alpha, 
                      int fy_oy, int fy_on, int fn_oy, int fn_on,
                      double &odds_cl, double &odds_cu) {
   double cv_normal_l, cv_normal_u, a, b;
  
   if(odds == bad_data_double ||
      fy_oy == 0 || fy_on == 0 || fn_oy == 0 || fn_on == 0) {
      odds_cl = odds_cu = bad_data_double;
      return;
   }
   
   //
   // Compute the upper and lower critical values from the 
   // normal distribution.
   //
   cv_normal_l = normal_cdf_inv(alpha/2.0, 0.0, 1.0);
   cv_normal_u = normal_cdf_inv(1.0 - (alpha/2.0), 0.0, 1.0);
   
   //
   // Compute the upper and lower bounds of the confidence interval
   //
   a = exp(cv_normal_l*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on));
   b = exp(cv_normal_u*sqrt(1.0/fy_oy + 1.0/fy_on + 1.0/fn_oy + 1.0/fn_on));
   
   odds_cl = odds * a;
   odds_cu = odds * b;
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Compute a confidence interval for the Hanssen and Kuipers
// discriminant.
// Taken from Statistical Methods in the Atmospheric Sciences 2nd Ed.
// Daniel S. Wilks, page 328
//
////////////////////////////////////////////////////////////////////////

void compute_hk_ci(double hk, double alpha, 
                   int fy_oy, int fy_on, int fn_oy, int fn_on,
                   double &hk_cl, double &hk_cu) {
   double cv_normal, stdev;
   double h, f, h_var, f_var;
   int n, h_n, f_n;
  
   //
   // Get the counts
   //
   n   = fy_oy + fy_on + fn_oy + fn_on;
   h_n = fy_oy + fn_oy;
   f_n = fn_on + fy_on;
   
   if(hk == bad_data_double || h_n == 0 || f_n == 0) {
      hk_cl = hk_cu = bad_data_double;
      return;
   }
   
   //
   // Compute the critical value for the normal distribution based 
   // on the sample size
   //
   cv_normal = normal_cdf_inv(alpha/2.0, 0.0, 1.0);
   
   //
   // Compute the hit rate and false alarm rate
   //
   h = (double) fy_oy/h_n;
   f = (double) fy_on/f_n;
         
   //
   // Compute a variance for H and F
   //      
   h_var = sqrt(h*(1.0-h)/h_n + cv_normal*cv_normal/(4.0*h_n*h_n))
         / (1.0 + cv_normal*cv_normal/h_n);

   f_var = sqrt(h*(1.0-h)/f_n + cv_normal*cv_normal/(4.0*f_n*f_n))
         / (1.0 + cv_normal*cv_normal/f_n); 
   
   //
   // Compute the standard deviation for HK
   //
   stdev = sqrt(h_var*h_var + f_var*f_var);
   
   //
   // Compute the upper and lower bounds of the confidence interval
   //   
   hk_cl = hk - cv_normal*stdev;
   hk_cu = hk + cv_normal*stdev;

   return;
}

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