In this sub-project, we compute usual and special functions: gamma, gammmaRatio, betaRatio,... More...

Collaboration diagram for The Analysis sub-project:

Classes
class	STK::ISerie< Serie >
	Interface base class for Series. More...
Functions
template<class Serie >
home iovleff Developpement workspace stkpp projects STatistiK include STK_Algo h Real	STK::sumAlternateSerie (const ISerie< Serie > &f, Integer const &n=0)
	Sum an alternating serie using the Chebichev polynomials.
template<class Serie >
Real	STK::sumSerie (const ISerie< Serie > &f, Integer const &iter_max=10)
	Sum a serie using the epsilon acceleration process.
template<class Seriea , class Serieb >
Real	STK::continuedFraction (const ISerie< Seriea > &a, const ISerie< Serieb > &b, Integer const &iter_max=100)
	Evaluate a continued fraction.
Real	STK::Const::pi ()
	compute pi
Real	STK::Const::pi_2 ()
	Compute $\pi/2$ using pi()
Real	STK::Const::euler ()
	Compute the Euler constant.
Real	STK::Funct::betaRatio (Real const &a, Real const &b, Real const &x, bool lower_tail=true)
	Compute the incomplete beta function ratio.
Real	STK::Funct::factorial (Integer const &n)
	This function computes for Integer argument.
Real	STK::Funct::factorial (Real const &z)
	This function computes when z is an integer in a Real format.
Real	STK::Funct::factorialLn (Integer const &n)
	This function computes $\ln(n!)$ for Integer argument.
Real	STK::Funct::factorialLn (Real const &z)
	This function computes $\ln(z!)$ when z is an integer in a Real fromat.
Real	STK::Funct::gamma (Real const &z)
	This function computes the function $\Gamma(z)$ .
Real	STK::Funct::gammaLn (Real const &z)
	This function computes the function $\ln(\Gamma(z))$ .
Real	STK::Funct::gammaLnStirlingError (Real const &z)
	Compute the error when we compute $\ln(\Gamma(z))$ using the Stirling's formula.
Real	STK::Funct::gammaLnStirlingError (Integer const &z)
	Compute the error when we compute $\ln(\Gamma(z))$ using the Stirling's formula and z is an Integer.
void	STK::Funct::stirlingCoefficients (STK::Vector &A)
	This function computes the n first coefficients of the Stirling's serie.
Real	STK::Funct::gammaRatio (Real const &a, Real const &x, const bool &lower_tail)
	Compute the incomplete gamma functions ratio.
Real	STK::Funct::gammaRatioQ (Real const &a, Real const &x)
	Compute the incomplete gamma function ratio Q(a,x).
Real	STK::Funct::gammaRatioP (Real const &a, Real const &x)
	Compute the incomplete gamma function ratio P(a,x).
Real	STK::Funct::poisson_pdf_raw (Real const &x, Real const &lambda)
	Compute the poisson density.
Real	STK::Funct::poisson_pdf_raw (Integer const &x, Real const &lambda)
	Compute the poisson density.
Real	STK::Funct::erf_raw (Real const &a)
	Compute the error fonction erf(a)
Real	STK::Funct::erfc_raw (Integer const &a)
	Compute the complementary error function erfc(a)
Real	STK::Funct::normal_cdf_raw (Real const &x)
	Compute the cumulative distribution function of the normal density.
Real	STK::Funct::normal_pdf_raw (Real const &x)
	compute the probability distribution function of the normal density
Real	STK::Funct::g0 (Real const &x)
	compute the partial deviance $g_0(x) = x\log(x)+ 1 - x$ .
Real	STK::Funct::dev0 (Real const &a, Real const &b)
	compute the partial deviance $d_0(a,b) = a\log(a/b)+ b - a$ .
Real	STK::Funct::log1p (Real const &x)
	compute the fonction $\log(1+x)$ .
Real	STK::Funct::expm1 (Real const &x)
	compute the fonction $\exp(x)-1$ .
static Real	STK::Funct::betaRatio_cf (Real const &a, Real const &b, Real const &x, bool lower_tail=true, Integer const &iterMax=1000)
	Compute the continued fraction of the beta function.
static Real	STK::Funct::betaRatio_ae (Real const &a, Real const &b, Real const &x, bool xm1, bool lower_tail)
	Compute the incomplete beta function ratio I_x(a,b) using its asymptotic expansion.
static Real	STK::Funct::serie_up (Real const &s, Real const &a, Real const &x, Integer const &n)
	Compute the recurrence formula of the incomplete beta ratio function.
static Real	STK::Funct::betaRatio_up (Real const &a, Real const &b, Real const &x, bool xm1, bool lower_tail)
	Compute the incomplete beta function ratio I_x(a,b) using its recurrence formula and its asymptotic expansion.
static Real	STK::Funct::betaRatio_sr (Real const &a, Real const &b, Real const &x, bool lower_tail)
	Compute the incomplete beta function ratio I_x(a,b) using its series representation.
static Real	STK::Funct::coefs_odd_se (Real const &std, Real const &qmp, Vector &A)
	compute the odd coefs of the beta Ratio function serie expansion.
static Real	STK::Funct::coefs_even_se (Real const &std, Real const &qmp, Vector &A)
	compute the even coefs of the beta Ratio function serie expansion.
static Real	STK::Funct::betaRatio_se (Real const &a, Real const &b, Real const &x, bool xm1, bool lower_tail, Integer const &iterMax=20)
	Compute the incomplete beta function ratio I_x(a,b) using its serie expansion.
Real	STK::Funct::erfc_raw (Real const &a)
	Compute the function $\mathrm{erfc}(a) = \frac{2}{\sqrt{\pi}} \int_a^{+\infty} e^{-t^2} dt$ .
static Real	STK::Funct::lanczosSerie (Real const &z)
	Compute the Lanzcos correction serie for the gamma function with n = 21 terms.
static Real	STK::Funct::gammaLanczos (Real const &z)
	Compute the gamma function using the Lanzcos expansion using n = 21 terms and r= 22.618910.
static double	STK::Funct::stirlingSerie (Real const &z)
	Compute the Stirling's serie for the gammaLn function.
static Real	STK::Funct::gammaStirling (Real const &z)
	This function computes the gamma function using the Stirling approximation.
static Real	STK::Funct::gammaLnStirling (Real const &z)
	This function computes the log gamma function using the Stirling approximation.
static Real	STK::Funct::apois (Real const &a, Real const &b)
	Compute the poisson density up to a factor.
static Real	STK::Funct::gammaRatio_dl (Real const &a, Real const &x, const bool &lower_tail)
	Compute the incomplete gamma function ratio Q(a,x) using the Taylor serie development representation $Q(a, x) = 1- \frac{x^a}{\Gamma(a+1)} - \frac{x^a}{\Gamma(a)} \sum_{n=1}^\infty (-1)^n\frac{x^n}{(a+n)n!}$ .
static Real	STK::Funct::gammaRatio_cf (Real const &a, Real const &x, const bool lower_tail)
	Compute the incomplete gamma function ratio Q(a, x) using its continued fraction representation.
static Real	STK::Funct::gammaRatio_sr (Real const &a, Real const &x, const bool &lower_tail)
	Compute the incomplete gamma function ratio P(a,x) using the serie development representation.
static Real	STK::Funct::gammaRatio_ae (Real const &a, Real const &x, const bool &lower_tail)
	Compute the incomplete gamma function ratio Q(a,x) using its asymptotic expansion.
static Real	STK::Funct::poisson_ae (Real const &a1, Real const &apd, const bool &lower_tail=true)
	Compute the incomplete gamma function ratio P(a,x) using the Poisson asymptotic expansion.
Variables
static const Real	STK::Funct::factorialArray [51]
	array for the 51th fisrt factorial elements.
static const Real	STK::Funct::factorialLnArray [51]
	array for the 51th fisrt ln factorial elements.
static const Real	STK::Funct::factorialHalvesArray [50]
	array for the 51th fisrt halves factorial elements.
static const Real	STK::Funct::factorialLnHalvesArray [50]
	array for the 51th fisrt halves ln factorial elements.
static const Real	STK::Funct::gammaLnStirlingErrorArray [100]
	array of the gammaLnStirlingError approximation for the values $n=1, 2, 3, \ldots, 99$ .
static const Real	STK::Funct::gammaLnStirlingErrorHalvesArray [100]
	array of the gammaLnStirlingError for the values $z=0.5, 1.5, , \ldots, 99.5$ .
static const Real	STK::Funct::lanczosCoefArray [21]
	array of the Lanzcos coefficients.
static const Real	STK::Funct::stirlingCoefArray [9]
	array of the Stirling coefficients.

Detailed Description

In this sub-project, we compute usual and special functions: gamma, gammmaRatio, betaRatio,...

Function Documentation

template<class Serie >

home iovleff Developpement workspace stkpp projects STatistiK include STK_Algo h Real STK::sumAlternateSerie	(	const ISerie< Serie > &	f,
		Integer const &	n = `0`
	)

Sum an alternating serie using the Chebichev polynomials.

Compute

$S = \sum_{k=0}^{+\infty} (-1)^k a_k$

where the coefficient of the series are given or computed using the parameter ISerie. A Serie must have the following protoype

 class Serie : public ISerie<Serie>
 {
   // return the first coefficient
   inline Integer first();
   // return the first coefficient
   inline Integer next();
 }

The number of of iterations is the first number such that

$\frac{2}{(3+\sqrt(8))^n} < \epsilon,$

where $\epsilon$ is the precision of the type Real.

Parameters:

f	the ISerie giving the terms of the serie
n	the number of iterations

Definition at line 78 of file STK_Algo.h.

References STK::Const::_LN2_, STK::Const::_SQRT2_, and STK::sum().

Referenced by STK::Funct::gammaRatio_dl(), and STK::Funct::log1p().

{
  // Compute the rate of convergence
  Real rate = 3.0 + 2.0 * Const::_SQRT2_;
  // compute the number of iterations if needed
  Integer niter = n;
  if (niter==0)
  {
    niter = (Integer )( log(Arithmetic<Real>::epsilon())
                  / (Const::_LN2_-log(rate))
                  ) + 1;
  }
  // initialization
  rate   = pow(rate, niter);
  Real b = -2.0/(rate+1.0/rate), c = -1.0, sum = 0.0;
  // first iteration (k = 0)
  sum += (c = b-c) * f.first();
  b   *= (- 2.*niter*niter);
  // iterations
  for (Integer k = 1; k<niter; k++)
  {
    sum += (c = b-c) * f.next();
    Real aux = k+1.;
    b   *=   ((k + niter)/(k + aux)) * ((k - niter)/aux) * 2.;
  }
  return sum;
}

template<class Serie >

Real STK::sumSerie	(	const ISerie< Serie > &	f,
		Integer const &	iter_max = `10`
	)

Sum a serie using the epsilon acceleration process.

Compute

$S = \sum_{k=0}^{+\infty} a_k$

where the coefficients of the series are given or computed by the templated parameter Serie.

 class Serie : public ISerie<Serie>
 {
   // return the first coefficient
   inline Integer first();
   // return the first coefficient
   inline Integer next();
 }

The series should be monotone in absolute value. We use the epsilon algorithm acceleration process with 6 epsilon.

Parameters:

f	the functor giving the terms of the serie
iter_max	the number of iterations

Definition at line 131 of file STK_Algo.h.

References STK::sum().

{
  Real e0[7]; // original serie
  Real e1[6];
  Real e2[5];
  Real e3[4];
  Real e4[3];
  Real e5[2];
  Real delta, sum; // e6

  e0[0] = f.first();        //f[0];
  e0[1] = e0[0] + f.next();
  e0[2] = e0[1] + f.next();
  e0[3] = e0[2] + f.next();
  e0[4] = e0[3] + f.next();
  e0[5] = e0[4] + f.next();
  e0[6] = e0[5] + f.next();
  
  // first epsilon e_{-1}[n] = 0
  e1[0] = 1/(e0[1] - e0[0]);
  e1[1] = 1/(e0[2] - e0[1]);
  e1[2] = 1/(e0[3] - e0[2]);
  e1[3] = 1/(e0[4] - e0[3]);
  e1[4] = 1/(e0[5] - e0[4]);
  e1[5] = 1/(e0[6] - e0[5]);
  
  // second epsilon
  e2[0] = e0[0] + 1/(e1[1] - e1[0]);
  e2[1] = e0[1] + 1/(e1[2] - e1[1]);
  e2[2] = e0[2] + 1/(e1[3] - e1[2]);
  e2[3] = e0[3] + 1/(e1[4] - e1[3]);
  e2[4] = e0[4] + 1/(e1[5] - e1[4]);

  // third epsilon
  e3[0] = e1[0] + 1/(e2[1] - e2[0]);
  e3[1] = e1[1] + 1/(e2[2] - e2[1]);
  e3[2] = e1[2] + 1/(e2[3] - e2[2]);
  e3[3] = e1[3] + 1/(e2[4] - e2[3]);

  // fourth epsilon
  e4[0] = e2[0] + 1/(e3[1] - e3[0]);
  e4[1] = e2[1] + 1/(e3[2] - e3[1]);
  e4[2] = e2[2] + 1/(e3[3] - e3[2]);

  // fifth epsilon
  e5[0] = e3[0] + 1/(e4[1] - e4[0]);
  e5[1] = e3[1] + 1/(e4[2] - e4[1]);

  // sum e6[0]
  sum = e4[0] + (delta = 1/(e5[1] - e5[0]));

  // main loop
  for (int n=1; n<=iter_max; n++)
  { // roll s0
    e0[4] = e0[5];
    e0[5] = e0[6];
    e0[6] += f.next();
  
    // roll first epsilon
    e1[3] = e1[4];
    e1[4] = e1[5];
    if ( (delta = (e0[6] - e0[5])) == 0) break;
    e1[5] = 1./delta;
  
    // roll second epsilon
    e2[2] = e2[3];
    e2[3] = e2[4];
    if ( (delta = (e1[5] - e1[4])) == 0) break;
    e2[4] = e0[4] + 1./delta;

    // roll third epsilon
    e3[1] = e3[2];
    e3[2] = e3[3];
    if ( (delta = (e2[4] - e2[3])) == 0) break;
    e3[3] = e1[3] + 1./delta;

    // roll fourth epsilon
    e4[0] = e4[1];
    e4[1] = e4[2];
    if ( (delta = (e3[3] - e3[2])) == 0) break;
    e4[2] = e2[2] + 1./delta;

    // roll fifth epsilon
    e5[0] = e5[1];
    if ( (delta = (e4[2] - e4[1])) == 0) break;
    e5[1] = e3[1] + 1./delta;

    // roll sixth epsilon and compute sum
    if ( (delta = (e5[1] - e5[0])) == 0) break;
    if (!Arithmetic<Real>::isFinite(delta = 1./delta)) break;
    // sum
    Real sum1 = e4[0] + delta;
    if (sum1 == sum) break;
    sum = sum1;
  }
  
  return sum;
}

template<class Seriea , class Serieb >

Real STK::continuedFraction	(	const ISerie< Seriea > &	a,
		const ISerie< Serieb > &	b,
		Integer const &	iter_max = `100`
	)

Evaluate a continued fraction.

Compute

$S = \frac{a_1}{b_1+}\frac{a_2}{b_2+}\frac{a_3}{b_3+}\frac{a_4}{b_4+} \ldots$

where the coefficients of the series are given or computed by the templated parameter Serie.

 class Serie : public ISerie<Serie>
 {
   // return the first coefficient
   inline Integer first();
   // return the first coefficient
   inline Integer next();
 }

Parameters:

iter_max	the number of iterations
a	Denominator serie
b	Numerator serie

Definition at line 253 of file STK_Algo.h.

References STK::abs(), STK::ISerie< Serie >::first(), and STK::ISerie< Serie >::next().

{
  // initialize a_n
  Real an = a.first(); // a_1
  // note a_1 =0 => cf = 0
  if (an == 0.) return 0.;
  // initialize b_n
  Real bn = b.first(); // b_1
  // initialize numerator
  Real Nold = 0, Ncur = an; // = a_1
  // initialize denominator
  Real Dold = 1, Dcur = bn; // = b_1
  // result
  Real cf = 0.;
  // iterations
  for (Integer n=1; n<=iter_max; n++)
  {
    // compute a_n
    an = a.next();
    // check trivial case
    // (a_n =0) and (D_n = 0) => cf = +\infty
    if (an == 0.) return Ncur/Dcur;
    // update b_n
    bn = b.next();
    // compute numerator and denominator
    Real Nnew = bn * Ncur + an * Nold;
    Real Dnew = bn * Dcur + an * Dold;
    // update old numerator and denominator
    Nold = Ncur;
    Dold = Dcur;
    // update current numerator and denominator
    Ncur = Nnew;
    Dcur = Dnew;
    // normalize if necessary with 2^32
    if (abs(Dcur) > 4294967296.0)
    {
      Ncur /= Dcur;
      Nold /= Dcur;
      Dold /= Dcur;
      Dcur = 1.;
    }
    // normalize if necessary with 2^32
    if (abs(Ncur) > 4294967296.0)
    {
      Ncur /= 4294967296.0;
      Nold /= 4294967296.0;
      Dold /= 4294967296.0;
      Dcur /= 4294967296.0;
    }
    // if D_n not to small check cv
    if (abs(Dcur) != 0)
    {
      Real cfnew = Ncur/Dcur;
      // check cv
      if (abs(cf - cfnew) < abs(cfnew)*Arithmetic<Real>::epsilon())
        return cfnew;
      else
        cf = cfnew;
    }
  }
  // avoid warning at compilation
  return cf;
}

Real STK::Const::pi ( )

compute pi

Compute $\pi$ using the Bailey et al.

(Bailey et al. 1997, Adamchik and Wagon 1997) formula.

$\pi = \sum_{k=0}^{+\infty} \left(\frac{1}{16}\right)^k \left( \frac{4}{8k+1} - \frac{2}{8k+4} - \frac{1}{8k+5} - \frac{1}{8k+6} \right)$

Definition at line 118 of file STK_Const_Math.cpp.

Referenced by STK::Const::pi_2().

{
  SeriePi f;
  return sumSerie<SeriePi>(f);
}

Real STK::Const::pi_2 ( ) [inline]

Compute $\pi/2$ using pi()

Definition at line 108 of file STK_Const_Math.h.

References STK::Const::pi().

{ return pi()/2.0;}

Real STK::Const::euler ( )

Compute the Euler constant.

Compute the euler constant $\gamma$ using the formula

$\gamma = \frac{\log(2)}{2} + \frac{1}{\log(2)} \sum_{k=2}^{+\infty} \frac{\log(k)}{k}$

.

Definition at line 167 of file STK_Const_Math.cpp.

References STK::Const::_LN2_.

{
  SerieEuler f;
  return Const::_LN2_/2.0 + sumAlternateSerie<SerieEuler>(f)/Const::_LN2_;
}

Real STK::Funct::betaRatio	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	lower_tail
	)

Compute the incomplete beta function ratio.

Compute the beta ratio function.

$I_x(a,b) = \int_0^x u^{a-1} (1-u)^{b-1} du$

for $0\leq x \leq 1$ .

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 716 of file STK_Funct_betaRatio.cpp.

                                   : 1.;
  if (x>=1) return lower_tail ? 1. : 0.;
  // compute 1-x
  Real y = (0.5 - x) + 0.5;
  // parameters
  Real p=a/(a+b);
  Real q=b/(a+b);
  // small a or b
  if (min(a,b)<=1)
  {
    // case b<=1 and a>=15
    if (a>=15)
      return betaRatio_ae(a, b, x, false, lower_tail);
    // case a<=1 and b>=15
    if (b>=15)
      return betaRatio_ae(b, a, x, true, !lower_tail);
    // ((a<15) and (b<15) and (min(a,b)<=1))
    return (x<=0.5) ? betaRatio_sr(a, b, x, lower_tail)
                    : betaRatio_sr(b, a, y, !lower_tail);
  }
  // (1<a<=15) and (1<b<=15)
  if (max(a,b)<=15)
  {
    // general case
    return (x < p) ? betaRatio_cf(a, b, x,  lower_tail)
                   : betaRatio_cf(b, a, y, !lower_tail);
  }
  // ((1<a<=15) and (15<b)) or ((15<a) and (1<b<=15))
  if (min(a,b)<=15)
  {
    return (a<=15) ? betaRatio_up(a, b, x, false, lower_tail)
                   : betaRatio_up(b, a, x, true, !lower_tail);
  }

//  Real x0, y0, a0, b0, p0, q0;
//  bool flag;
//  if (x>p)
//  { x0=y; y0=x; a0=b; b0=a; flag = !lower_tail; p0=q; q0=p;}
//  else
//  { x0=x; y0=y; a0=a; b0=b; p0=p; q0=q; flag = lower_tail;}
//  // center region for x
//  if ((x0>1.03*p0)&&(y0<0.97*q0))
//    return betaRatio_se(a0-1, b0-1, x0, false, flag);
//  return betaRatio_cf(a0, b0, x0, flag);
  // (a>15) and (b>15) and (a<b)
  if (a<b)
  {
    // center region for x in [p , q]
    if ((x>=p-0.01*(q-p))&&(y<=q+0.01*(q-p)))
      return betaRatio_se(a-1, b-1, x, false, lower_tail);
    // extreme values
    return (x<p) ? betaRatio_cf(a, b, x,  lower_tail)
                 : betaRatio_cf(b, a, y, !lower_tail);
  }
  // (a>15) and (b>15) and (b<a)
  // center region for x in [q, p]
  if ((x>=q-0.01*(p-q))&&(y<=p+0.01*(p-q)))
    return betaRatio_se(b-1, a-1, x, true, !lower_tail);
  // extreme values
  return (x<p) ? betaRatio_cf(a, b, x,  lower_tail)
               : betaRatio_cf(b, a, y, !lower_tail);
}

} // namespace Funct

} // namespace STK

Real STK::Funct::factorial ( Integer const & n )

This function computes $ n! $ for Integer argument.

Compute $n!=1\times 2\times \ldots \times n$ using the values stored in factorialArray for n<51 and using the gamma function for n>50.

Parameters:

n	given value for the factorial function

Definition at line 671 of file STK_Funct_gamma.cpp.

References STK::Funct::factorialArray, and STK::Funct::gamma().

{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(n)) return n;
  // Negative reals argument not allowed
  if (n < 0)
    throw std::domain_error("Funct::factorial(n) "
                       "Negative Integer argument");
  // if n is less than 51 return a stored value
  if (n < 51) return factorialArray[n];
  else        return gamma(n + 1.);
}

Real STK::Funct::factorial ( Real const & z )

This function computes $ z! $ when z is an integer in a Real format.

Compute $z!=1\times 2\times \ldots \times z$ using the values stored in factorialArray for n<51 and using the gamma function for n>50.

Parameters:

z	given value for the factorial function

Definition at line 690 of file STK_Funct_gamma.cpp.

References STK::Funct::factorialArray, and STK::Funct::gamma().

{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  Real n = floor(z);
  // Negative integers or reals arguments not allowed
  if ((n < 0)||(n != z))
    throw std::domain_error("Funct::factorial(x) "
                       "negative or not discrete argument");

  // if n is a small number return a stored value
  if (n < 51) return factorialArray[(Integer)n];
  else        return gamma(z + 1.);
}

Real STK::Funct::factorialLn ( Integer const & n )

This function computes $\ln(n!)$ for Integer argument.

Compute $\ln(n!)=\ln(2)+ \ldots \ln(n)$ using the values stored in factorialLnArray for n<51 and using the gamma function for n>50.

Parameters:

n	given value for the factorial function

Definition at line 711 of file STK_Funct_gamma.cpp.

References STK::Funct::factorialLnArray, and STK::Funct::gammaLn().

{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(n)) return n;
  // Negative reals argument not allowed
  if (n < 0)
    throw std::domain_error("Funct::factorial(n) "
                       "Negative Integer argument");
  // if n is a small number return a constant
  if (n < 51) return factorialLnArray[n];
  else        return gammaLn( n + 1.);
}

Real STK::Funct::factorialLn ( Real const & z )

This function computes $\ln(z!)$ when z is an integer in a Real fromat.

Compute $\ln(z!)=\ln(2)+ \ldots \ln(z)$ using the values stored in factorialLnArray for n<51 and using the gammaLn function for n>50.

Parameters:

z	given value for the factorial function

Definition at line 730 of file STK_Funct_gamma.cpp.

References STK::Funct::factorialLnArray, and STK::Funct::gammaLn().

{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  Real n = floor(z);
  // Negative integers or reals arguments not allowed
  if ((n < 0)||(n != z))
    throw std::domain_error("Funct::factorial(x) "
                       "negative or not discrete argument");
  // if n is a small number return a constant
  if (n < 51) return factorialLnArray[(Integer )n];
  else        return gammaLn(z + 1.);

}

Real STK::Funct::gamma ( Real const & z )

This function computes the function $\Gamma(z)$ .

The gamma function is valid when z is non zero nor a negative integer.

The negative part is computed using the reflection formula

$\Gamma(z) \Gamma(1-z) = \frac{\pi}{\sin(\pi z)}.$

if |z| <8 we use the gamma Lanczos method, else we use the Stirling approxiamtion method.

Parameters:

z	given value for the gamma function

Definition at line 755 of file STK_Funct_gamma.cpp.

References STK::Const::_PI_, STK::abs(), STK::Funct::factorialArray, STK::Funct::factorialHalvesArray, STK::Funct::gammaLanczos(), STK::Funct::gammaStirling(), and STK::isEven().

Referenced by STK::Funct::factorial(), and STK::Funct::gammaLn().

{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if (( z<=0 && z == floor(z))||(Arithmetic<Real>::isInfinite(z)))
    throw std::domain_error("Funct::gamma(z) "
                       "Null or Negative integer argument");
  // If z is an integer (z integer and z<0 is not possible)
  if ((z == floor(z)))
  {
    // if n is a small number return a constant else stirling is accurate
    if (z < 51) return factorialArray[(Integer )z-1];
    else        return gammaStirling(z);
  }
  // compute the absolute value of x -> y
  Real y = abs(z);
  Real value;
  // compute the sign of the gamma function
  Real signgam = (z<0 && isEven((Integer )floor(y))) ? -1 : 1;
  // if y is an integer and a half, use reflection formula
  if (y == floor(y)+0.5)
  {
    // if n+0.5 is a small number return a constant else stirling is accurate
    if (y<50) value = factorialHalvesArray[(Integer )(y)];
    else      value = gammaStirling(y);
  }
  else // normal case
  {
    if (y <= 8)
    {
      Integer n = (Integer )y; // convert y to Integer value
      Real r = y-n;             // r in [0,1]
      value = gammaLanczos(r);  // use Lanzcos approximation
      // scale result
      for (Real i=0.; i<n; i+=1.) value *= (r+i);
    }
    else
    {
      if (y <=16)
      {
        Real r = y+6.;
        value = gammaStirling(r);
        for (Real i=5.; i>=0.; i-=1.) value /= (y+i);
      }
      else
        value = gammaStirling(y);
    }
  }
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflection formula
  Real den = y*sin(Const::_PI_ *y)*value;
  // overflow
  if (den == 0.0) return signgam*Arithmetic<Real>::infinity();
  // normal case
  return -Const::_PI_/den;
}

Real STK::Funct::gammaLn ( Real const & z )

This function computes the function $\ln(\Gamma(z))$ .

The log gamma function is valid when z is non zero nor a negative integer.

if |z| <16 we use the gamma Lanczos method, else we use the Stirling approxiamtion method.

Parameters:

z	given value for the gamma function

Definition at line 821 of file STK_Funct_gamma.cpp.

References STK::Const::_LNSQRTPI_2_, STK::Const::_PI_, STK::abs(), STK::Funct::factorialLnArray, STK::Funct::factorialLnHalvesArray, STK::Funct::gamma(), STK::Funct::gammaLnStirling(), and STK::Funct::stirlingSerie().

Referenced by STK::Funct::factorialLn(), STK::Funct::gammaLnStirlingError(), and STK::Funct::gammaRatio_dl().

{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if (( z<=0 && z == floor(z))||(Arithmetic<Real>::isInfinite(z)))
    throw std::domain_error("Funct::gammaLn(z) "
                       "Null or Negative integer argument");
  // compute the absolute value of x -> y
  Real y = abs(z);
  Real value;
  // If x is an integer
  if (y == floor(z))
  { if (y<=50) value = factorialLnArray[(Integer )y];
    else       value = gammaLnStirling(y);
  }
  // If x is an integer and a half
  if (y == floor(y)+0.5)
  { if (y<=50) value = factorialLnHalvesArray[(Integer )y];
    else       value = gammaLnStirling(y);
  }
  // small values -> use
  if (y <= 16)  return log(abs(gamma(z)));
  else          value = gammaLnStirling(y);
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflectiong formula
  Real sinpiy = abs(sin(Const::_PI_ *y));
  // overflow
  if (sinpiy == 0.0) return -Arithmetic<Real>::infinity();
  // normal case
  return Const::_LNSQRTPI_2_+(z-0.5)*log(y)-z-log(sinpiy)+stirlingSerie(z);
}

Real STK::Funct::gammaLnStirlingError ( Real const & z )

Compute the error when we compute $\ln(\Gamma(z))$ using the Stirling's formula.

Computes the ln of the error term in Stirling's formula.

For z <100, integers or half-integers, use stored values. For z >= 100, uses the stirling serie $ 1/12n - 1/360n^3 + ... $ For other z < 100, uses gammaLn directly (don't use this to write gammaLn!)

$\ln(SE(z)) = \ln(\Gamma(z)) - (z - 1/2) \ln(z) + z - \ln(2\pi)/2$

Parameters:

z	given value for the gamma function

Definition at line 867 of file STK_Funct_gamma.cpp.

References STK::Const::_LNSQRT2PI_, STK::Funct::gammaLn(), STK::Funct::gammaLnStirlingErrorArray, STK::Funct::gammaLnStirlingErrorHalvesArray, and STK::Funct::stirlingSerie().

Referenced by STK::Funct::apois(), STK::Funct::betaRatio_cf(), STK::Funct::betaRatio_sr(), STK::Funct::betaRatio_up(), and STK::Funct::poisson_pdf_raw().

{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if ( z<=0 || Arithmetic<Real>::isInfinite(z))
    throw std::domain_error("Funct::gammaLnStirlingError(z) "
                       "Null or Negative integer argument");
  // z is a discrete value
  if (z == floor(z))
  {
    if (z<100.0) return gammaLnStirlingErrorArray[(Integer )z];
    else         return stirlingSerie(z);
  }
  // z is a discrete value halves
  if (z== floor(z) + 0.5)
  {
    if (z<100.0) return(gammaLnStirlingErrorHalvesArray[(Integer )z]);
    else         return stirlingSerie(z);
  }
  // other cases
  if (z > 16) return stirlingSerie(z);
  // z <16
  return (gammaLn(z) - (z-0.5)*log(z) + z - Const::_LNSQRT2PI_);
}

Real STK::Funct::gammaLnStirlingError ( Integer const & z )

Compute the error when we compute $\ln(\Gamma(z))$ using the Stirling's formula and z is an Integer.

Computes the ln of the error term in Stirling's formula.

For z <100, integers or half-integers, use stored values. For z >= 100, uses the stirling serie $ 1/12n - 1/360n^3 + ... $ For other z < 100, uses gammaLn directly (don't use this to write gammaLn!)

$\ln(SE(z)) = \ln(\Gamma(z)) - (z - 1/2) \ln(z) + z - \ln(2\pi)/2$

Parameters:

z	given value for the gamma function

Definition at line 905 of file STK_Funct_gamma.cpp.

References STK::Funct::gammaLnStirlingErrorArray, and STK::Funct::stirlingSerie().

{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(z)) return Arithmetic<Real>::NA();
  // Negative integer argument not allowed
  if ( z<=0 || Arithmetic<Integer>::isInfinite(z))
    throw std::domain_error("Funct::gammaLnStirlingError(z) "
                       "Null or Negative integer argument");
  // z is a discrete value
  if (z<100.0) return gammaLnStirlingErrorArray[z];
  else         return stirlingSerie(z);
}

void STK::Funct::stirlingCoefficients ( STK::Vector & A )

This function computes the n first coefficients of the Stirling's serie.

The Stirling Coefficients are computed recursively.

The number of coefficient have to be known first : it is given by the size of A.

Parameters:

A	the Vector where the coefficients will be stored, the size of A give the number of coefficients.

Definition at line 925 of file STK_Funct_gamma.cpp.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), and STK::ITContainer1D< TYPE, TContainer1D >::size().

{
  Integer shift = A.first()-1;
  Integer n = A.size();

  // initialization
  for (Integer j=1; j<=n; j++)
    A[shift+j] = pow(4., -j)/(2.*j + 1.);

  // for each coefficients
  Real aux = A[shift+1];
  for (Integer k=2; k<=n; k++)
  {
    // Update remaining terms
    Real coef = (k-0.5) * (k-1.) / 6.;
    for (Integer l=k, r=0; l<=n; l++, r++)
    {
      A[shift+l] -= aux* coef;
      coef *= ((l/(r+2.)) * ((l+0.5)/(r+2.5))) / 4.;
    }
    aux = A[shift+k];
    A[shift+k] /= (2.*k-1.);
  }
}

Real STK::Funct::gammaRatio	(	Real const &	a,
		Real const &	x,
		const bool &	lower_tail
	)

Compute the incomplete gamma functions ratio.

Compute the incomplete gamma function ratio P(a,x)

$P(a, x) = \frac{1}{\Gamma(a)} \int_0^x e^{-t} t^{a-1} dt$

.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 424 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::gammaRatio_ae(), STK::Funct::gammaRatio_cf(), STK::Funct::gammaRatio_dl(), STK::Funct::gammaRatio_sr(), and STK::Funct::poisson_ae().

Referenced by STK::Funct::gammaRatioP(), and STK::Funct::gammaRatioQ().

{
  // Check if a and x are available
  if (Arithmetic<Real>::isNA(a)||Arithmetic<Real>::isNA(x)) return a;
  // Negative parameter not allowed
  if (a<=0)
    throw std::domain_error("Funct::gammaRatio(a,x,lower_tail) "
                       "Negative parameter a");
  // trivial case
  if (x<=0) return lower_tail ? 0 : 1.;
  if (Arithmetic<Real>::isInfinite(x)) return lower_tail ? 1. : 0.;

  // small values of x
  if (x <= 1.) return gammaRatio_dl(a, x, lower_tail);

  // large a compared to x
  if (0.75*a> x) return gammaRatio_sr(a, x, lower_tail);

  // large x compared to a 
  if (0.75*x > a)
  { // if a>1 we can use the a.e. before computing
    if (a>1) return gammaRatio_ae(a, x, lower_tail);
    else     return gammaRatio_cf(a, x, lower_tail);
  }

  // (a<100) and (x ~ a)
  if (a<100.)
  {
    // a greater than x
    if (a>x) return gammaRatio_sr(a, x, lower_tail);
    else
    { // if a>1 we can use the a.e. before computing
      if (a>1) return gammaRatio_ae(a, x, lower_tail);
      else     return gammaRatio_cf(a, x, lower_tail);}
    }
  // (x ~ a) and (a>=100)
  return poisson_ae(a-1, x, lower_tail);
}

Real STK::Funct::gammaRatioQ	(	Real const &	a,
		Real const &	x
	)

Compute the incomplete gamma function ratio Q(a,x).

Compute the incomplete gamma function ratio Q(a,x)

$Q(a, x) = \frac{1}{\Gamma(a)} \int_x^\infty e^{-t} t^{a-1} dt$

.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function

Definition at line 399 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::gammaRatio().

Referenced by STK::Funct::betaRatio_ae().

{
  return gammaRatio(a, x, false);
}

Real STK::Funct::gammaRatioP	(	Real const &	a,
		Real const &	x
	)

Compute the incomplete gamma function ratio P(a,x).

Compute the incomplete gamma function ratio P(a,x)

$P(a, x) = \frac{1}{\Gamma(a)} \int_0^x e^{-t} t^{a-1} dt$

.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function

Definition at line 411 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::gammaRatio().

Referenced by STK::Funct::betaRatio_ae().

{
  return gammaRatio(a, x, true);
}

Real STK::Funct::poisson_pdf_raw	(	Real const &	x,
		Real const &	lambda
	)

Compute the poisson density.

Compute the function:

$p(x, \lambda) = e^{-\lambda} \lambda^x/\Gamma(x+1)$

with good accuracy using the partial deviance.

This is the verision when x is Real.

See also:: http://www.herine.net/stat/software/dbinom.html

Parameters:

x	value to evaluate the function
lambda	value of the parameter

Definition at line 64 of file STK_Funct_poisson_raw.cpp.

References STK::Const::_SQRT2PI_, STK::Funct::dev0(), and STK::Funct::gammaLnStirlingError().

Referenced by STK::Funct::betaRatio_ae(), STK::Funct::gammaRatio_ae(), STK::Funct::gammaRatio_dl(), and STK::Funct::gammaRatio_sr().

{
  // check trivial values
  if (x<0.) return( 0. );
  // if lambda is 0, we have P(X=0) = 1
  if (lambda==0.) return( (x==0.) ? 1. : 0. );
  // special value
  if (x==0.) return( exp(-lambda) );
  // stirling approximation and deviance
  return( exp(-gammaLnStirlingError(x)-dev0(x, lambda))
        / (Const::_SQRT2PI_*sqrt(x))
        );
}

Real STK::Funct::poisson_pdf_raw	(	Integer const &	x,
		Real const &	lambda
	)

Compute the poisson density.

Compute the function:

$p(x, \lambda) = e^{-\lambda} \lambda^x/\Gamma(x+1)$

with good accuracy using the partial deviance.

This is the verision when x is an Integer.

See also:: http://www.herine.net/stat/software/dbinom.html

Parameters:

x	value to evaluate the function
lambda	value of the parameter

Definition at line 91 of file STK_Funct_poisson_raw.cpp.

References STK::Const::_SQRT2PI_, STK::Funct::dev0(), and STK::Funct::gammaLnStirlingError().

{
  // check trivial values
  if (x<0) return( 0. );
  // if lambda is 0, we have P(X=0) = 1
  if (lambda==0) return( (x==0) ? 1. : .0 );
  // special value
  if (x==0) return( exp(-lambda) );
  // stirling approximation and deviance
  return( exp(-gammaLnStirlingError(x)-dev0(x, lambda))
        / (Const::_SQRT2PI_*sqrt(x))
        );
}

Real STK::Funct::erf_raw ( Real const & a )

Compute the error fonction erf(a)

Compute the function

$\mathrm{erf}(a) = \frac{2}{\sqrt{\pi}} \int_0^{a} e^{-t^2} dt$

.

Parameters:

[in] a value to evaluate the function

Definition at line 192 of file STK_Funct_erf_raw.cpp.

References STK::abs(), STK::Funct::p1evl(), STK::Funct::polevl(), STK::Funct::T, and STK::Funct::U.

Referenced by STK::Funct::erfc_raw().

{
  if ( abs(a)> 1.0) return ( 1.0 - erfc(a) );
  Real z = a * a;
  // result
  return ( a * polevl( z, T, 4)/ p1evl( z, U, 5 ));
}

Real STK::Funct::erfc_raw ( Integer const & a )

Compute the complementary error function erfc(a)

Real STK::Funct::normal_cdf_raw ( Real const & x )

Compute the cumulative distribution function of the normal density.

Compute the cumulative distribution function of the normal density

$\mathrm{\Phi}(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-t^2/2} dt = ( 1 + \mathrm{erf}(z) ) / 2 = \mathrm{erfc}(z) / 2$

where $ z = x/sqrt(2) $ .

Computation use functions

erf

and

 erfc

.

Parameters:

x	value to evaluate the function

Definition at line 81 of file STK_Funct_normal_raw.cpp.

References STK::Const::_1_SQRT2_, and STK::abs().

Referenced by STK::Funct::betaRatio_se(), and STK::Funct::poisson_ae().

{
  Real t = x * Const::_1_SQRT2_;
  Real z = abs(t);

  if ( z < Const::_1_SQRT2_) return 0.5 + 0.5 * erf(t);
  Real y = 0.5 * erfc(z);
  return ( t > 0) ? 1. - y : y;
}

Real STK::Funct::normal_pdf_raw ( Real const & x )

compute the probability distribution function of the normal density

$\mathrm{\phi}(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$

Parameters:

x	value to evaluate the function

Definition at line 62 of file STK_Funct_normal_raw.cpp.

References STK::Const::_1_SQRT2PI_.

Referenced by STK::Funct::poisson_ae().

{
  return Const::_1_SQRT2PI_ * exp(-0.5 * x * x);
}

Real STK::Funct::g0 ( Real const & x )

compute the partial deviance $g_0(x) = x\log(x)+ 1 - x$ .

Compute the function:

$g_0(x) = x\log(x)+ 1 - x$

with good accuracy for x near 1.

using the Taylor serie of

$\log\left( \frac{a-x}{a+x}\right)$

x should be >= 0, no check is done.

See also:: http://www.herine.net/stat/software/dbinom.html

Parameters:

[in] x of the partial deviance

Definition at line 61 of file STK_Funct_g0.cpp.

References STK::abs(), and STK::sum().

{
  // special values
  if (x == 0.) return 1.;
  if (x == 1.) return 0.;
  // general case
  Real v = (x-1.)/(x+1.);
  // if 7/9 < x < 9/7  use Taylor serie of log((1+v)/(1-v))
  if (v < 0.125)
  {
    Real sum  = (x-1)*v;
    Real ej   = 2*x*v;
    Real term;
    v = v*v;
    Integer n =0;
    do
    {
      ej *= v;
      sum += (term = (ej/(2*(++n)+1)));
    }
    while (abs(term) > abs(sum) * Arithmetic<Real>::epsilon());
    return sum;
  }
  // else compute directly
  return x*log(double(x))-x + 1.;
}

Real STK::Funct::dev0	(	Real const &	a,
		Real const &	b
	)

compute the partial deviance $d_0(a,b) = a\log(a/b)+ b - a$ .

Compute the partial deviance:

$d_0(a, b) = a\log(a/b)+ b - a$

with good accuracy for a/b is near 1 using the Taylor series of

$\log\left( \frac{a-b}{a+b}\right)$

a should be >=0 and b should be > 0, no check is done.

See also:: http://www.herine.net/stat/software/dbinom.html

Parameters:

a	first parameter
b	second paramater

Returns:: the evaluation of the function d_0

Definition at line 62 of file STK_Funct_dev0.cpp.

References STK::abs(), and STK::sum().

Referenced by STK::Funct::apois(), STK::Funct::betaRatio_cf(), STK::Funct::betaRatio_se(), STK::Funct::betaRatio_sr(), STK::Funct::betaRatio_up(), STK::Funct::poisson_ae(), and STK::Funct::poisson_pdf_raw().

{
  // special values
  if (a == 0.) return b;
  if (b == a)  return 0.;
  // ratio
  Real v = (a-b)/(a+b);
  // small v -> use dl
  if (abs(v)<0.125)
  {
    Real sum  = (a-b)*v;
    Real ej   = 2*a*v;
    Real term;
    v = v*v;
    Integer n =0;
    do
    {
      ej *= v;
      sum += (term = (ej/(2*(++n)+1)));
    }
    while (abs(term) > abs(sum) * Arithmetic<Real>::epsilon());
    return sum;
  }
  // else compute directly the result
  return (a*log(double(a/b))+b-a);
}

Real STK::Funct::log1p ( Real const & x )

compute the fonction $\log(1+x)$ .

When |x| is less than 0.375 the function use a Taylor serie otherwise the default log function is used.

Parameters:

x	value to evaluate the function

Definition at line 88 of file STK_Funct_log1p.cpp.

References STK::abs(), and STK::sumAlternateSerie().

Referenced by STK::Funct::betaRatio_ae().

{
  // trivial values
  if (x == 0.)  return 0.;
  if (x == -1.) return(-Arithmetic<Real>::infinity());
  // check domain
  if (x  < -1)
     throw std::domain_error("Funct::log1p(x) "
                        "Negative Real argument");
  // small values
  if (abs(x) < .375)
  {
    // create functor and compute the alternate serie
    Serielog1p f(x);
    return x * (1. - x * sumAlternateSerie(f));
  }
  // large values : compute directly the result
  return log(1. + x);
}

Real STK::Funct::expm1 ( Real const & x )

compute the fonction $\exp(x)-1$ .

Compute the function

$exp(x)-1.$

with accuracy using it's first order Taylor serie when x is near 0.

Improve the result using a Newton Raphson step.

Definition at line 54 of file STK_Funct_expm1.cpp.

References STK::Const::_LN2_, STK::abs(), and STK::sum().

Referenced by STK::Funct::gammaRatio_dl().

{
  // small values of x
  if (abs(x) < Const::_LN2_)
  {
    // a + 1 != 1 -> use compute Taylor serie else use first order
    // Taylor expansion
    if (abs(x) > Arithmetic<Real>::epsilon())
    {
      // Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ...
      Real term = x, sum =x, n =1.;
      do
      {
        sum += (term *= x/(++n));
      }
      while (abs(term) > abs(sum) * Arithmetic<Real>::epsilon()) ;
      
      return sum ;
    }
    else return (x / 2. + 1.) * x;
  }
  else return exp(double(x))-1.;
}

static Real STK::Funct::betaRatio_cf	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	lower_tail = `true`,
		Integer const &	iterMax = `1000`
	)		`[static]`

Compute the continued fraction of the beta function.

Compute the continued fraction:

$\frac{a_1}{b_1+}\frac{a_2}{b_2+}\frac{a_3}{b_3+}\ldots$

with $ a_1 = 1 $ , $a_{n+1} c_n d_n /((a+2n-2)(a+2n))$ and $b_{n+1} = (1-(a+b)x) + c_n + d_{n+1}$ , where

$c_n = \frac{n(b-n)x}{a+2n-1} \qquad d_n = \frac{(a+n-1)(a+b+n-1)x}{a+2n-1}.$

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
lower_tail	`true` if we want the lower tail, `false` otherwise
iterMax	maximal number of iteration

Returns:: the value of the beta ratio function using the continued fraction method

Definition at line 80 of file STK_Funct_betaRatio.cpp.

References STK::Const::_1_SQRT2PI_, STK::abs(), STK::Funct::dev0(), STK::Funct::gammaLnStirlingError(), and STK::sum().

{
  // Constants
  Real sum = a+b, y = (0.5 -x) + 0.5, res;
  /* Compute the Function:
   *  B(a,b,x,y) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} x^{a} (1-x)^{b}
   */
  Real bt = ( Const::_1_SQRT2PI_*sqrt((a*b)/sum)
        * exp( ( gammaLnStirlingError(sum)
                -gammaLnStirlingError(a)
                -gammaLnStirlingError(b)
                )
              -(dev0(a, sum*x) + dev0(b, y*sum))
              )
        );
  if (bt == 0.) return lower_tail ? 0. : 1.;

  // initialize b_n
  Real bn = (1+a*y-b*x)/(a+1); // b_1
  // initialize numerator
  Real Nold = 0, Ncur = 1.; // = a_1 = 1
  // initialize denominator
  Real Dold = 1, Dcur = bn; // = b_1

  // normalize if necessary
  if (abs(Dcur) > 4294967296.0)
  {
    Ncur /= Dcur;
    Dold /= Dcur;
    Dcur = 1.;
  }

  // auxiliary constant variable
  const Real cte= (0.5-sum*x)+0.5;
  // dn (can be reused)
  Real dn = (a*sum*x)/(a+1);
  // result
  Real cf = 1/Dcur;
  // auxiliary variables
  Real ap2n = a+2, apn = a+1;
  // iterations
  for (Integer n=1; /*n<=iterMax*/; n++, apn++, ap2n += 2)
  {
    // compute auxiliaries variables
    Real nx = n*x, cn = n*(b*x-nx)/(ap2n-1);
    // compute a_n
    Real an = (cn/ap2n)*(dn/(ap2n-2));
    // check trivial case)
    if (an == 0.) break;
    // update d_n
    dn = apn*(sum*x+nx)/(ap2n+1);
    // compute bn
    bn = cte + cn + dn;
    // compute numerator and denominator
    Real anew = bn * Ncur + an * Nold;
    Nold = Ncur;
    Ncur = anew;
    anew = bn * Dcur + an * Dold;
    Dold = Dcur;
    Dcur = anew;
    // normalize if necessary
    if (abs(Dcur) > 4294967296.0)
    {
      Ncur /= Dcur;
      Nold /= Dcur;
      Dold /= Dcur;
      Dcur = 1.;
    }
    // normalize if necessary with 2^32
    if (abs(Ncur) > 4294967296.0)
    {
      Ncur /= 4294967296.0;
      Nold /= 4294967296.0;
      Dold /= 4294967296.0;
      Dcur /= 4294967296.0;
    }
    // test D_n not too small
    if (abs(Dcur) != 0)
    {
      Real cfold = cf;
      cf = Ncur/Dcur;
      // check cv
      if (abs(cf - cfold) < abs(cf)*Arithmetic<Real>::epsilon())
      {
        res = lower_tail ? bt*cf/a : 1-bt*cf/a;
        break;
      }
    }
  }
#ifdef STK_DEBUG
  // result
  res = lower_tail ? bt*cf/a : 1-bt*cf/a;
  if (!Arithmetic<Real>::isFinite(res))
  {
    std::cout << "in cf : bt= " << bt
              << " cf = " << cf
              << " a  = " << a
              << "lt =  " << lower_tail
              << "\n";
  }
  return res;
#else
  return lower_tail ? bt*cf/a : 1-bt*cf/a;
#endif

}

/* @ingroup Analysis
 *  @brief compute the coefs of the beta Ratio function in
 *  its asymptotic expansion.
 *

static Real STK::Funct::betaRatio_ae	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	xm1,
		bool	lower_tail
	)		`[inline, static]`

Compute the incomplete beta function ratio I_x(a,b) using its asymptotic expansion.

Compute the incomplete beta function ratio I_x(a,b) using it' asymptotic expansion.

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
xm1	true if we want to evaluate the function at 1-x
lower_tail	`true` if we want the lower tail, `false` otherwise

Returns:: the value of the beta ratio function using its asymptotic expansion

Definition at line 256 of file STK_Funct_betaRatio.cpp.

References d2, d3, d4, d5, d6, d7, d8, STK::Funct::gammaRatioP(), STK::Funct::gammaRatioQ(), STK::Funct::log1p(), and STK::Funct::poisson_pdf_raw().

Referenced by STK::Funct::betaRatio_up().

{
  // Compute b-1
  Real bm1 = b-1;
  // compute \nu = a+(b-1)/2
  Real nu = a+bm1/2.;
  // compute D = -\nu*log(x)
  Real D = (xm1 ? -log1p(-x) : -log(x)) * nu;
  // check trivial case
  if (D == 0.) return lower_tail ? 1. : 0.;
  if (Arithmetic<Real>::isInfinite(D)) return lower_tail ? 0. : 1.;
  // compute 12*\nu^2
  nu *= 12*nu;
  // variables for the series
  Real term = 1.;
  Real ak_term = 1.;
  Real a_k = 1.;
  // numerator and denominator
  Real den = term;
  Real num = 0.;
  // update term
  term *= bm1;
  // Generalized bernouilli coefs
  // coef 1 (d_1 = 1/2)
  term *= (b)*(b+1)/nu;
  Real coef = term/2; // d_1 = 1/2
  den += coef;
  a_k += (ak_term *= D/(b+1));
  num += coef * a_k;
  // coef 2
  term *= (b+2)*(b+3)/nu;
  coef = d2(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+2));
  a_k += (ak_term *= D/(b+3));
  num += coef * a_k;
  // coef 3
  term *= (b+4)*(b+5)/nu;
  coef = d3(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+4));
  a_k += (ak_term *= D/(b+5));
  num += coef * a_k;
  // coef 4
  term *= (b+6)*(b+7)/nu;
  coef = d4(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+6));
  a_k += (ak_term *= D/(b+7));
  num += coef * a_k;
  // coef 5
  term *= (b+8)*(b+9)/nu;
  coef = d5(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+8));
  a_k += (ak_term *= D/(b+9));
  num += coef * a_k;
  // coef 6
  term *= (b+10)*(b+11)/nu;
  coef = d6(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+10));
  a_k += (ak_term *= D/(b+11));
  num += coef * a_k;
  // coef 7
  term *= (b+12)*(b+13)/nu;
  coef = d7(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+12));
  a_k += (ak_term *= D/(b+13));
  num += coef * a_k;
  // coef 8
  term *= (b+14)*(b+15)/nu;
  coef = d8(bm1)*term;
  den += coef;
  a_k += (ak_term *= D/(b+14));
  a_k += (ak_term *= D/(b+15));
  num += coef * a_k;
  // result (P or Q ?)
  return lower_tail ? gammaRatioQ(b, D)
                    + (num / den) * poisson_pdf_raw(b, D)
                    : gammaRatioP(b, D)
                    - (num / den) * poisson_pdf_raw(b, D);
}

static Real STK::Funct::serie_up	(	Real const &	s,
		Real const &	a,
		Real const &	x,
		Integer const &	n
	)		`[inline, static]`

Compute the recurrence formula of the incomplete beta ratio function.

Compute the recurrence formula of the incomplete beta ratio function I_x(a,b).

Parameters:

s	a+b
a	first parameter
x	value to evaluate the function
n	number of step

Definition at line 358 of file STK_Funct_betaRatio.cpp.

References STK::sum().

Referenced by STK::Funct::betaRatio_up().

{
  // initilize d0 and sum
  Real sum = 1./a;
  Real di  = 1./a;
  // compute \sum_{i=0}^{n-1} \frac{d_i x^i}{a}
  for (Integer i=1; i<n; i++)
    sum += (di *= x*(s-1+i)/(a+i));
  // return result
  return sum;
}

static Real STK::Funct::betaRatio_up	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	xm1,
		bool	lower_tail
	)		`[inline, static]`

Compute the incomplete beta function ratio I_x(a,b) using its recurrence formula and its asymptotic expansion.

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
xm1	true if we want to evaluate the function at 1-x
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 387 of file STK_Funct_betaRatio.cpp.

References STK::Const::_1_SQRT2PI_, STK::Funct::betaRatio_ae(), STK::Funct::dev0(), STK::Funct::gammaLnStirlingError(), and STK::Funct::serie_up().

{
  // number of iterations
  Integer n = Integer(a);
  // compute residual
  Real a0 = a-n;
  // zero case
  if (!a0) { --n; a0=1.;}
  // sum
  Real s= a0+b;
  // compute x and 1-x
  Real x0, y;
  if (xm1)
  { x0 = (0.5 -x) + 0.5; y =x;}
  else
  { x0 = x; y = (0.5 -x) + 0.5;}
  /* Compute the constant
   *  B(a,b,x,y) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} x^a (1-x)^b
   * using saddle-point expansion
   */
  Real bt = ( Const::_1_SQRT2PI_*sqrt((a0*b)/s)
            * exp( ( gammaLnStirlingError(s)
                    -gammaLnStirlingError(a0)
                    -gammaLnStirlingError(b)
                   )
                   -(dev0(a0, s*x0) + dev0(b, y*s))
                 )
            );
  // check trivial case
  if (!bt)
    return betaRatio_ae(b, a0, x, !xm1, !lower_tail);
  // compute serie
    bt *= serie_up(s, a0, x0, n);
  // general case
  return  lower_tail ? betaRatio_ae(b, a0, x, !xm1, !lower_tail) - bt
                     : betaRatio_ae(b, a0, x, !xm1, !lower_tail) + bt;
}

static Real STK::Funct::betaRatio_sr	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	lower_tail
	)		`[static]`

Compute the incomplete beta function ratio I_x(a,b) using its series representation.

Compute the incomplete beta function ratio I_x(a,b) using it's serie representation.

$\frac{x^a}{B(a,b)} \left( \frac{1}{a} + \sum_{n=1}^{\infty} \frac{(1-b)(2-b)\ldots(n-b)}{n! (a+n)} \right)$

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
lower_tail	`true` if we want the lower tail, `false` otherwise

Returns:: the value of the beta ratio function using its series representation

Definition at line 448 of file STK_Funct_betaRatio.cpp.

References STK::Const::_1_SQRT2PI_, STK::abs(), STK::Funct::dev0(), STK::Funct::gammaLnStirlingError(), and STK::sum().

{
  // Constant
  Real sum = a+b;
  // compute B(a,b,x,y) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} x^{a}
  Real bt = ( Const::_1_SQRT2PI_*sqrt((a*b)/sum)
            * exp( sum*x
                 + ( gammaLnStirlingError(sum)
                    -(gammaLnStirlingError(a)+gammaLnStirlingError(b))
                   )
                 -(dev0(a, sum*x) + dev0(b, sum))
                 )
            );
  // check factor
  if (!bt) return lower_tail ? 0. : 1.;
  // parameters
  Real term, n= 1., c = 1.;
  // sum = \sum_{n=1}^\infty (1-b/n) x^n /(a+n)
  sum = 0.0;
  do
  {
    sum += (term =  (c *= (1-b/n)*x) / (a+n));
    ++n;
  }
  while (abs(term) > abs(sum) * Arithmetic<Real>::epsilon());
  // return result
  return lower_tail ? bt*(1/a + sum)
                    : (1-bt/a)-(bt*sum);
}

static Real STK::Funct::coefs_odd_se	(	Real const &	std,
		Real const &	qmp,
		Vector &	A
	)		`[static]`

compute the odd coefs of the beta Ratio function serie expansion.

Given the 2t first coefs, compute and return the (2t+1)-th coefficient. The value is stored to the back of the Vector.

Parameters:

A	vector of dimension 1:n, n>=2, of the coefs
std	the binomial standard deviation $\sqrt{pq}$
qmp	the difference q-p

Definition at line 493 of file STK_Funct_betaRatio.cpp.

References STK::ITContainer1D< TYPE, TContainer1D >::push_back(), and STK::sum().

Referenced by STK::Funct::betaRatio_se().

  {
    // compute result
    Real sum = (qmp/3+std)*(qmp/3-std)/(4*std);
    // save it in A
    A.push_back(sum);
    // and return it
    return sum;
  }
  // initialize the sums
  Real sum, sum1 =0, sum2 =0;
  // update sum
  for (Integer k=1, k1=2; k<t; k++, k1++)
  {
    sum1 += A[k1]*A[n+1-k];
    sum2 += A[k] *A[n-k];
  }
  // compute the result
  sum = ( ((qmp*A[n] - A[t]*A[t])/2 - sum2)/(t+1)
          - sum1- A[t+1]*A[t+1]/2
        )/std;
  // save it in A
  A.push_back(sum);
  // and return it
  return sum;
}

static Real STK::Funct::coefs_even_se	(	Real const &	std,
		Real const &	qmp,
		Vector &	A
	)		`[static]`

compute the even coefs of the beta Ratio function serie expansion.

Given the 2t+1 first coefs, compute and return the (2t+2)-th coefficient. The value is stored to the back of the Vector.

Parameters:

A	vector of dimension 1:n, n>=3, of the coefs
std	the binomial standard deviation $\sqrt{pq}$
qmp	the difference q-p

Definition at line 537 of file STK_Funct_betaRatio.cpp.

References STK::sum().

Referenced by STK::Funct::betaRatio_se().

  {
    // compute result
    Real sum = -qmp*(1+A[3]/std)*2/15;
    // save it in A
    A.push_back(sum);
    // and return it
    return sum;
  }
  // initialize the sums
  Real sum, sum1 =0, sum2 =0;
  // update sum
  for (Integer k=1, k1=2; k<=t; k++, k1++)
  {
    sum1 += A[k1]*A[n+1-k];
    sum2 += A[k] *A[n-k];
  }
  // compute the result
  sum = ( (qmp*A[n] - 2*sum2)/(n+2) - sum1)/std;
  // save it in A
  A.push_back(sum);
  // and return it
  return sum;
}

static Real STK::Funct::betaRatio_se	(	Real const &	a,
		Real const &	b,
		Real const &	x,
		bool	xm1,
		bool	lower_tail,
		Integer const &	iterMax = `20`
	)		`[static]`

Compute the incomplete beta function ratio I_x(a,b) using its serie expansion.

Compute the incomplete beta function ratio I_x(a,b) using it's serie expansion.

Parameters:

a	first parameter, must be >0
b	second parameter, must be >0
x	value to evaluate the function
xm1	`true` if x is in [q,p], `false` otherwise
lower_tail	`true` if we want the lower tail, `false` otherwise
iterMax	Maximal number of iteration (default is 20)

Definition at line 581 of file STK_Funct_betaRatio.cpp.

References STK::Const::_1_SQRT2PI_, _T, STK::abs(), STK::Funct::coefs_even_se(), STK::Funct::coefs_odd_se(), STK::Funct::dev0(), dnorm(), STK::Funct::normal_cdf_raw(), pnorm(), STK::IArray1DBase< TYPE, PTRELT, TArray1D >::reserve(), and stk_cout.

{
  // parameters
  Real s  = a+b;
  Real p  = a/s, q  = b/s;
  Real sx = s*x;
  Real sy = s-sx;
  //Real D  = a*log(x/p)+b*log(y/q);
  Real D  = (xm1 ? dev0(a, sy)+dev0(b, sx) : dev0(a, sx)+dev0(b, sy));

  Real z2 = 2 * D;
  Real z  = (((x<p)&&!xm1)||((x>q)&&xm1)) ? -sqrt(z2) : sqrt(z2);

  //  Compute normal cdf
  Real pnorm = lower_tail ?  normal_cdf_raw( z) : normal_cdf_raw(-z);
#ifdef STK_DEBUG
  if (!Arithmetic<Real>::isFinite(pnorm))
  {
    stk_cout << _T(" a= ") << a
             << _T(" b= ") << b
             << _T(" x= ") << x
             << _T(" p= ") << p
             << _T(" q= ") << q
             << _T(" D= ") << D
             << _T("\n");
  }
#endif
  //  Compute normal pdf
  Real dnorm = Const::_1_SQRT2PI_ * exp(-D);
  // check large values of z
  if (!dnorm) return pnorm;
  // epsilon value for the serie
  Real eps = Arithmetic<Real>::epsilon()*pnorm;

  // auxiliary variables
  Real std = sqrt(a)*sqrt(b)/s, qmp = q - p;

  // serie coefs
  Vector A(2);
  A.reserve(iterMax); // reserve enough space
  A[1] = std;
  A[2] = qmp/3;
  // variables for the series
  Real odd_term  = 1;
  Real even_term = 2/sqrt(s);

  Real sum1_term = z, sum1 = z;
  Real sum2_term = 1, sum2 = 1;

  // variables for the numerator and the denominator
  Real num = qmp*even_term/3; // = A[2]*even_term;
  Real den = std;             // = A[1];

  // compute ratio
  Real res = num/den;

  // computation
  for (Integer k=1; k<iterMax; k++)
  {
    // odd_term : (2k+1)!!/(a+b)^k
    odd_term *= (2*k+1)/s;
    // even_term (2*k)!!/(a+b)^{k+1/2}
    even_term *= 2*(k+1)/s;

    // compute a_{2k+1}
    Real a_n = coefs_odd_se(std, qmp, A);
    // check result
    if (Arithmetic<Real>::isFinite(a_n))
    {
      // denominator \sum_0^k (2j+1)!!a_{2j+1}/(a+b)^j
      den += odd_term * a_n;
      // num+=\sum_1^k (2k+1)!!a_{2k+1}/(a+b)^k (\sum_1^k z^{2j-1}/(2j-1)!!)
      num += odd_term * a_n * sum1;

      // compute a_{2k+2}
      a_n = coefs_even_se(std, qmp, A);
      // check result
      if (Arithmetic<Real>::isFinite(a_n))
      {
        // update sum2_term = z^{2k}/(2k)!!
        // update sum2      = \sum_0^k z^{2j}/(2j)!!
        sum2 += (sum2_term *= z2/(2*k));
        // num+=\sum_0^k (2k+2)!!a_{2k+2}/(a+b)^k (\sum_1^k z^{2j}/(2j)!!)
        num += even_term * a_n * sum2;
      }
    }
    // test cv
    Real old_res = res;
    res = num/den;
    if (abs(old_res - res)*dnorm < eps)
    { break;}

    // update sum1_term = z^{2k+1}/(2k+1)!!
    // update sum1      = \sum_1^{k+1} z^{2j-1}/(2j-1)!!
    sum1 += (sum1_term *= z2/(2*k+1));
  }
#ifdef STK_DEBUG
  // result
  Real result = lower_tail ? pnorm - res * dnorm
                           : pnorm + res * dnorm;
  if (!Arithmetic<Real>::isFinite(result))
  {
    std::cout << "Error in betRatio_se\n"
              << "pnorm = " << pnorm
              << " res = " << res
              << " dnorm  = " << dnorm
              << "lt =  " << lower_tail
              << "\n";
  }
  return result;
#else
  // result
    return lower_tail ? pnorm - res * dnorm
                      : pnorm + res * dnorm;
#endif
}

Real STK::Funct::erfc_raw ( Real const & a )

Compute the function

$\mathrm{erfc}(a) = \frac{2}{\sqrt{\pi}} \int_a^{+\infty} e^{-t^2} dt$

.

Parameters:

[in] a value to evaluate the function

Definition at line 157 of file STK_Funct_erf_raw.cpp.

References STK::abs(), STK::Funct::erf_raw(), STK::Funct::P, STK::Funct::p1evl(), STK::Funct::polevl(), STK::Funct::Q, STK::Funct::R, and STK::Funct::S.

{
  Real p, q, y;

  Real x = abs(a);
  if ( x < 1.0)  return ( 1.0 - erf_raw(a) );

  Real z = exp(-a * a);

  if ( x < 8.0)
  {
    p = polevl( x, P, 8);
    q = p1evl( x, Q, 8);
  }
  else
  {
    p = polevl( x, R, 5);
    q = p1evl( x, S, 6);
  }
  y = (z * p)/q;

  if ( a < 0)  y = 2.0 - y;
  // result
  return (y);
}

static Real STK::Funct::lanczosSerie ( Real const & z ) [static]

Compute the Lanzcos correction serie for the gamma function with n = 21 terms.

Parameters:

z	given value for the lanzcos Serie

Definition at line 592 of file STK_Funct_gamma.cpp.

References STK::Funct::lanczosCoefArray, and STK::sum().

Referenced by STK::Funct::gammaLanczos().

{
  Real sum = 0.0;
  for (Integer k=20; k>=0; k--)
    sum += lanczosCoefArray[k]/(z+k);
  return 2.0240434640140357514731512432760e-10 + sum;
}

static Real STK::Funct::gammaLanczos ( Real const & z ) [static]

Compute the gamma function using the Lanzcos expansion using n = 21 terms and r= 22.618910.

Parameters:

z	given value for the gamma function

Definition at line 605 of file STK_Funct_gamma.cpp.

References STK::Funct::lanczosSerie().

Referenced by STK::Funct::gamma().

{
  // 2 * sqrt(e/pi) = 1.86038273...
  return 1.8603827342052657173362492472666631120594218414085774528900013
       * exp((z-0.5)*(log(z+22.118910)-1.))
       * lanczosSerie(z);
}

static double STK::Funct::stirlingSerie ( Real const & z ) [static]

Compute the Stirling's serie for the gammaLn function.

Parameters:

z	given value for the stirling Serie

Definition at line 617 of file STK_Funct_gamma.cpp.

References STK::Funct::stirlingCoefArray.

Referenced by STK::Funct::gammaLn(), STK::Funct::gammaLnStirling(), STK::Funct::gammaLnStirlingError(), and STK::Funct::gammaStirling().

{
  // use stirling formula
  const Real z2 = z * z;
  if (z <= 50)
    return ( ( stirlingCoefArray[0]
             + ( stirlingCoefArray[1]
               + ( stirlingCoefArray[2]
                 + stirlingCoefArray[3]/z2
                 )/z2
               )/z2
             )/z
           );
    return ( ( stirlingCoefArray[0]
             + ( stirlingCoefArray[1]
                 + stirlingCoefArray[2]/z2
               )/z2
             )/z
           );
}

static Real STK::Funct::gammaStirling ( Real const & z ) [static]

This function computes the gamma function using the Stirling approximation.

This approximation is valid for large values of z.

Parameters:

z	given value for the gamma function

Definition at line 645 of file STK_Funct_gamma.cpp.

References STK::Const::_SQRT2PI_, and STK::Funct::stirlingSerie().

Referenced by STK::Funct::gamma().

{
  return ( Const::_SQRT2PI_
         * exp((z-0.5)*(log(z)-1.)+stirlingSerie(z)-0.5)
         );
}

static Real STK::Funct::gammaLnStirling ( Real const & z ) [static]

This function computes the log gamma function using the Stirling approximation.

This approximation is valid for large values of z.

Parameters:

z	given value for log gamma function

Definition at line 659 of file STK_Funct_gamma.cpp.

References STK::Const::_LNSQRT2PI_, and STK::Funct::stirlingSerie().

Referenced by STK::Funct::gammaLn().

{
  return ( Const::_LNSQRT2PI_ + (z-0.5)*log(z) - z + stirlingSerie(z) );
}

static Real STK::Funct::apois	(	Real const &	a,
		Real const &	b
	)		`[static]`

Compute the poisson density up to a factor.

Compute the function:

$R(a, b) = e^{-b} \frac{b^a}{\Gamma(a)} = a p(a, b)$

where p(a,b) denote the poisson density. We use the partial deviance.

See also:: http://www.herine.net/stat/software/dbinom.html b should be greater than zero.

Parameters:

a	first parameter
b	second parameter

Definition at line 64 of file STK_Funct_gammaRatio.cpp.

References STK::Const::_1_SQRT2PI_, STK::Funct::dev0(), and STK::Funct::gammaLnStirlingError().

Referenced by STK::Funct::gammaRatio_cf(), and STK::Funct::gammaRatio_dl().

{
  // special value
  if (b==0.) return( 0. );
  // stirling approximation and deviance
  return( sqrt(a)*exp(-gammaLnStirlingError(a)-dev0(a, b))
        * (Const::_1_SQRT2PI_)
        );
}

static Real STK::Funct::gammaRatio_dl	(	Real const &	a,
		Real const &	x,
		const bool &	lower_tail
	)		`[static]`

Compute the incomplete gamma function ratio Q(a,x) using the Taylor serie development representation

$Q(a, x) = 1- \frac{x^a}{\Gamma(a+1)} - \frac{x^a}{\Gamma(a)} \sum_{n=1}^\infty (-1)^n\frac{x^n}{(a+n)n!}$

.

See also:: Abramovitz and Stegun [6.5.29] This representation is well behaved for x<=1.

In order to avoid loss of significant digit near 1, we evaluate $1- \exp\{a\ln(x)-\ln(\Gamma(a+1))\}$ using the function expm1 if we want to compute the upper tail.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 124 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::apois(), STK::Funct::expm1(), STK::Funct::gammaLn(), STK::Funct::poisson_pdf_raw(), and STK::sumAlternateSerie().

Referenced by STK::Funct::gammaRatio().

{
  Seriedl f(a, x);
  if (lower_tail)
     return (1.- a*sumAlternateSerie(f))*poisson_pdf_raw(a, x)*exp(x);
  else
    return apois(a, x) * sumAlternateSerie(f) * exp(x)
         - expm1(a*log(x)-gammaLn(a+1.));
}

static Real STK::Funct::gammaRatio_cf	(	Real const &	a,
		Real const &	x,
		const bool	lower_tail
	)		`[static]`

Compute the incomplete gamma function ratio Q(a, x) using its continued fraction representation.

Compute the incomplete gamma function ratio Q(a, x) using the continued fraction representation

$Q(a, x) = \frac{e^{-x} x^a}{\Gamma(a)} \left( \frac{1}{x+1-a-} \frac{1.(1-a)}{x+3-a-} \frac{2.(2-a)}{x+5-a-} \ldots \right)$

This cf is well behaved for almost all value of a and x.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 154 of file STK_Funct_gammaRatio.cpp.

References STK::abs(), and STK::Funct::apois().

Referenced by STK::Funct::gammaRatio(), and STK::Funct::gammaRatio_ae().

{
  // initialize b_n
  Real bn = x + 1. -a; // b_1
  // initialize numerator
  Real Nold = 0, Ncur = 1.; // = a_1 = 1
  // initialize denominator
  Real Dold = 1, Dcur = bn; // = b_1

  // normalize if necessary by 2^32
  if (abs(Dcur) > 4294967296.0)
  {
    Ncur /= Dcur;
    Dold /= Dcur;
    Dcur = 1.;
  }

  // auxiliary variables
  Real aminusn = a;
  // result
  Real cf = Ncur/Dcur;
  // iterations
  for (Integer n=1; ; n++)
  {
    // compute a_n
    Real an = n*(--aminusn);
    // check trivial case
    // (note a_n =0 => b_n != 0 for x>0 and thus cf!=N_n/D_n exists)
    if (an == 0.) break;
    // update b_n
    bn += 2.;
    // compute numerator and denominator
    Real anew = bn * Ncur + an * Nold;
    Nold = Ncur;
    Ncur = anew;
    anew = bn * Dcur + an * Dold;
    Dold = Dcur;
    Dcur = anew;
    // normalize if necessary with 2^32
    if (abs(Dcur) > 4294967296.0)
    {
      Ncur /= Dcur;
      Nold /= Dcur;
      Dold /= Dcur;
      Dcur = 1.;
    }
    // normalize if necessary with 2^32
    if (abs(Ncur) > 4294967296.0)
    {
      Ncur /= 4294967296.0;
      Nold /= 4294967296.0;
      Dold /= 4294967296.0;
      Dcur /= 4294967296.0;
    }
    // test D_n not too small
    if (abs(Dcur) != 0)
    {
      Real cfold = cf;
      cf = Ncur/Dcur;
      // check cv: relative; absolute for small cf
      if (abs(cf - cfold) < abs(cf)*Arithmetic<Real>::epsilon())
        break;
    }
  }
  // cv
  return (lower_tail) ? (1 - apois(a, x)*cf)
                      : (apois(a, x)*cf);
}

static Real STK::Funct::gammaRatio_sr	(	Real const &	a,
		Real const &	x,
		const bool &	lower_tail
	)		`[static]`

Compute the incomplete gamma function ratio P(a,x) using the serie development representation.

Compute the incomplete gamma function ratio P(a,x) using its serie development representation

$P(a, x) = \frac{e^{-x} x^a}{\Gamma(a+1)} \sum_{n=0}^\infty \frac{1}{(a+1) \ldots (a+n)} x^n$

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 240 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::poisson_pdf_raw(), and STK::sum().

Referenced by STK::Funct::gammaRatio(), and STK::Funct::gammaRatio_ae().

{
  Real y = a+1., term = x / (a+1.),  sum = term;
  /* sum =  1+\sum_{n=1}^\infty x^n / (a+1)*(a+2)*...*(a+n))
   */
  do
  {
    y++;
    sum += (term *= x / y);
  }
  while (term > sum * Arithmetic<Real>::epsilon());
  sum+=1;

  if (lower_tail)
    return poisson_pdf_raw(a, x)*sum;
  else
    return 1.-poisson_pdf_raw(a, x)*sum;
}

static Real STK::Funct::gammaRatio_ae	(	Real const &	a,
		Real const &	x,
		const bool &	lower_tail
	)		`[static]`

Compute the incomplete gamma function ratio Q(a,x) using its asymptotic expansion.

Compute the incomplete gamma function ratio Q(a,x) using its asymptotic expansion

$Q(a, x) = \frac{e^{-x} x^{a-1}}{\Gamma(a)} \left[ 1+ \sum_{n=1}^{N-1} \frac{(a-1)\ldots (a-n)}{x^n} +\frac{(a+1) \ldots (a+N)\theta_N}{x^N} \right]$

with $\theta_N < \frac{x}{x-(a_N)}$ . We must have a>1 and x>=a : no check is done.

Parameters:

a	parameter of the gamma ratio function
x	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 282 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::gammaRatio_cf(), STK::Funct::gammaRatio_sr(), STK::Funct::poisson_pdf_raw(), and STK::sum().

Referenced by STK::Funct::gammaRatio().

{
  Real term = 1, sum = 0, b=a-1;
  // sum =  1+\sum_{n=1}^\infty (a-1)*...*(a - n) / x^n
  while ((b > 1.) && (term > sum * Arithmetic<Real>::epsilon()))
  {
    sum += (term *= b / x);
    b--;
  }
  sum+=1.;
  // achieved cv
  if (b>1.)
    return lower_tail ? (1.-poisson_pdf_raw(a-1, x)*sum)
                      : (   poisson_pdf_raw(a-1, x)*sum);
  // case (b>x) and (b <=1)
  if (b>x)
    return lower_tail ? (gammaRatio_sr(b, x, lower_tail)-poisson_pdf_raw(a-1, x)*sum)
                      : (gammaRatio_sr(b, x, lower_tail)+poisson_pdf_raw(a-1, x)*sum);
  // case (b <x) and (b <=1)
  return lower_tail ? gammaRatio_cf(b, x, lower_tail)-poisson_pdf_raw(a-1, x)*sum
                    : gammaRatio_cf(b, x, lower_tail)+poisson_pdf_raw(a-1, x)*sum;
}

static Real STK::Funct::poisson_ae	(	Real const &	a1,
		Real const &	apd,
		const bool &	lower_tail = `true`
	)		`[inline, static]`

Compute the incomplete gamma function ratio P(a,x) using the Poisson asymptotic expansion.

Compute the incomplete gamma function ratio P(a,x) using the asymptotic expansion of the poisson distribution

$\Gamma(a+1, a+d)$

Parameters:

a1	parameter (a-1) of the gamma ratio function
apd	value to evaluate the gamma ratio function
lower_tail	`true` if we want the lower tail, `false` otherwise

Definition at line 321 of file STK_Funct_gammaRatio.cpp.

References STK::Funct::dev0(), STK::Funct::normal_cdf_raw(), and STK::Funct::normal_pdf_raw().

Referenced by STK::Funct::gammaRatio().

{
  // odd coefficients
  static const Real coefs_i[8] =
  {
    0., /* placeholder used for 1-indexing */
          2/3.,
         -4/135.,
          8/2835.,
         16/8505.,
      -8992/12629925.,
    -334144/492567075.,
     698752/1477701225.
  };
  // stirling coefficients
  static const Real coefs_s[8] =
  {
    0., /* placeholder used for 1-indexing */
             1/12.,
             1/288.,
          -139/51840.,
          -571/2488320.,
        163879/209018880.,
       5246819/75246796800.,
    -534703531/902961561600.
  };
  // compute D = -x(log(1+d/x) - d/x) = xlog(x/(x+d))+(x+d)-x 
  Real D   = dev0(a1, apd);
  // compute sqrt(2D)
  Real sqrt2D = sqrt (2. * D);
  // compute D/x
  Real D_x = D/a1;

  // treat negative difference
  if (apd - a1 < 0) sqrt2D = -sqrt2D;

  // variables for the numerator
  Real num = 0;
  Real sum1, num1_term, sum2, num2_term;
  num1_term = sum1 = sqrt(a1);
  num2_term = sum2 = sqrt2D;

  // variables for the denominator
  Real den = a1;
  Real den_term = 1;
  // computation
  for (Integer i = 1; i < 8; i++)
  {
    // compute the numerator
    num += num1_term * coefs_i[i];
    num += num2_term * coefs_s[i];
    // first sum
    sum1 *= (D_x / i) ;  // sqrt(x) * D^n/(n! x^n)
    num1_term = num1_term/a1 + sum1;
    // second sum
    sum2 *= (D_x / (i + 0.5)); // sqrt(2D) * D^n/((1+0.5)*...(n+0.5)*x^n)
    num2_term = num2_term/a1 + sum2;
    // compute the denominator
    den += den_term * coefs_s[i];
    den_term /= a1;
  }
  // result (P or Q ?)
   return lower_tail ? normal_cdf_raw( sqrt2D)
                     - (num / den) * normal_pdf_raw(sqrt2D)
                     : normal_cdf_raw(-sqrt2D)
                     + (num / den) * normal_pdf_raw(sqrt2D);
}

Variable Documentation

const Real STK::Funct::factorialArray[51] [static]

array for the 51th fisrt factorial elements.

The coefficients $n! =n 1 \times 2 \times ... \times n$ . have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; fact.bc

Definition at line 57 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::factorial(), and STK::Funct::gamma().

const Real STK::Funct::factorialLnArray[51] [static]

array for the 51th fisrt ln factorial elements.

The coefficients $\ln(n!) = \ln(1)+ \ln(2)+ \ldots + \ln(n)$ have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; fact.bc

Definition at line 120 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::factorialLn(), and STK::Funct::gammaLn().

const Real STK::Funct::factorialHalvesArray[50] [static]

array for the 51th fisrt halves factorial elements.

The coefficients $\Gamma(n+1/2) = \frac{(2n-1)!!}{2^n} \sqrt(\pi)$ . have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; facthalf.bc

Definition at line 184 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::gamma().

const Real STK::Funct::factorialLnHalvesArray[50] [static]

array for the 51th fisrt halves ln factorial elements.

The coefficients

$ln(\Gamma(n+1/2)) = \ln\left(\frac{(2n-1)!!}{2^n} \right)+\ln(\pi)/2$

. have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; facthalf.bc

Definition at line 249 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::gammaLn().

const Real STK::Funct::gammaLnStirlingErrorArray[100] [static]

array of the gammaLnStirlingError approximation for the values $n=1, 2, 3, \ldots, 99$ .

The value are computed using the formula

$\ln(SE(n)) = \ln(\Gamma(n)) - (n - 1/2) \ln(n) + n - \ln(2\pi)/2$

These coefficient have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; ln_stirl_error.bc

Definition at line 317 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::gammaLnStirlingError().

const Real STK::Funct::gammaLnStirlingErrorHalvesArray[100] [static]

array of the gammaLnStirlingError for the values $z=0.5, 1.5, , \ldots, 99.5$ .

These coefficient have been computed using the infinite precision calculator bc.

See also:: http://www.gnu.org/software/bc/; ln_stirr_error.bc

Definition at line 430 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::gammaLnStirlingError().

const Real STK::Funct::lanczosCoefArray[21] [static]

Initial value:

{
   1.5333183020199267370932516012553e0,  
  -1.1640274608858812982567477805332e1,  
   4.0053698000222503376927701573076e1,  
  -8.2667863469173479039227422723581e1,  
   1.1414465885256804336106748692495e2,  
  -1.1135645608449754488425056563075e2,  
   7.9037451549298877731413453151252e1,  
  -4.1415428804507353801947558814560e1,  
   1.6094742170165161102085734210327e1,  
  -4.6223809979028638614212851576524e0,  
   9.7030884294357827423006360746167e-1, 
  -1.4607332380456449418243363858893e-1, 
   1.5330325530769204955496334450658e-2, 
  -1.0773862404547660506042948153734e-3, 
   4.7911128916072940196391032755132e-5, 
  -1.2437781042887028450811158692678e-6, 
   1.6751019107496606112103160490729e-8, 
  -9.7674656970897286097939311684868e-11,
   1.8326577220560509759575892664132e-13,
  -6.4508377189118502115673823719605e-17,
   1.3382662604773700632782310392171e-21 
}

array of the Lanzcos coefficients.

These coefficient have been found in the thesis of Pugh, Glendon (2004)

See also:: http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf

Definition at line 541 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::lanczosSerie().

const Real STK::Funct::stirlingCoefArray[9] [static]

Initial value:

{
 +0.0833333333333333333333333333333333333333333333333333333333333333333,
 -0.0027777777777777777777777777777777777777777777777777777777777777777,
 +0.000793650793650793650793650793650793650793650793650793650793650,
 -0.0005952380952380952380952380952380952380952380952380952380952380,
 +0.000841750841750841750841750841750841750841750841750841750841750,
 -0.0019175269175269175269175269175269175269175269175269175269175269,
 +0.006410256410256410256410256410256410256410256410256410256410256410,
 -0.029550653594771241830065359477124183006535947712418300653594771241,
 +0.17964437236883057316493849001588939669435025472177
}

array of the Stirling coefficients.

These coefficients have been computed by Horace S. Uhler (1942)

See also:: Horace S. Uhler (1942), The Coefficients of Stirling's Series for log Gamma (x), Proceeding of the National Academy of Science of United States.

Definition at line 574 of file STK_Funct_gamma.cpp.

Referenced by STK::Funct::stirlingSerie().

The Analysis sub-project

Classes

Functions

Variables

Detailed Description

Function Documentation

Variable Documentation