ltqnorm.h

/**************************************************************************
 * Copyright (c) 2004-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 * GNU General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _LTQNORM_H
#define _LTQNORM_H

#ifdef HAVE_GSL
#include "gsl/gsl_cdf.h"
#endif


extern "C" 
{
  double ltqnorm(double p);

  double gaussianPvalue(double r, double mean = 0, double sigma = 1);
}

const double LTQNORM_LARGEST_Z = 8.20954;
// to keep the parallel with Corban's code.
// const double LTQNORM_LARGEST_Z = 30;
const double LTQNORM_SMALLEST_INVERSE = 1.11022e-16;

// largestZ = 8.20954;
// smallestInverse = 1.11022e-16;

inline MyNT computeInverseNormalPvalueTransform(MyNT x)
{
#ifdef HAVE_GSL
  return(gsl_cdf_ugaussian_Qinv(x));
#else
  // \todo HACKHACKHACK. this check should not be based on running
  // ltqnormThreshold() and manually entering the values in ltqhorm.h.
  if (x <= LTQNORM_SMALLEST_INVERSE)
    return(LTQNORM_LARGEST_Z);
  return(ltqnorm(1 - x));
#endif  
}

struct inverseNormalPvalueTransform : public unary_function< MyNT, MyNT > 
{
  MyNT operator()(const MyNT __x) const
    {
      return(computeInverseNormalPvalueTransform(__x));
    }
  
};

inline void ltqnormThreshold(MyNT &x, MyNT &inverse, unsigned int num = 100)
{
  x = 1;
  while (true)
    {
      // pass just x because computeInverseNormalPvalueTransform invokes ltqnorm(1 - x).
      inverse = computeInverseNormalPvalueTransform(x);
//      inverse = ltqnorm(1 - x);
      cout << "x = " << x << ", inverse = " << inverse << endl;
//       if (inverse != inf)
        x = x/2;
//       else
//         break;
        num--;
        if (0 == num)
          break;
    }
  

}



#endif // _LTQNORM_H