STK_Law_Normal.cpp

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
 
    This program 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 Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/*
 * Project: stkpp::STatistik::Law
 * Purpose: implementation of the Normal Distribution
 * Author:  Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef IS_RTKPP_LIB
 
#include "../include/STK_Law_Normal.h"
#include <Analysis/include/STK_Funct_raw.h>
 
const double a1 = -39.69683028665376;
const double a2 = 220.9460984245205;
const double a3 = -275.9285104469687;
const double a4 = 138.3577518672690;
const double a5 =-30.66479806614716;
const double a6 = 2.506628277459239;
 
const double b1 = -54.47609879822406;
const double b2 = 161.5858368580409;
const double b3 = -155.6989798598866;
const double b4 = 66.80131188771972;
const double b5 = -13.28068155288572;
 
const double c1 = -0.007784894002430293;
const double c2 = -0.3223964580411365;
const double c3 = -2.400758277161838;
const double c4 = -2.549732539343734;
const double c5 = 4.374664141464968;
const double c6 = 2.938163982698783;
 
const double d1 = 0.007784695709041462;
const double d2 = 0.3224671290700398;
const double d3 = 2.445134137142996;
const double d4 = 3.754408661907416;
 
#endif
 
namespace STK
{
 
namespace Law
{
 
 
#ifndef IS_RTKPP_LIB
 
/*  Generate a pseudo Normal random variate. */
Real Normal::rand() const
{ return (sigma_ == 0.) ? mu_ : generator.randGauss(mu_, sigma_);}
 
/*
 *  Give the value of the pdf at x.
 */
Real Normal::pdf(Real const& x) const
{
  // check parameter
  if ( !Arithmetic<Real>::isFinite(x))
    STKDOMAIN_ERROR_1ARG(Normal::pdf,x,invalid argument);
  // trivial case
  if (sigma_ == 0.) return (x==mu_) ? 1. : 0.;
  // compute pdf
  const Real y = (x - mu_)/sigma_;
  return Const::_1_SQRT2PI_ * std::exp(-0.5 * y * y)  / sigma_;
}
 
/*
 * Give the value of the log-pdf at x.
 */
Real Normal::lpdf(Real const& x) const
{
  // check parameter
  if ( !Arithmetic<Real>::isFinite(x))
    STKDOMAIN_ERROR_1ARG(Normal::lpdf,x,invalid argument);
  // trivial case
  if (sigma_ == 0)
    return (x==mu_) ? 0. : -Arithmetic<Real>::infinity();
  // compute lpdf
  const Real y = (x - mu_)/sigma_;
  return -(Const::_LNSQRT2PI_ + std::log(double(sigma_)) + 0.5 * y * y);
}
 
/*
 * The cumulative distribution function at t.
 */
Real Normal::cdf(Real const& t) const
{
  // check parameter
  if (!Arithmetic<Real>::isFinite(t))
    STKDOMAIN_ERROR_1ARG(Normal::cdf,t,invalid argument);
  // trivial case
  if (sigma_ == 0) return (t < mu_) ? 0. : 1.;
  // change of variable
  return (0.5*Funct::erfc_raw(-((t - mu_) / sigma_)*Const::_1_SQRT2_));
}
 
/* The inverse cumulative distribution function at p.*/
Real Normal::icdf(Real const& p) const
{
  // check parameter
  if ( (!Arithmetic<Real>::isFinite(p)) || (p > 1.) || (p < 0.) )
    STKDOMAIN_ERROR_1ARG(Normal::icdf,p,invalid argument);
 // trivial cases
 if (p == 0.)  return -Arithmetic<Real>::infinity();
 if (p == 1.)  return  Arithmetic<Real>::infinity();
 if (p == 0.5) return  mu_;
 
 double p_low =  0.02425, p_high = 1 - p_low;
 double x = 0.0;
 //Rational approximation for lower region.
 if ( p < p_low)
 {
   double q = sqrt(-2*log(p));
   x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
 }
 //Rational approximation for central region.
 if (p_low <= p && p <= p_high)
 {
   double r = (p-0.5)*(p-0.5);
   x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*(p-0.5)/(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1);
 }
 //Rational approximation for upper region.
 if (p_high < p )
 {
   double q = sqrt(-2*log(1-p));
   x = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
 }
 
 /* The relative error of the approximation has absolute value less
    than 1.15e-9.  One iteration of Halley's rational method (third
    order) gives full machine precision... */
// if(( 0 < p)&&(p < 1))
// {
//     double e = cdf(u) - p;
//     double u = e * Const::_SQRT2PI_ * exp(x*x/2);
//     x = x - u/(1 + x*u/2);
// }
 return  mu_ + sigma_ * x;
}
 
/*  Generate a pseudo Normal random variate with the specified parameters.
 *  (static)
 */
Real Normal::rand(Real const& mu, Real const& sigma)
{
#ifdef STK_DEBUG
  // check parameters
  if ( !Arithmetic<Real>::isFinite(mu) || !Arithmetic<Real>::isFinite(sigma) || sigma < 0)
    STKDOMAIN_ERROR_2ARG(Normal::rand,mu,sigma,invalid argument);
#endif
  // return variate
  return (sigma == 0.) ? mu : generator.randGauss(mu, sigma);
}
 
/*
 *  Give the value of the pdf at x.
 */
Real Normal::pdf(Real const& x, Real const& mu, Real const& sigma)
{
#ifdef STK_DEBUG
  // check parameters
  if ( !Arithmetic<Real>::isFinite(mu) || !Arithmetic<Real>::isFinite(sigma) || sigma < 0)
    STKDOMAIN_ERROR_2ARG(Normal::rand,mu,sigma,invalid argument);
#endif
  // check parameter
  if ( !Arithmetic<Real>::isFinite(x))
    STKDOMAIN_ERROR_1ARG(Normal::pdf,x,invalid argument);
  // trivial case
  if (sigma == 0.) return (x==mu) ? 1. : 0.;
  // compute pdf
  const Real y = (x - mu)/sigma;
  return Const::_1_SQRT2PI_ * std::exp(-0.5 * y * y)  / sigma;
}
/*
 * Give the value of the log-pdf at x.
 */
Real Normal::lpdf(Real const& x, Real const& mu, Real const& sigma)
{
#ifdef STK_DEBUG
  // check parameters
  if ( !Arithmetic<Real>::isFinite(mu) || !Arithmetic<Real>::isFinite(sigma) || sigma < 0)
    STKDOMAIN_ERROR_2ARG(Normal::rand,mu,sigma,invalid argument);
#endif
  // check parameter
  if ( !Arithmetic<Real>::isFinite(x))
    STKDOMAIN_ERROR_1ARG(Normal::lpdf,x,invalid argument);
  // trivial case
  if (sigma == 0)
    return (x==mu) ? 0. : -Arithmetic<Real>::infinity();
  // compute lpdf
  const Real y = (x - mu)/sigma;
  return -(Const::_LNSQRT2PI_ + std::log(double(sigma)) + 0.5 * y * y);
}
 
/*
 * The cumulative distribution function at t.
 */
Real Normal::cdf(Real const& t, Real const& mu, Real const& sigma)
{
  // check parameter
  if (!Arithmetic<Real>::isFinite(t))
    STKDOMAIN_ERROR_1ARG(Normal::cdf,t,invalid argument);
  // trivial case
  if (sigma == 0) return (t < mu) ? 0. : 1.;
  // change of variable
  return (0.5*Funct::erfc_raw(-((t - mu) / sigma)*Const::_1_SQRT2_));
}
 
/* The inverse cumulative distribution function at p.*/
Real Normal::icdf(Real const& p, Real const& mu, Real const& sigma)
{
  // check parameter
  if ( (!Arithmetic<Real>::isFinite(p)) || (p > 1.) || (p < 0.) )
    STKDOMAIN_ERROR_1ARG(Normal::icdf,p,invalid argument);
 // trivial cases
 if (p == 0.)  return -Arithmetic<Real>::infinity();
 if (p == 1.)  return  Arithmetic<Real>::infinity();
 if (p == 0.5) return  mu;
 
  double p_low =  0.02425, p_high = 1 - p_low;
  double x = 0.0;
  //Rational approximation for lower region.
  if ( p < p_low)
  {
    double q = sqrt(-2*log(p));
    x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
  }
  //Rational approximation for central region.
  if (p_low <= p && p <= p_high)
  {
    double r = (p-0.5)*(p-0.5);
    x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*(p-0.5) / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1);
  }
  //Rational approximation for upper region.
  if (p_high < p )
  {
    double q = sqrt(-2*log(1-p));
    x = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
  }
 
  /* The relative error of the approximation has absolute value less
     than 1.15e-9.  One iteration of Halley's rational method (third
     order) gives full machine precision... */
  //Pseudo-code algorithm for refinement
  // if(( 0 < p)&&(p < 1))
  // {
  //     double e = cdf(u) - p;
  //     double u = e * Const::_SQRT2PI_ * exp(x*x/2);
  //     x = x - u/(1 + x*u/2);
  // }
  return  mu + sigma * x;
}
 
#endif
} // namespace law
 
} // namespace STK