calc_rhos.cpp 3.65 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
#include<iostream>
#include<stdio.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_math.h>
#include<gsl/gsl_roots.h>
#include<gsl/gsl_cdf.h>
#include<gsl/gsl_integration.h>
#include "calc_rhos.hpp"

using namespace std;


// function and structure declarations
double integrate_gaussian(double rho0, void *parameters);
struct Qfunction_params { double xth, rho; };
struct int_gauss_params{ double p1, p2, sec_mom; };




////////////////////////////////////////////////////
// calc_rho_given_alpha
// need to calculate value of rho, the correlation
// between two Gaussian random variables
// so that the covariance of the resulting Bernoulli
// variables will be alpha*p1*p2
// solve this 1D equation for rho numerically
// (The variance of each Gaussian is determined by 
// the requirement that the Bernoulli obtained by
// thresholding at 1 has probability pj of being 1.)
/////////////////////////////////////////////////////

double calc_rho_given_alpha(double p1, double p2, double alpha,
			    int &status) {

  // if correlation or a probability is zero, return rho=0
  if(alpha==0 || p1==0 || p2==0) {
    status = 0;
    return 0.0;
  }
  if(alpha < -1) {
    cerr << "alpha < -1, cannot calc rho" << endl;
    status = -1;
    return 0.0;
  }
  
  // desired second moment among bernoulli random variables
  double b_sec_mom = p1*p2*(1+alpha);

  // set up gsl function FF to point to function
  // integrate_gaussian with
  // parameters determined by arguments p and bcorr
  struct int_gauss_params para;
  para.p1 = p1;
  para.p2 = p2;
  para.sec_mom = b_sec_mom;
  gsl_function FF;
  FF.function = &integrate_gaussian;
  FF.params = &para;



  int iter=0, max_iter=1000;
  double rho=0;
  
  // we know rho has to be in interval [x_lo,x_hi]
  double x_lo = -1, x_hi=1;
  if (alpha > 0)
   x_lo=0;
  else x_hi=0;

  // Initialize solver and iterate until converge to solution
  const gsl_root_fsolver_type *T=gsl_root_fsolver_brent;
  gsl_root_fsolver *s = gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(s, &FF, x_lo, x_hi);

  do {
    iter++;
    status = gsl_root_fsolver_iterate (s);
    if(status)
      break;

    rho = gsl_root_fsolver_root (s);
    x_lo = gsl_root_fsolver_x_lower (s);
    x_hi = gsl_root_fsolver_x_upper (s);
    status = gsl_root_test_interval (x_lo, x_hi, 0, 0.00001);
  }
  while (status == GSL_CONTINUE && iter < max_iter);
  
  gsl_root_fsolver_free (s);

  return rho;

}


// Q function, when integrated from xth to infinity,
// gives the second moment of the Bernoulli random variables
double Qfunction(double x, void *parameters){

  struct Qfunction_params *para = (struct Qfunction_params *)parameters;
  double rho = para->rho;
  double xth = para->xth;
  double Q=gsl_cdf_gaussian_Q(xth-rho*x, sqrt(1-rho*rho));
  double g= 1/sqrt(2*M_PI)*exp(-0.5*x*x)*Q;
  return g;

}


// Take the integral of the Gaussian that corresponds to the
// second moment of the Bernoulli random variables.
// Subtract off their required value so that function is 
// zero when found correct value of rho
double integrate_gaussian(double rho0, void *parameters){

  struct int_gauss_params * para = (struct int_gauss_params *) parameters;
  double p1=para->p1;
  double p2=para->p2;
  double sec_mom=para->sec_mom;
  double rho=rho0;
  double xth1=gsl_cdf_ugaussian_Qinv(p1);
  double xth2=gsl_cdf_ugaussian_Qinv(p2);

  struct Qfunction_params qparams;
  qparams.xth=xth1;
  qparams.rho=rho;
  gsl_function G;
  G.function = &Qfunction;
  G.params = &qparams;

  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
  double res, err;
  gsl_integration_qagiu(&G, xth2, 1e-7, 1e-7, 1000, ww, &res, &err);
  gsl_integration_workspace_free (ww);

  res -= sec_mom;
  return res;

}