STK++ 0.9.13
Namespaces | Classes | Functions | Variables
Special functions tools

In this project we compute usual and special functions. More...

Namespaces

namespace  STK::Funct
 The namespace Funct enclose all usual and special functions.
 
namespace  STK::Algo
 The namespace Algo enclose all generic algorithms.
 

Classes

class  STK::Funct::Serielog1p
 This Series computes. More...
 
class  STK::IFunction< Function >
 Interface base class for functions. More...
 
class  STK::ISerie< Serie >
 Interface base class for Series. More...
 
class  STK::Funct::Seriedl
 This Serie computes. More...
 

Functions

template<class Serie >
Real STK::sumAlternateSerie (const ISerie< Serie > &f, int const &n=0)
 Sum an alternating series using the Chebichev polynomials.
 
template<class Serie >
Real STK::sumSerie (const ISerie< Serie > &f, int 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, int const &iter_max=100)
 Evaluate a continued fraction.
 
template<class Function >
Real STK::Algo::SecantMethod (IFunction< Function > const &f, Real const &x0, Real const &x1, Real tol)
 apply the secant method for finding the zero of a function.
 
template<class Function >
Real STK::Algo::BrentMethod (IFunction< Function > const &f, Real const &x0, Real const &x1, Real tol)
 perform the brent's algorithm.
 
template<class Function >
Real STK::Algo::findZero (IFunction< Function > const &f, Real const &x0, Real const &x1, Real tol)
 find the zero of a function.
 
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_sr (Real const &a, Real const &b, Real x, bool xm1, bool lower_tail)
 Compute the incomplete beta function ratio I_x(a,b) using its series representation.
 
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.
 
Real STK::Funct::betaRatio_cf (Real const &a, Real const &b, Real x, bool xm1, bool lower_tail=true)
 Compute the incomplete beta function ratio using the continued fraction method.
 
Real STK::Funct::betaRatio_se (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 the serie expansion method.
 
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 the asymptotic expansion method.
 
Real STK::Funct::betaRatio (Real const &a, Real const &b, Real const &x, bool lower_tail=true)
 Compute the incomplete beta function ratio Compute the beta ratio function.
 
Real STK::Funct::lanczosSerie (Real const &z)
 Compute the Lanzcos correction series for the gamma function with n = 21 terms.
 
Real STK::Funct::gammaLanczos (Real const &z)
 Compute the gamma function using the Lanzcos expansion using n = 21 terms and r= 22.618910.
 
double STK::Funct::stirlingSerie (Real const &z)
 Compute the Stirling's series for the lgamma function.
 
Real STK::Funct::gammaStirling (Real const &z)
 This function computes the gamma function using the Stirling approximation.
 
Real STK::Funct::lgammaStirling (Real const &z)
 This function computes the log gamma function using the Stirling approximation.
 
Real STK::Funct::lgammaStirlingError (Real const &z)
 Compute the error when we compute $ \ln(\Gamma(z)) $ using the Stirling's formula.
 
Real STK::Funct::lgammaStirlingError (int n)
 Compute the error when we compute $ \ln(\Gamma(z)) $ using the Stirling's formula and z is an integer.
 
Real STK::Funct::factorial (int n)
 This function computes $ n! $ for integer argument.
 
Real STK::Funct::factorial (Real const &z)
 This function computes $ z! $ when z is an integer in a Real format.
 
Real STK::Funct::lfactorial (int n)
 This function computes $ \ln(n!) $ for integer argument.
 
Real STK::Funct::lfactorial (Real const &z)
 This function computes $ \ln(z!) $ when z is an integer in Real format.
 
Real STK::Funct::gamma (Real const &z)
 This function computes the function $ \Gamma(z) $.
 
Real STK::Funct::lgamma (Real const &z)
 This function computes the function $ \ln(\Gamma(z)) $.
 
Real STK::Funct::gammaRatio (Real const &a, Real const &x, 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::dev0 (Real const &a, Real const &b)
 compute the partial deviance $ d_0(a,b) = a\log(a/b)+ b - a $.
 
Real STK::Funct::b1 (Real const &a, Real const &b, Real const &x, bool xm1)
 Compute the function.
 
Real STK::Funct::g0 (Real const &x)
 compute the partial deviance $g_0(x) = x\log(x)+ 1 - x$.
 
Real STK::Funct::log1p (Real const &x)
 compute the function $ \log(1+x) $.
 
Real STK::Funct::expm1 (Real const &x)
 compute the function $ \exp(x)-1 $.
 
Real STK::Funct::beta_pdf_raw (Real const &x, Real const &a, Real const &b)
 Compute the beta density function.
 
Real STK::Funct::binomial_pdf_raw (Real const &x, Real const &n, Real const &p)
 Compute the generalized binomial probability mass function.
 
Real STK::Funct::binomial_pdf_raw (int x, int n, Real const &p)
 Compute the binomial probability mass function.
 
Real STK::Funct::binomial_lpdf_raw (Real const &x, Real const &n, Real const &p)
 Compute the generalized binomial log-probability mass function.
 
Real STK::Funct::binomial_lpdf_raw (int x, int n, Real const &p)
 Compute the binomial log-probability mass function.
 
Real STK::Funct::poisson_pdf_raw (Real const &x, Real const &lambda)
 Compute the Poisson density.
 
Real STK::Funct::poisson_pdf_raw (int x, Real const &lambda)
 Compute the poisson density with integer value.
 
Real STK::Funct::poisson_lpdf_raw (Real const &x, Real const &lambda)
 Compute the log-poisson density.
 
Real STK::Funct::poisson_lpdf_raw (int x, Real const &lambda)
 Compute the log-poisson density with integer value.
 
Real STK::Funct::erf_raw (Real const &a)
 Compute the error function erf(a) Compute the function.
 
Real STK::Funct::erfc_raw (Real const &a)
 Compute the complementary error function erfc(a) Compute the function.
 
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::psi_raw (Real x)
 Compute the psi function.
 
template<int N>
Real STK::Funct::evalPolynomial (Real x, const Real *P)
 Polynomial evaluator.
 
template<int N>
Real STK::Funct::evalPolynomial1 (Real x, const Real *P)
 Polynomial evaluator.
 

Variables

const Real STK::Const::factorialArray [51]
 array for the 51th fisrt factorial elements.
 
const Real STK::Const::doubleFactorialOddArray [19]
 array for the double factorial of odd numbers.
 
const Real STK::Const::doubleFactorialEvenArray [20]
 array for the double factorial of even numbers.
 
const Real STK::Const::doubleFactorialArray [39]
 array for the double factorial.
 
const Real STK::Const::bernouilliNumbersArray [21]
 Array of the 40th first Bernouilli numbers Bernouilli(n) n=0, 2, 4, ... ,40.
 

Detailed Description

In this project we compute usual and special functions.

The Analysis provide all the tools necessary to the computation of the usual and special functions like gamma, beta, gamma ratio, beta ration functions. It provide generic algorithms and the usual mathematical constants.

Function Documentation

◆ b1()

Real STK::Funct::b1 ( Real const &  a,
Real const &  b,
Real const &  x,
bool  xm1 
)
inline

Compute the function.

\[ B_1(a,b,x) = \frac{ x^{a} (1-x)^{b}}{B(a,b)} \]

using the partial deviance $ (a+b) * (p*log(x/p)+q*log((1-x)/q)) $.

Parameters
a,b,xparameters of the generalized beta
xm1true if x is to be taken as 1-x

Definition at line 103 of file STK_Funct_raw.h.

104{
105 if (x == 0) return 0;
106 if (x == 1) return 0;
107 Real s = a+b, sx = s*x, sy = s*(1.-x);
108 return ( std::exp(- Const::_LNSQRT2PI_
109 + 0.5 * ( a<b ? std::log(a) + log1p(-a/s) : std::log(b) + log1p(-b/s))
110 + (lgammaStirlingError(s)-lgammaStirlingError(a)-lgammaStirlingError(b))
111 - (xm1 ? dev0(a, sy)+dev0(b, sx) : dev0(a, sx)+dev0(b, sy))
112 )
113 );
114}
STK::Funct::dev0
Real dev0(Real const &a, Real const &b)
compute the partial deviance .
Definition STK_Funct_raw.h:74
STK::Funct::log1p
Real log1p(Real const &x)
compute the function .
Definition STK_Funct_raw.h:145
STK::Funct::lgammaStirlingError
Real lgammaStirlingError(Real const &z)
Compute the error when we compute using the Stirling's formula.
Definition STK_Funct_gamma.h:149
STK::Real
double Real
STK fundamental type of Real values.
Definition STK_Typedefs.h:118

References STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), and STK::Funct::log1p().

◆ beta_pdf_raw()

Real STK::Funct::beta_pdf_raw ( Real const &  x,
Real const &  a,
Real const &  b 
)
inline

Compute the beta density function.

Compute the function (beta pdf)

\[ B(a,b,x) = \frac{ x^{a-1} (1-x)^{b-1}}{B(a,b)}. \]

Parameters
a,b,xparameters of the beta density

Definition at line 191 of file STK_Funct_raw.h.

192{
193 // trivial cases
194 if (x < 0 || x > 1) return(0);
195 if (x == 0)
196 {
197 if(a > 1) return(0);
198 if(a < 1) return(Arithmetic<Real>::infinity());
199 return(b); // a == 1
200 }
201 if (x == 1)
202 {
203 if(b > 1) return(0);
204 if(b < 1) return(Arithmetic<Real>::infinity());
205 return(a); // b == 1
206 }
207 // limit cases with x \in (0,1)
208 if (a == 0 || b == 0) { return(0);}
209 // general case
210 Real s = a+b, sx = s*x, sy = s*(1.-x);
211 return ( std::exp(- Const::_LNSQRT2PI_
212 + 0.5 *( a<b ? std::log(a) + log1p(-a/s) : std::log(b) + log1p(-b/s))
213 + lgammaStirlingError(s)-lgammaStirlingError(a)-lgammaStirlingError(b)
214 - dev0(a, sx) - dev0(b, sy)
215 - std::log(x) - log1p(-x)
216 )
217 );
218}

References STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), and STK::Funct::log1p().

◆ betaRatio()

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

Compute the incomplete beta function ratio Compute the beta ratio function.

\[ I_x(a,b) = \frac{\int_0^x u^{a-1} (1-u)^{b-1}}{\int_0^\infty u^{a-1} (1-u)^{b-1}} du \]

for $ 0\leq x \leq 1$.

Parameters
a,bfirst and second parameters, must be >0
xvalue to evaluate the function
lower_tailtrue if we want the lower tail, false otherwise

Definition at line 547 of file STK_Funct_betaRatio.cpp.

548{
549 // Check if a and x are available
550 if ( isNA(a) || isNA(b) || isNA(x) ) return Arithmetic<Real>::NA();
551 // Negative parameter not allowed
552 if ((a<=0)||(b<=0))
553 STKDOMAIN_ERROR_2ARG(Funct::betaRatio,a,b,negative or null parameters not allowed);
554 return betaRatio_raw(a,b,x,lower_tail);
555}
STKDOMAIN_ERROR_2ARG
#define STKDOMAIN_ERROR_2ARG(Where, Arg1, Arg2, Error)
Definition STK_Macros.h:147
STK::isNA
bool isNA(Type const &x)
utility method allowing to know if a value is a NA (Not Available) value
Definition STK_Arithmetic.h:214
STK::Funct::betaRatio_raw
Real betaRatio_raw(Real const &a, Real const &b, Real const &x, bool lower_tail)
Definition STK_Funct_betaRatio.cpp:556

References STK::Funct::betaRatio(), STK::Funct::betaRatio_raw(), STK::isNA(), STK::Arithmetic< Type >::NA(), and STKDOMAIN_ERROR_2ARG.

Referenced by STK::Funct::betaRatio(), STK::Law::Beta::cdf(), and STK::Law::Beta::cdf().

◆ betaRatio_ae()

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 the asymptotic expansion method.

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,bfirst and second parameters, must be >0
xvalue to evaluate the function
xm1true if we want to evaluate the function at 1-x, false otherwise
lower_tailtrue if we want the lower tail, false otherwise
Returns
the value of the beta ratio function using its asymptotic expansion

Definition at line 75 of file STK_Funct_betaRatio.cpp.

78{
79#ifdef STK_BETARATIO_DEBUG
80 stk_cout << _T("BetaRatio_ae(") << a << _T(", ") << b << _T(", ") << x << _T(")\n");
81#endif
82 // Compute b-1
83 Real bm1 = b-1;
84 // compute \nu = a+(b-1)/2
85 Real nu = a+bm1/2.;
86 // compute D = -\nu*log(x)
87 Real D = (xm1 ? -log1p(-x) : -log(x)) * nu;
88 // check trivial case
89 if (D == 0.) return lower_tail ? 1. : 0.;
90 if (Arithmetic<Real>::isInfinite(D)) return lower_tail ? 0. : 1.;
91 // compute 12*\nu^2
92 nu *= 12*nu;
93 // variables for the series
94 Real term = 1.;
95 Real ak_term = 1.;
96 Real a_k = 1.;
97 // numerator and denominator
98 Real den = term;
99 Real num = 0.;
100 // update term
101 term *= bm1;
102 // Generalized bernouilli coefs
103 // coef 1 (d_1 = 1/2)
104 term *= (b)*(b+1)/nu;
105 Real coef = term/2; // d_1 = 1/2
106 den += coef;
107 a_k += (ak_term *= D/(b+1));
108 num += coef * a_k;
109 // coef 2
110 term *= (b+2)*(b+3)/nu;
111 coef = d2(bm1)*term;
112 den += coef;
113 a_k += (ak_term *= D/(b+2));
114 a_k += (ak_term *= D/(b+3));
115 num += coef * a_k;
116 // coef 3
117 term *= (b+4)*(b+5)/nu;
118 coef = d3(bm1)*term;
119 den += coef;
120 a_k += (ak_term *= D/(b+4));
121 a_k += (ak_term *= D/(b+5));
122 num += coef * a_k;
123 // coef 4
124 term *= (b+6)*(b+7)/nu;
125 coef = d4(bm1)*term;
126 den += coef;
127 a_k += (ak_term *= D/(b+6));
128 a_k += (ak_term *= D/(b+7));
129 num += coef * a_k;
130 // coef 5
131 term *= (b+8)*(b+9)/nu;
132 coef = d5(bm1)*term;
133 den += coef;
134 a_k += (ak_term *= D/(b+8));
135 a_k += (ak_term *= D/(b+9));
136 num += coef * a_k;
137 // coef 6
138 term *= (b+10)*(b+11)/nu;
139 coef = d6(bm1)*term;
140 den += coef;
141 a_k += (ak_term *= D/(b+10));
142 a_k += (ak_term *= D/(b+11));
143 num += coef * a_k;
144 // coef 7
145 term *= (b+12)*(b+13)/nu;
146 coef = d7(bm1)*term;
147 den += coef;
148 a_k += (ak_term *= D/(b+12));
149 a_k += (ak_term *= D/(b+13));
150 num += coef * a_k;
151 // coef 8
152 term *= (b+14)*(b+15)/nu;
153 coef = d8(bm1)*term;
154 den += coef;
155 a_k += (ak_term *= D/(b+14));
156 a_k += (ak_term *= D/(b+15));
157 num += coef * a_k;
158 // result (P or Q ?)
159 return lower_tail ? gammaRatioQ(b, D)
160 + (num / den) * poisson_pdf_raw(b, D)
161 : gammaRatioP(b, D)
162 - (num / den) * poisson_pdf_raw(b, D);
163}
d6
#define d6(z)
Definition STK_Funct_betaRatio.cpp:48
d8
#define d8(z)
Definition STK_Funct_betaRatio.cpp:52
d5
#define d5(z)
Definition STK_Funct_betaRatio.cpp:46
d7
#define d7(z)
Definition STK_Funct_betaRatio.cpp:50
d4
#define d4(z)
Definition STK_Funct_betaRatio.cpp:45
d2
#define d2(z)
Definition STK_Funct_betaRatio.cpp:43
d3
#define d3(z)
Definition STK_Funct_betaRatio.cpp:44
stk_cout
#define stk_cout
Standard stk output stream.
Definition STK_Typedefs.h:73
_T
#define _T(x)
Let x unmodified.
Definition STK_Typedefs.h:71
STK::Funct::poisson_pdf_raw
Real poisson_pdf_raw(Real const &x, Real const &lambda)
Compute the Poisson density.
Definition STK_Funct_raw.h:321
STK::Funct::gammaRatioQ
Real gammaRatioQ(Real const &a, Real const &x)
Compute the incomplete gamma function ratio Q(a,x).
Definition STK_Funct_gammaRatio.h:85
STK::Funct::gammaRatioP
Real gammaRatioP(Real const &a, Real const &x)
Compute the incomplete gamma function ratio P(a,x).
Definition STK_Funct_gammaRatio.h:97

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

Referenced by STK::Funct::betaRatio_raw(), and STK::Funct::betaRatio_up().

◆ betaRatio_cf()

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

Compute the incomplete beta function ratio using the continued fraction method.

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}. \]

We are using a new formulation of the

Parameters
a,bfirst and second parameters, must be >0
xvalue to evaluate the function
xm1true if we are looking for 1-x rather than x, false otherwise
lower_tailtrue if we want the lower tail, false otherwise
Returns
the value of the beta ratio function using the continued fraction method

Definition at line 302 of file STK_Funct_betaRatio.cpp.

305{
306#ifdef STK_BETARATIO_DEBUG
307 stk_cout << _T("BetaRatio_cf(") << a << _T(", ") << b << _T(", ") << x << _T(")\n");
308#endif
309 // parameters
310 Real s = a+b, sx = s*x, sy = s-sx;
311 /* Compute the function:
312 * B(a,b,x,y) = \frac{ x^{a} (1-x)^{b}}{B(a,b)}
313 * using (a+b) * (p*log(x/p)+q*log(y/q))
314 */
315 Real bt = b1(a, b, x, xm1);
316 if (bt*Arithmetic<Real>::epsilon() == 0.) return lower_tail ? 0. : 1.;
317
318 // initialize numerator
319 Real Nold = 0, Ncur = 1.; // = a_1 = 1
320 // initialize denominator
321 Real Dold = 1., Dcur = xm1 ? 1.-sy/(a+1.) : 1.-sx/(a+1.); // = b_1
322 // normalize if necessary
323 if (std::abs(Dcur) > 4294967296.)
324 {
325 Ncur /= Dcur;
326 Dold /= Dcur;
327 Dcur = 1.;
328 }
329 // auxiliary constant variables
330 const Real x_4 = xm1 ? (0.125 - x/4.) + 0.125 : x/4.;
331 const Real cte_b1 = xm1 ? x/2. + 0.5 : (0.5 - x/2.) + 0.5;
332 const Real cte_b2 = xm1 ? (a-1.)*(b+s-1.)*((0.25-x/2.)+0.25) : (a-1.)*(b+s-1.)*x/2.;
333 const Real cte_a1 = a*(a-2.)*(b+s)*(b+s-2.), cte_a2 = (a-1.)*(b+s-1.)*x_4;
334 // result
335 Real cf = 1/Dcur;
336 // auxiliary variables
337 Real ap2n = a+2., ap2nm1 = a+1.,ap2nm2 = a;
338 // iterations
339 for (int n=1; /*n<=iterMax*/; n++, ap2n += 2., ap2nm1 += 2., ap2nm2 += 2.)
340 {
341 // compute a_n
342 Real a_n = (cte_a2/ap2nm1 - x_4) * (x_4 + cte_a2/ap2nm1) - cte_a1*(x_4/ap2n)*(x_4/ap2nm2) ;
343 // check trivial case
344 if (a_n*Arithmetic<Real>::epsilon() == 0.) break;
345 // compute b_n
346 Real b_n = cte_b1 - cte_b2/(ap2nm1*(ap2n+1));
347 // compute numerator and denominator
348 Real anew = b_n * Ncur + a_n * Nold;
349 Nold = Ncur;
350 Ncur = anew;
351 // check for underflow
352 if (Ncur*Arithmetic<Real>::epsilon() == 0.) { cf =0.; break;}
353 anew = b_n * Dcur + a_n * Dold;
354 Dold = Dcur;
355 Dcur = anew;
356 // normalize if necessary
357 if (std::abs(Dcur) > 4294967296.)
358 {
359 Ncur /= Dcur;
360 Nold /= Dcur;
361 Dold /= Dcur;
362 Dcur = 1.;
363 }
364 // normalize if necessary with 2^32
365 if (std::abs(Ncur) > 4294967296.)
366 {
367 Ncur /= 4294967296.;
368 Nold /= 4294967296.;
369 Dold /= 4294967296.;
370 Dcur /= 4294967296.;
371 }
372 // test D_n not too small
373 if (std::abs(Dcur) != 0)
374 {
375 Real cfold = cf;
376 cf = Ncur/Dcur;
377 // check cv
378 if (std::abs(cf - cfold) < std::abs(cf)*Arithmetic<Real>::epsilon()) { break;}
379 }
380 }
381 return lower_tail ? bt*cf/a : 1-bt*cf/a;
382}
b1
const double b1
Definition STK_Law_Normal.cpp:48

References _T, b1, and stk_cout.

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

◆ betaRatio_se()

Real STK::Funct::betaRatio_se ( 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 the serie expansion method.

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

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

Parameters
a,bfirst and second parameters, must be >0
xvalue to evaluate the function
xm1true if we want 1-x value, false otherwise
lower_tailtrue if we want the lower tail, false otherwise

Definition at line 486 of file STK_Funct_betaRatio.cpp.

489{
490#ifdef STK_BETARATIO_DEBUG
491 stk_cout << _T("BetaRatio_se(") << a << _T(", ") << b << _T(", ") << x << _T(")\n");
492#endif
493 // parameters
494 Real s = a+b, p = a/s, q = b/s, sx = s*x, sy = s-sx;
495 //Real z2 = 2 * (a+b) (p*log(x/p)+q*log(y/q));
496 Real z2 = 2. * (xm1 ? dev0(a, sy)+dev0(b, sx) : dev0(a, sx)+dev0(b, sy));
497 Real z = (((x<p)&&!xm1)||((x>q)&&xm1)) ? -std::sqrt(z2) : std::sqrt(z2);
498 // Compute normal cdf
499 Real pnorm = lower_tail ? normal_cdf_raw( z) : normal_cdf_raw(-z);
500 // Compute normal pdf
501 Real dnorm = Const::_1_SQRT2PI_ * exp(-z2/2.);
502 // check large values of z
503 if (dnorm < Arithmetic<Real>::epsilon()) return pnorm;
504
505 // auxiliary variables
506 Real a1 = std::sqrt(a)*std::sqrt(b)/s, a2 = (b - a)/(3.*s), a2_a1 = a2/a1, sqrts = std::sqrt(s);
507
508 // series coefs
509 std::vector<Real> A, B;
510 A.reserve(10); // reserve enough space
511 A.push_back(p); // a0
512 A.push_back(a1); // a1
513 A.push_back(a2); // a2
514 A.push_back(a1 * (a2_a1/2.+0.5) * (a2_a1/2.-0.5)); //a3
515 A.push_back(-0.4*a2*((a2_a1/2.+0.5) * (a2_a1/2.-0.5)+1.)); // a4
516 // compute the numerator and the denominator
517 Real num = (a2 * 2.
518 + (A[3] * 3. * z
519 + (A[4] * 4. * (2+ z2)
520 + (coefs_odd_se(A) * 5. * z * (3 + z2)
521 + (coefs_even_se(A) * 6. * (2*4 + z2 * (4. + z2))
522 + (coefs_odd_se(A) * 7. * z *(5*3 + z2*(5 + z2))
523 + (coefs_even_se(A) * 8. * (2*4*6 + z2 * (4*6 + z2 *(6 + z2)))
524 + (coefs_odd_se(A) * 9. * z * (7*5*3 + z2 * (7*5 + z2 * (7+z2)))
525 + (coefs_even_se(A) * 10. * (2*4*6*8 + z2 * (4*6*8 + z2 *(6*8 + z2 * (8 + z2))))
526 + (coefs_odd_se(A) * 11. * z * (9*7*5*3 + z2 * (9*7*5 + z2 * (9*7+z2 * (9 + z2))))
527 )/sqrts)/sqrts)/sqrts)/sqrts)/sqrts)/sqrts)/sqrts)/sqrts)/sqrts)/sqrts;
528 Real den = a1 + (A[3] + (A[5] + (A[7] + (A[9] + A[11] * 11./s) * 9./s)*7./s)*5./s)*3./s;
529
530 // compute ratio
531 Real ratio = num/den;
532 // result
533 return lower_tail ? pnorm - ratio * dnorm : pnorm + ratio * dnorm;
534}
a2
const double a2
Definition STK_Law_Normal.cpp:42
a1
const double a1
Definition STK_Law_Normal.cpp:41
STK::Funct::normal_cdf_raw
Real normal_cdf_raw(Real const &x)
Compute the cumulative distribution function of the normal density.
Definition STK_Funct_raw.h:498

References _T, a1, a2, STK::Funct::dev0(), STK::Funct::normal_cdf_raw(), and stk_cout.

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

◆ betaRatio_sr()

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

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 series 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)} x^n \right) \]

Parameters
a,bfirst and second parameters, must be >0
xvalue to evaluate the function
xm1true if we want to compute the function at 1-x
lower_tailtrue if we want the lower tail, false otherwise
Returns
the value of the beta ratio function using its series representation

Definition at line 180 of file STK_Funct_betaRatio.cpp.

183{
184#ifdef STK_BETARATIO_DEBUG
185 stk_cout << _T("BetaRatio_sr(") << a << _T(", ") << b << _T(", ") << x << _T(")\n");
186#endif
187 // constants
188 Real s = a+b, sx = s*x, sy = s*(1-x);
189 // compute B(a,b,x,y) = \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} x^{a}
190 Real bt = ( Const::_1_SQRT2PI_*sqrt((a*b)/s)
191 * exp( ( lgammaStirlingError(s)
192 -(lgammaStirlingError(a)+lgammaStirlingError(b))
193 )
194 + (xm1 ? sy- (dev0(a, sy) + dev0(b, s)) : sx- (dev0(a, sx) + dev0(b, s)))
195 )
196 );
197 Real y = (0.5 - x) + 0.5;
198 if (xm1) std::swap(x,y);
199 // check factor
200 if (!bt) return lower_tail ? 0. : 1.;
201 // parameters
202 int n = 1;
203 Real sum = 0.0, c= 1., term;
204 do
205 {
206 sum += (term = (c *= (1.-b/n)*x) / (a+n));
207 ++n;
208 }
209 while (std::abs(term) > std::abs(sum) * Arithmetic<Real>::epsilon());
210 // return result
211 return lower_tail ? bt*(1./a + sum) : (1.-bt/a-bt*sum);
212}
STK::sum
Arrays::SumOp< Lhs, Rhs >::result_type sum(Lhs const &lhs, Rhs const &rhs)
convenience function for summing two arrays
Definition STK_Array2D_Functors.h:554

References _T, STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), stk_cout, and STK::sum().

Referenced by STK::Funct::betaRatio_raw(), and STK::Funct::betaRatio_up().

◆ betaRatio_up()

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.

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

Parameters
a,bfirst and second parameters, must be >0
xvalue to evaluate the function
xm1true if we want to evaluate the function at 1-x
lower_tailtrue if we want the lower tail, false otherwise

Definition at line 264 of file STK_Funct_betaRatio.cpp.

267{
268#ifdef STK_BETARATIO_DEBUG
269 stk_cout << _T("BetaRatio_up(") << a << _T(", ") << b << _T(", ") << x << _T(")\n");
270#endif
271 // number of iterations
272 int n = int(a);
273 // compute residual
274 Real a0 = a-n;
275 // zero case
276 if (!a0) { --n; a0=1.;}
277 // sum
278 Real c= serie_up(a0, b, x, xm1, n) * (lower_tail ? - 1. : 1.);
279 return (b<=15) ? betaRatio_sr(a0, b, x, xm1, lower_tail) +c
280 : betaRatio_ae(b, a0, x, !xm1, !lower_tail) +c;
281}
STK::Funct::betaRatio_ae
Real 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 the asymptotic expansion method.
Definition STK_Funct_betaRatio.cpp:75
STK::Funct::betaRatio_sr
Real betaRatio_sr(Real const &a, Real const &b, Real x, bool xm1, bool lower_tail)
Compute the incomplete beta function ratio I_x(a,b) using its series representation.
Definition STK_Funct_betaRatio.cpp:180

References _T, STK::Funct::betaRatio_ae(), STK::Funct::betaRatio_sr(), and stk_cout.

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

◆ binomial_lpdf_raw() [1/2]

Real STK::Funct::binomial_lpdf_raw ( int  x,
int  n,
Real const &  p 
)
inline

Compute the binomial log-probability mass function.

Compute the function

\[ B(n,p,x) = x\log(p) (n-x)\log(1-p) \log\left(\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)}\right) \]

. The function assume that there is no NA or infinite values and that 0<=p<=1, n>=0, 0<=x<=n

Parameters
x,n,pparameters of the binomial density

Definition at line 292 of file STK_Funct_raw.h.

293{
294 // trivial and limit cases
295 if (p == 0) return((x == 0) ? 0 : -Arithmetic<Real>::infinity());
296 if (p == 1) return((x == n) ? 0 : -Arithmetic<Real>::infinity());
297 if (x == 0) return((n == 0) ? 0 : n*log1p(-p));
298 if (x == n) return(n*std::log(p));
299 // other cases
300 int y = n-x;
301 return( - Const::_LNSQRT2PI_ - 0.5 *( std::log(x) + log1p(-x/(Real)n))
302 + lgammaStirlingError(n)-lgammaStirlingError(x)-lgammaStirlingError(y)
303 - dev0(x, n*p) - dev0(y, n - n*p)
304 );
305}

References STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), and STK::Funct::log1p().

◆ binomial_lpdf_raw() [2/2]

Real STK::Funct::binomial_lpdf_raw ( Real const &  x,
Real const &  n,
Real const &  p 
)
inline

Compute the generalized binomial log-probability mass function.

Compute the function

\[ B(n,p,x) = x\log(p) (n-x)\log(1-p) \log\left(\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)}\right) \]

. The function assume that there is no NA or infinite values and that 0<=p<=1, n>=0, 0<=x<=n

Parameters
x,n,pparameters of the binomial density

Definition at line 266 of file STK_Funct_raw.h.

267{
268 // trivial and limit cases
269 if (p == 0) return((x == 0) ? 0 : -Arithmetic<Real>::infinity());
270 if (p == 1) return((x == n) ? 0 : -Arithmetic<Real>::infinity());
271 if (x == 0) return((n == 0) ? 0 : n*log1p(-p));
272 if (x == n) return(n*std::log(p));
273 // other cases
274 Real y = n-x;
275 return ( - Const::_LNSQRT2PI_ - 0.5 *( std::log(x) + log1p(-x/n))
276 + lgammaStirlingError(n)-lgammaStirlingError(x)-lgammaStirlingError(y)
277 - dev0(x, n*p) - dev0(y, n- n*p)
278 );
279}

References STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), and STK::Funct::log1p().

Referenced by STK::Law::Binomial::lpdf(), and STK::Law::Binomial::lpdf().

◆ binomial_pdf_raw() [1/2]

Real STK::Funct::binomial_pdf_raw ( int  x,
int  n,
Real const &  p 
)
inline

Compute the binomial probability mass function.

Compute the function

\[ B(n,p,x) = p^{x} (1-p)^{n-x} \binom{n}{x} \]

. The function assume that there is no NA or infinite values and that 0<=p<=1, n>=0, 0<=x<=n

Parameters
x,n,pparameters of the binomial density

Definition at line 253 of file STK_Funct_raw.h.

254{ return binomial_pdf_raw(Real(x), Real(n), p);}
STK::Funct::binomial_pdf_raw
Real binomial_pdf_raw(Real const &x, Real const &n, Real const &p)
Compute the generalized binomial probability mass function.
Definition STK_Funct_raw.h:228

References STK::Funct::binomial_pdf_raw().

◆ binomial_pdf_raw() [2/2]

Real STK::Funct::binomial_pdf_raw ( Real const &  x,
Real const &  n,
Real const &  p 
)
inline

Compute the generalized binomial probability mass function.

Compute the function

\[ B(n,p,x) = p^{x} (1-p)^{n-x} \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \]

. The function assume that there is no NA or infinite values and that 0<=p<=1, n>=0, 0<=x<=n

Parameters
x,n,pparameters of the binomial density

Definition at line 228 of file STK_Funct_raw.h.

229{
230 // trivial and limit cases
231 if (p == 0) return ((x == 0) ? 1 : 0);
232 if (p == 1) return ((x == n) ? 1 : 0);
233 if (x == 0) return ((n == 0) ? 1 : std::exp(n*log1p(-p)) );
234 if (x == n) return std::pow(p,n);
235 // other cases
236 Real y = n-x;
237 return ( std::exp(- Const::_LNSQRT2PI_ - 0.5 *( std::log(x) + log1p(-x/n))
238 + lgammaStirlingError(n)-lgammaStirlingError(x)-lgammaStirlingError(y)
239 - dev0(x, n*p) - dev0(y, n- n*p)
240 )
241 );
242}

References STK::Funct::dev0(), STK::Funct::lgammaStirlingError(), and STK::Funct::log1p().

Referenced by STK::Funct::binomial_pdf_raw(), STK::Law::Binomial::pdf(), and STK::Law::Binomial::pdf().

◆ BrentMethod()

template<class Function >
Real STK::Algo::BrentMethod ( IFunction< Function > const &  f,
Real const &  x0,
Real const &  x1,
Real  tol 
)

perform the brent's algorithm.

Parameters
fthe function
x0the first point of the algorithm
x1the second point of the algorithm
tolthe tolerance of the method
Returns
the new value

Definition at line 148 of file STK_Algo_FindZero.h.

149{
150 Real a = x0, b = x1, fa = f(x0), fb = f(x1);
151 // set low and high values
152 if (std::abs(fa) < std::abs(fb))
153 {
154 std::swap(a, b);
155 std::swap(fa, fb);
156 }
157 // trivial case
158 if (std::abs(fb)<tol) return b;
159 // if root is not bracketed secant method have to be used
160 if (fa * fb >0.) { return SecantMethod(f, a, b, tol);}
161 Real bkm1 = a, fbkm1 = fa;
162 bool mflag = true;
163 Real bkm2 = bkm1; // value of bkm2 is not used at the first iteration as mflag = true
164
165 int iter =0;
166#ifdef STK_ANALYSIS_VERBOSE
167 stk_cout << _T("iter = ") << iter <<_T(". bkm1, a, b = ") << bkm1 << _T(" ") << a << _T(" ") << b << _T("\n");
168 stk_cout << _T("iter = ") << iter <<_T(". fbkm1, fak, fbk = ") << fbkm1 << _T(" ") << fa << _T(" ") << fb << _T("\n");
169#endif
170 // start Brent iterations
171 while( std::abs(b-a)>tol )
172 {
173 iter++;
174 Real s;
175 if ((fa != fbkm1) && (fb != fbkm1)) // inverse quadratic approximation
176 {
177 Real dab = fa - fb, dac = fa -fbkm1, dbc = fb -fbkm1;
178 s = a * fb * fbkm1 / dab / dac
179 + b * fa * fbkm1 / -dab / dbc
180 + bkm1 * fa * fb / dac / dbc;
181 }
182 else // secante method
183 { s = b - fb * (b - a) / (fb - fa);}
184 // check if we shall fall-back to dichotomy
185 Real tmp = (3. * a + b) / 4., diff1 =std::abs(b - bkm1), diff2 = std::abs(bkm1-bkm2);
186 if (!( ((s > tmp) && (s < b)) || ((s < tmp) && (s > b)) ) // s not between b and tmp
187 ||( mflag &&( (std::abs(s-b) >= diff1 / 2.)||(diff1 < tol)))
188 ||(!mflag &&( (std::abs(s-b) >= diff2 / 2.)||(diff2 < tol)))
189 ||(s<=f.xmin())||(s>=f.xmax())
190 )
191 {
192 s = (a + b) / 2.;
193 mflag = true;
194 }
195 else { mflag = false;}
196
197 Real fs = f(s);
198 // save history
199 bkm2 = bkm1; bkm1 = b; fbkm1 = fb;
200 // check for bracketing the root
201 if (fa * fs < 0) { b = s; fb = fs; }
202 else { a = s; fa = fs; }
203 // if |f(a)| < |f(b)| then swap (a, b)
204 if (std::abs(fa) < std::abs(fb))
205 {
206 std::swap(a, b);
207 std::swap(fa, fb);
208 }
209#ifdef STK_ANALYSIS_VERBOSE
210 stk_cout << _T("iter = ") << iter <<_T(". bkm1, ak, b = ") << bkm1 << _T(" ") << ak << _T(" ") << b << _T("\n");
211 stk_cout << _T("iter = ") << iter <<_T(". fbkm1, fa, fb = ") << fbkm1 << _T(" ") << fa << _T(" ") << fb << _T("\n");
212#endif
213 if (std::abs(fb) < tol ) return b;
214 if (iter > MAX_ITER) return Arithmetic<Real>::NA();
215 }
216 return b;
217}
MAX_ITER
#define MAX_ITER
Definition STK_Algo_FindZero.h:42
STK::Algo::SecantMethod
Real SecantMethod(IFunction< Function > const &f, Real const &x0, Real const &x1, Real tol)
apply the secant method for finding the zero of a function.
Definition STK_Algo_FindZero.h:79

References _T, MAX_ITER, STK::Arithmetic< Type >::NA(), STK::Algo::SecantMethod(), and stk_cout.

Referenced by STK::Algo::findZero(), and STK::Algo::SecantMethod().

◆ continuedFraction()

template<class Seriea , class Serieb >
Real STK::continuedFraction ( const ISerie< Seriea > &  a,
const ISerie< Serieb > &  b,
int 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 template parameter Serie.

Parameters
iter_maxthe number of iterations
aDenominator serie
bNumerator serie

Definition at line 223 of file STK_Algo.h.

227{
228 // initialize a_n
229 Real an = a.first(); // a_1
230 // note a_1 =0 => cf = 0
231 if (an == 0.) return 0.;
232 // initialize b_n
233 Real bn = b.first(); // b_1
234 // initialize numerator
235 Real Nold = 0, Ncur = an; // = a_1
236 // initialize denominator
237 Real Dold = 1, Dcur = bn; // = b_1
238 // result
239 Real cf = 0.;
240 // iterations
241 for (int n=1; n<=iter_max; n++)
242 {
243 // compute a_n
244 an = a.next();
245 // check trivial case
246 // (a_n =0) and (D_n = 0) => cf = +\infty
247 if (an == 0.) return Ncur/Dcur;
248 // update b_n
249 bn = b.next();
250 // compute numerator and denominator
251 Real Nnew = bn * Ncur + an * Nold;
252 Real Dnew = bn * Dcur + an * Dold;
253 // update old numerator and denominator
254 Nold = Ncur;
255 Dold = Dcur;
256 // update current numerator and denominator
257 Ncur = Nnew;
258 Dcur = Dnew;
259 // normalize if necessary with 2^32
260 if (std::abs(Dcur) > 4294967296.0)
261 {
262 Ncur /= Dcur;
263 Nold /= Dcur;
264 Dold /= Dcur;
265 Dcur = 1.;
266 }
267 // normalize if necessary with 2^32
268 if (std::abs(Ncur) > 4294967296.0)
269 {
270 Ncur /= 4294967296.0;
271 Nold /= 4294967296.0;
272 Dold /= 4294967296.0;
273 Dcur /= 4294967296.0;
274 }
275 // if D_n not to small check cv
276 if (std::abs(Dcur) != 0)
277 {
278 Real cfnew = Ncur/Dcur;
279 // check cv
280 if (std::abs(cf - cfnew) < std::abs(cfnew)*Arithmetic<Real>::epsilon())
281 return cfnew;
282 else
283 cf = cfnew;
284 }
285 }
286 // avoid warning at compilation
287 return cf;
288}

◆ dev0()

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

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

Definition at line 74 of file STK_Funct_raw.h.

75{
76 // special values
77 if (a == 0.) return b;
78 if (b == a) return 0.;
79 // ratio
80 Real v = (a-b)/(a+b);
81 // small v -> use dl
82 if (std::abs(v)<0.125)
83 {
84 Real sum = (a-b)*v, ej = 2*a*v, term;
85 v = v*v;
86 int n =0;
87 do
88 { sum += (term = ((ej *= v) / ( ((++n)<<1) +1 ) ));}
89 while (std::abs(term) > std::abs(sum) * Arithmetic<Real>::epsilon());
90 return sum;
91 }
92 // else compute directly the result
93 return (a*log(double(a/b))+b-a);
94}

References STK::sum().

Referenced by STK::Funct::b1(), STK::Funct::beta_pdf_raw(), STK::Funct::betaRatio_se(), STK::Funct::betaRatio_sr(), STK::Funct::binomial_lpdf_raw(), STK::Funct::binomial_lpdf_raw(), STK::Funct::binomial_pdf_raw(), STK::Law::Poisson::lpdf(), STK::Law::Poisson::lpdf(), STK::Law::Poisson::pdf(), STK::Law::Poisson::pdf(), STK::Funct::poisson_lpdf_raw(), and STK::Funct::poisson_pdf_raw().

◆ erf_raw()

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

Compute the error function erf(a) Compute the function.

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

Parameters
[in]avalue to evaluate the function

Definition at line 401 of file STK_Funct_raw.h.

402{
403 Real T[5] =
404 { 9.60497373987051638749E0,
405 9.00260197203842689217E1,
406 2.23200534594684319226E3,
407 7.00332514112805075473E3,
408 5.55923013010394962768E4
409 };
410 Real U[5] =
411 {
412 /* 1.00000000000000000000E0,*/
413 3.35617141647503099647E1,
414 5.21357949780152679795E2,
415 4.59432382970980127987E3,
416 2.26290000613890934246E4,
417 4.92673942608635921086E4
418 };
419 if ( std::abs(a)> 1.0) return ( 1.0 - erfc(a) );
420 Real z = a * a;
421 return a * evalPolynomial<4>( z, T)
422 / evalPolynomial1<5>( z, U);
423}

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

◆ erfc_raw()

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

Compute the complementary error function erfc(a) Compute the function.

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

Parameters
[in]avalue to evaluate the function

Definition at line 433 of file STK_Funct_raw.h.

434{
435 Real P[9] =
436 { 2.46196981473530512524E-10,
437 5.64189564831068821977E-1,
438 7.46321056442269912687E0,
439 4.86371970985681366614E1,
440 1.96520832956077098242E2,
441 5.26445194995477358631E2,
442 9.34528527171957607540E2,
443 1.02755188689515710272E3,
444 5.57535335369399327526E2
445 };
446 Real Q[8] =
447 {
448 /* 1.00000000000000000000E0,*/
449 1.32281951154744992508E1,
450 8.67072140885989742329E1,
451 3.54937778887819891062E2,
452 9.75708501743205489753E2,
453 1.82390916687909736289E3,
454 2.24633760818710981792E3,
455 1.65666309194161350182E3,
456 5.57535340817727675546E2
457 };
458 Real R[6] =
459 { 5.64189583547755073984E-1,
460 1.27536670759978104416E0,
461 5.01905042251180477414E0,
462 6.16021097993053585195E0,
463 7.40974269950448939160E0,
464 2.97886665372100240670E0
465 };
466 Real S[6] =
467 {
468 /* 1.00000000000000000000E0,*/
469 2.26052863220117276590E0,
470 9.39603524938001434673E0,
471 1.20489539808096656605E1,
472 1.70814450747565897222E1,
473 9.60896809063285878198E0,
474 3.36907645100081516050E0
475 };
476 Real x = std::abs(a);
477 if ( x < 1.0) return ( 1.0 - erf_raw(a) );
478 Real z = exp(-a * a);
479 Real y = ( x < 8.0) ? (z * evalPolynomial<8>(x, P))/evalPolynomial1<8>(x, Q)
480 : (z * evalPolynomial<5>(x, R))/evalPolynomial1<6>(x, S);
481 return (a < 0) ? 2.0 - y : y;
482}
STK::Funct::erf_raw
Real erf_raw(Real const &a)
Compute the error function erf(a) Compute the function.
Definition STK_Funct_raw.h:401

References STK::Funct::erf_raw().

Referenced by STK::Law::Normal::cdf(), and STK::Law::Normal::cdf().

◆ euler()

Real STK::Const::euler ( )

Compute the Euler constant.

◆ evalPolynomial()

template<int N>
Real STK::Funct::evalPolynomial ( Real  x,
const Real *  P 
)
inline

Polynomial evaluator.

Evaluate the quantity

\[ P(x) = P[0] x^n + P[1] x^(n-1) + \ldots + P[n] \]

Note
Coefficients are stored in reverse order. There is no checks for out of bounds.

Definition at line 73 of file STK_Funct_Util.h.

74{ return P[N] + x * evalPolynomial<N-1>(x,P);}
STK::Funct::evalPolynomial
Real evalPolynomial(Real x, const Real *P)
Polynomial evaluator.
Definition STK_Funct_Util.h:73

References STK::Funct::evalPolynomial().

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

◆ evalPolynomial1()

template<int N>
Real STK::Funct::evalPolynomial1 ( Real  x,
const Real *  P 
)
inline

Polynomial evaluator.

Evaluate the quantity

\[ P(x) = 1. * x^n + P[0] x^(n-1) + \ldots + P[n-1] \]

Note
Coefficients are stored in reverse order. There is no checks for out of bounds.
Template Parameters
Ndegree of the polynomial

Definition at line 90 of file STK_Funct_Util.h.

91{ return P[N-1] + x * evalPolynomial1<N-1>(x,P);}
STK::Funct::evalPolynomial1
Real evalPolynomial1(Real x, const Real *P)
Polynomial evaluator.
Definition STK_Funct_Util.h:90

References STK::Funct::evalPolynomial1().

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

◆ expm1()

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

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

Definition at line 164 of file STK_Funct_raw.h.

165{
166 // small values of x
167 if (std::abs(x) < Const::_LN2_)
168 {
169 // a + 1 != 1 -> use compute Taylor serie else use first order
170 // Taylor expansion S = x + (1/2!) x^2 + (1/3!) x^3 + ...
171 if (std::abs(x) > Arithmetic<Real>::epsilon())
172 {
173 Real term = x, sum =x, n =1.;
174 do
175 { sum += (term *= x/(++n));}
176 while (std::abs(term) > std::abs(sum) * Arithmetic<Real>::epsilon()) ;
177 return sum ;
178 }
179 else return (x / 2. + 1.) * x;
180 }
181 return exp(double(x))-1.;
182}

References STK::sum().

◆ factorial() [1/2]

Real STK::Funct::factorial ( int  n)
inline

This function computes $ n! $ for integer argument.

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

Parameters
ngiven value for the factorial function

Definition at line 186 of file STK_Funct_gamma.h.

187{
188 // Check if z is Available and finite
189 if (!Arithmetic<int>::isFinite(n)) return Arithmetic<Real>::NA();
190 if (n < 0) { STKDOMAIN_ERROR_1ARG(factorial,n,"Negative argument");}
191 return factorial_raw(n);
192}
STKDOMAIN_ERROR_1ARG
#define STKDOMAIN_ERROR_1ARG(Where, Arg, Error)
Definition STK_Macros.h:165
STK::Funct::factorial_raw
Real factorial_raw(int)
Definition STK_Funct_gamma.h:193

References STK::Funct::factorial(), STK::Funct::factorial_raw(), STK::Arithmetic< Type >::NA(), and STKDOMAIN_ERROR_1ARG.

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

◆ factorial() [2/2]

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

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
zgiven value for the factorial function

Definition at line 204 of file STK_Funct_gamma.h.

205{
206 // Check if z is Available and finite
207 if (!Arithmetic<Real>::isFinite(z)) return z;
208 // Negative integers or reals arguments not allowed
209 if (z < 0 ||(z != std::floor(z))) { STKDOMAIN_ERROR_1ARG(factorial,z,"Negative or not discrete argument");}
210 return factorial_raw(z);
211}

References STK::Funct::factorial(), STK::Funct::factorial_raw(), and STKDOMAIN_ERROR_1ARG.

◆ findZero()

template<class Function >
Real STK::Algo::findZero ( IFunction< Function > const &  f,
Real const &  x0,
Real const &  x1,
Real  tol 
)

find the zero of a function.

Check if the initial values are inside the domain of definition of the function and call the Brent's method.

Parameters
fthe function
x0the first starting point of the algorithm
x1the second starting point of the algorithm
tolthe tolerance to apply
Returns
the zero of the function if any, a NA value otherwise

Definition at line 231 of file STK_Algo_FindZero.h.

232{
233 if (x0<f.xmin() || x0>f.xmax() || x1<f.xmin() || x1>f.xmax() || tol<=0 )
234 return Arithmetic<Real>::NA();
235 return BrentMethod(f, x0, x1, tol);
236}
STK::Algo::BrentMethod
Real BrentMethod(IFunction< Function > const &f, Real const &x0, Real const &x1, Real tol)
perform the brent's algorithm.
Definition STK_Algo_FindZero.h:148

References STK::Algo::BrentMethod(), and STK::Arithmetic< Type >::NA().

Referenced by STK::JointGammaModel< Array, WColVector >::computeParameters(), STK::ModelGamma_aj_bj< Data_, WColVector_ >::computeParameters(), STK::JointGammaModel< Array, WColVector >::computeParameters(), STK::ModelGamma_aj_bj< Data_, WColVector_ >::computeParameters(), STK::Gamma_a_bjk< Array >::run(), STK::Gamma_a_bk< Array >::run(), STK::Gamma_aj_bjk< Array >::run(), STK::Gamma_aj_bk< Array >::run(), STK::Gamma_ajk_b< Array >::run(), STK::Gamma_ajk_bj< Array >::run(), STK::Gamma_ajk_bjk< Array >::run(), STK::Gamma_ajk_bk< Array >::run(), STK::Gamma_ak_b< Array >::run(), STK::Gamma_ak_bj< Array >::run(), STK::Gamma_ak_bjk< Array >::run(), and STK::Gamma_ak_bk< Array >::run().

◆ g0()

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

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

Definition at line 119 of file STK_Funct_raw.h.

120{
121 // special values
122 if (x == 0.) return 1.;
123 if (x == 1.) return 0.;
124 // general case
125 Real v = (x-1.)/(x+1.);
126 // if 7/9 < x < 9/7 use Taylor serie of log((1+v)/(1-v))
127 if (v < 0.125)
128 {
129 Real sum = (x-1)*v, ej = 2*x*v, term;
130 v = v*v;
131 int n =0;
132 do
133 { sum += (term = ( (ej *= v) /( ((++n)<<1) +1)));}
134 while (std::abs(term) > std::abs(sum) * Arithmetic<Real>::epsilon());
135 return sum;
136 }
137 // else compute directly
138 return x*Real(log(double(x)))-x + 1.;
139}

References STK::sum().

◆ gamma()

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

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 approximation method.

Parameters
zgiven value for the gamma function

Definition at line 272 of file STK_Funct_gamma.h.

273{
274 // Check if z is Available and finite
275 if (!Arithmetic<Real>::isFinite(z)) return z;
276 // Negative integer argument not allowed
277 if ( z<=0 && z == std::floor(z))
278 { STKDOMAIN_ERROR_1ARG(Funct::gamma,z,"Negative integer or zero argument");}
279 return gamma_raw(z);
280}
STK::Funct::gamma_raw
Real gamma_raw(Real const &)
Definition STK_Funct_gamma.h:282

References STK::Funct::gamma(), STK::Funct::gamma_raw(), and STKDOMAIN_ERROR_1ARG.

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

◆ gammaLanczos()

Real STK::Funct::gammaLanczos ( Real const &  z)
inline

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

Parameters
zgiven value for the gamma function

Definition at line 86 of file STK_Funct_gamma.h.

87{
88 // 2 * sqrt(e/pi) = 1.86038273...
89 return 1.8603827342052657173362492472666631120594218414085774528900013
90 * exp((z-0.5)*(log(z+22.118910)-1.))
91 * lanczosSerie(z);
92}
STK::Funct::lanczosSerie
Real lanczosSerie(Real const &z)
Compute the Lanzcos correction series for the gamma function with n = 21 terms.
Definition STK_Funct_gamma.h:73

References STK::Funct::lanczosSerie().

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

◆ gammaRatio()

Real STK::Funct::gammaRatio ( Real const &  a,
Real const &  x,
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
aparameter of the gamma ratio function
xvalue to evaluate the gamma ratio function
lower_tailtrue if we want the lower tail, false otherwise

Definition at line 379 of file STK_Funct_gammaRatio.cpp.

380{
381 // Check if a and x are available
382 if (isNA(a)||isNA(x)) return a;
383 // Negative parameter not allowed
384 if (a<=0)
385 { STKDOMAIN_ERROR_2ARG(Funct::gammaRatioP,a,x,negative parameter);}
386 return gammaRatio_raw(a,x,lower_tail);
387}
STK::Funct::gammaRatio_raw
Real gammaRatio_raw(Real const &a, Real const &x, bool lower_tail)
Definition STK_Funct_gammaRatio.cpp:389

References STK::Funct::gammaRatio_raw(), STK::Funct::gammaRatioP(), STK::isNA(), and STKDOMAIN_ERROR_2ARG.

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

◆ gammaRatioP()

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

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
aparameter of the gamma ratio function
xvalue to evaluate the gamma ratio function

Definition at line 97 of file STK_Funct_gammaRatio.h.

98{ return gammaRatio(a, x, true);}
STK::Funct::gammaRatio
Real gammaRatio(Real const &a, Real const &x, bool lower_tail)
Compute the incomplete gamma functions ratio.
Definition STK_Funct_gammaRatio.cpp:379

References STK::Funct::gammaRatio().

Referenced by STK::Funct::betaRatio_ae(), STK::Law::Gamma::cdf(), STK::Law::Gamma::cdf(), and STK::Funct::gammaRatio().

◆ gammaRatioQ()

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

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
aparameter of the gamma ratio function
xvalue to evaluate the gamma ratio function

Definition at line 85 of file STK_Funct_gammaRatio.h.

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

References STK::Funct::gammaRatio().

Referenced by STK::Funct::betaRatio_ae(), STK::Law::Poisson::cdf(), and STK::Law::Poisson::cdf().

◆ gammaStirling()

Real STK::Funct::gammaStirling ( Real const &  z)
inline

This function computes the gamma function using the Stirling approximation.

This approximation is valid for large values of z.

Parameters
zgiven value for the gamma function

Definition at line 123 of file STK_Funct_gamma.h.

124{ return Const::_SQRT2PI_ * exp((z-0.5)*(log(z)-1.)+stirlingSerie(z)-0.5);}
STK::Funct::stirlingSerie
double stirlingSerie(Real const &z)
Compute the Stirling's series for the lgamma function.
Definition STK_Funct_gamma.h:98

References STK::Funct::stirlingSerie().

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

◆ lanczosSerie()

Real STK::Funct::lanczosSerie ( Real const &  z)
inline

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

Parameters
zgiven value for the lanzcos Series

Definition at line 73 of file STK_Funct_gamma.h.

74{
75 Real sum = 0.0;
76 for (int k=20; k>=0; k--)
77 sum += Const::lanczosCoefArray[k]/(z+k);
78 return 2.0240434640140357514731512432760e-10 + sum;
79}

References STK::sum().

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

◆ lfactorial() [1/2]

Real STK::Funct::lfactorial ( int  n)
inline

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
ngiven value for the factorial function

Definition at line 225 of file STK_Funct_gamma.h.

226{
227 // Check if z is Available and finite
228 if (!Arithmetic<int>::isFinite(n)) return n;
229 // Negative integers or reals arguments not allowed
230 if (n < 0) { STKDOMAIN_ERROR_1ARG(lfactorial,n,"Negative integer argument");}
231 // if n is a small number return a constant
232 return lfactorial_raw(n);
233}
STK::Funct::lfactorial_raw
Real lfactorial_raw(int)
Definition STK_Funct_gamma.h:234

References STK::Funct::lfactorial(), STK::Funct::lfactorial_raw(), and STKDOMAIN_ERROR_1ARG.

Referenced by STK::Funct::lfactorial(), and STK::Funct::lfactorial().

◆ lfactorial() [2/2]

Real STK::Funct::lfactorial ( Real const &  z)
inline

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

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

Parameters
zgiven value for the factorial function

Definition at line 245 of file STK_Funct_gamma.h.

246{
247 // Check if z is Available and finite
248 if (!Arithmetic<Real>::isFinite(z)) return z;
249 // Negative integers or real arguments not allowed
250 if ((z < 0)||(z != std::floor(z)))
251 { STKDOMAIN_ERROR_1ARG(lfactorial,z,"Negative integer or decimal argument");}
252 // if n is a small number return a constant
253 return lfactorial_raw(z);
254}

References STK::Funct::lfactorial(), STK::Funct::lfactorial_raw(), and STKDOMAIN_ERROR_1ARG.

◆ lgamma()

Real STK::Funct::lgamma ( Real const &  z)
inline

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
zgiven value for the gamma function

Definition at line 330 of file STK_Funct_gamma.h.

331{
332 // Check if z is Available and finite
333 if (!Arithmetic<Real>::isFinite(z)) return z;
334 // Negative integer argument not allowed
335 if ( z<=0 && z == std::floor(z))
336 { STKDOMAIN_ERROR_1ARG(Funct::lgamma,z,"Negative integer or zero argument");}
337 return lgamma_raw(z);
338}
STK::Funct::lgamma_raw
Real lgamma_raw(Real const &)
Definition STK_Funct_gamma.h:339

References STK::Funct::lgamma(), STK::Funct::lgamma_raw(), and STKDOMAIN_ERROR_1ARG.

Referenced by STK::Funct::lfactorial_raw(), STK::Funct::lfactorial_raw(), STK::Funct::lgamma(), STK::Funct::lgammaStirlingError(), and STK::GammaBase< Derived >::qValue().

◆ lgammaStirling()

Real STK::Funct::lgammaStirling ( Real const &  z)
inline

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

This approximation is valid for large values of z.

Parameters
zgiven value for log gamma function

Definition at line 133 of file STK_Funct_gamma.h.

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

References STK::Funct::stirlingSerie().

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

◆ lgammaStirlingError() [1/2]

Real STK::Funct::lgammaStirlingError ( int  n)
inline

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 series $ 1/12n - 1/360n^3 + ... $ For other z < 100, uses lgamma directly (don't use this to write lgamma!)

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

Parameters
ngiven value for the gamma function

Definition at line 175 of file STK_Funct_gamma.h.

176{ return (n<100.0) ? Const::lgammaStirlingErrorArray[n] : stirlingSerie(n);}

References STK::Funct::stirlingSerie().

◆ lgammaStirlingError() [2/2]

Real STK::Funct::lgammaStirlingError ( Real const &  z)
inline

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 lgamma directly (don't use this to write lgamma!)

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

Parameters
zgiven value for the gamma function

Definition at line 149 of file STK_Funct_gamma.h.

150{
151 int n = (int)std::floor(z);
152 // z is a discrete value
153 if (z == n)
154 { return (n<100) ? Const::lgammaStirlingErrorArray[n] : stirlingSerie(z);}
155 // z is a discrete value halves
156 if (z== n + 0.5)
157 { return (n<100) ? Const::lgammaStirlingErrorHalvesArray[n] : stirlingSerie(z);}
158 // other cases
159 return (z > 16) ? stirlingSerie(z) : lgamma(z) - (z-0.5)*log(z) + z - Const::_LNSQRT2PI_;
160}
STK::Funct::lgamma
Real lgamma(Real const &)
This function computes the function .
Definition STK_Funct_gamma.h:330

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

Referenced by STK::Funct::b1(), STK::Funct::beta_pdf_raw(), STK::Funct::betaRatio_sr(), STK::Funct::binomial_lpdf_raw(), STK::Funct::binomial_lpdf_raw(), STK::Funct::binomial_pdf_raw(), STK::Law::Poisson::lpdf(), STK::Law::Poisson::lpdf(), STK::Law::Poisson::pdf(), STK::Law::Poisson::pdf(), STK::Funct::poisson_lpdf_raw(), and STK::Funct::poisson_pdf_raw().

◆ log1p()

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

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

Parameters
xvalue to evaluate the function

Definition at line 145 of file STK_Funct_raw.h.

146{
147 // trivial values
148 if (x == 0.) return 0.;
149 if (x == -1.) return(-Arithmetic<Real>::infinity());
150 // small values
151 if (std::abs(x) < .375)
152 {
153 // create functor and compute the alternate serie
154 Serielog1p f(x);
155 return x * (1. - x * sumAlternateSerie(f));
156 }
157 // large values : compute directly the result
158 return log(1. + x);
159}
STK::sumAlternateSerie
Real sumAlternateSerie(const ISerie< Serie > &f, int const &n=0)
Sum an alternating series using the Chebichev polynomials.
Definition STK_Algo.h:65

References STK::sumAlternateSerie().

Referenced by STK::Funct::b1(), STK::Funct::beta_pdf_raw(), STK::Funct::betaRatio_ae(), STK::Funct::binomial_lpdf_raw(), STK::Funct::binomial_lpdf_raw(), and STK::Funct::binomial_pdf_raw().

◆ normal_cdf_raw()

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

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
STK::MultidimRegression
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Definition STK_MultidimRegression.h:52

and

erfc

.

Parameters
xvalue to evaluate the function

Definition at line 498 of file STK_Funct_raw.h.

499{
500 Real t = x * Const::_1_SQRT2_, z = std::abs(t);
501 if ( z < Const::_1_SQRT2_) return 0.5 + 0.5 * erf(t);
502 Real y = 0.5 * erfc(z);
503 return ( t > 0) ? 1. - y : y;
504}

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

◆ normal_pdf_raw()

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

compute the probability distribution function of the normal density.

compute the probability density function of the normal density

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

Parameters
xvalue to evaluate the function

Definition at line 514 of file STK_Funct_raw.h.

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

◆ pi()

Real STK::Const::pi ( )

compute pi

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

◆ pi_2()

Real STK::Const::pi_2 ( )
inline

Compute $ \pi/2 $ using pi()

Definition at line 108 of file STK_Const_Math.h.

108{ return pi()/2.0;}
STK::Const::pi
Real pi()
compute pi

References STK::Const::pi().

◆ poisson_lpdf_raw() [1/2]

Real STK::Funct::poisson_lpdf_raw ( int  x,
Real const &  lambda 
)
inline

Compute the log-poisson density with integer value.

Compute the function:

\[ p(x, \lambda) = -\lambda + x \log(\lambda) -\log(x!) \]

with good accuracy using the partial deviance. This is the version for x integer.

See also
http://www.herine.net/stat/software/dbinom.html
Parameters
xvalue to evaluate the function
lambdavalue of the parameter

Definition at line 390 of file STK_Funct_raw.h.

391{ return poisson_lpdf_raw(Real(x), lambda);}
STK::Funct::poisson_lpdf_raw
Real poisson_lpdf_raw(Real const &x, Real const &lambda)
Compute the log-poisson density.
Definition STK_Funct_raw.h:364

References STK::Funct::poisson_lpdf_raw().

◆ poisson_lpdf_raw() [2/2]

Real STK::Funct::poisson_lpdf_raw ( Real const &  x,
Real const &  lambda 
)
inline

Compute the log-poisson density.

Compute the function:

\[ p(x, \lambda) = -\lambda + x \log(\lambda) -\log(\Gamma(x+1)) \]

with good accuracy using the partial deviance. This is the version for x Real.

See also
http://www.herine.net/stat/software/dbinom.html
Parameters
xvalue to evaluate the function
lambdavalue of the parameter

Definition at line 364 of file STK_Funct_raw.h.

365{
366 // check trivial values
367 if (x<0.) return( -Arithmetic<Real>::infinity() );
368 // if lambda is 0, we have P(X=0) = 1
369 if (lambda==0.) return( (x==0.) ? 0 : -Arithmetic<Real>::infinity() );
370 // special value
371 if (x==0.) return( -lambda );
372 // stirling approximation and deviance
373 return( -lgammaStirlingError(x)-dev0(x, lambda)-Const::_LNSQRT2PI_-std::log(x)/2.);
374}

References STK::Funct::dev0(), and STK::Funct::lgammaStirlingError().

Referenced by STK::Law::Gamma::lpdf(), STK::Law::Gamma::lpdf(), and STK::Funct::poisson_lpdf_raw().

◆ poisson_pdf_raw() [1/2]

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

Compute the poisson density with integer value.

Compute the function:

\[ p(x, \lambda) = e^{-\lambda} \frac{\lambda^x}{x!} \]

with good accuracy using the partial deviance. This is the version for x integer.

See also
http://www.herine.net/stat/software/dbinom.html
Parameters
xvalue to evaluate the function
lambdavalue of the parameter

Definition at line 347 of file STK_Funct_raw.h.

348{ return poisson_pdf_raw(Real(x), lambda);}

References STK::Funct::poisson_pdf_raw().

◆ poisson_pdf_raw() [2/2]

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

Compute the Poisson density.

Compute the function:

\[ p(x, \lambda) = e^{-\lambda} \frac{\lambda^x}{\Gamma(x+1)} \]

with good accuracy using the partial deviance. This is the version for x Real.

See also
http://www.herine.net/stat/software/dbinom.html
Parameters
xvalue to evaluate the function
lambdavalue of the parameter

Definition at line 321 of file STK_Funct_raw.h.

322{
323 // check trivial values
324 if (x<0.) return( 0. );
325 // if lambda is 0, we have P(X=0) = 1
326 if (lambda==0.) return( (x==0.) ? 1. : 0. );
327 // special value
328 if (x==0.) return( Real(std::exp(-lambda)) );
329 // stirling approximation and deviance
330 return( std::exp(-lgammaStirlingError(x)-dev0(x, lambda))/(Const::_SQRT2PI_*std::sqrt(x)));
331}

References STK::Funct::dev0(), and STK::Funct::lgammaStirlingError().

Referenced by STK::Funct::betaRatio_ae(), STK::Law::Gamma::pdf(), STK::Law::Gamma::pdf(), and STK::Funct::poisson_pdf_raw().

◆ psi_raw()

Real STK::Funct::psi_raw ( Real  x)
inline

Compute the psi function.

Compute the psi function

\[ \psi(x) =\frac{d}{dx}(-\ln(\Gamma(x))) \]

the logarithmic derivative of the gamma function.

For integer x, we use the formula

\[ \psi(n) = -\gamma + \sum_{k=1}^{n-1} \frac{1}{k}. \]

This formula is used for 0 < n <= 20.

If x is negative, it is transformed to a positive argument by the reflection formula $ \psi(1-x) = \psi(x) + \pi \cot(\pi x) $.

For general positive x, the argument is made greater than 10 using the recurrence $ \psi(x+1) = \psi(x) + 1/x $. Then the following asymptotic expansion is applied:

\[ \psi(x) = \log(x) - \frac{1}{2} x - \sum_{k=1}^{\infty} \frac{B_{2k}}{2k x^{2k}} \]

where the $ B_{2k} $ are Bernoulli numbers.

Definition at line 542 of file STK_Funct_raw.h.

543{
544 // special value
545 if (x==1.) return -Const::_EULER_;
546 Real y, nz = 0.0, p = std::floor(x);
547 bool negative = false;
548 if( x < 0.0 ) // use reflection formula
549 {
550 negative = true;
551 // Remove the zeros of tan(pi x) by subtracting the nearest integer from x
552 if( (nz = x - p) != 0.5 )
553 {
554 if( nz > 0.5 ) { nz = x - p - 1.;}
555 nz = Const::_PI_/std::tan(Const::_PI_*nz);
556 }
557 x = 1.0 - x;
558 }
559 // check for positive integer up to 20
560 if( (x <= 20.0) && (x == p) )
561 {
562 y = 0.0;
563 for( int i=int(x)-1; i>1; i-- ) { y += 1.0/(Real)i; }
564 y += -Const::_EULER_ + 1.;
565 }
566 else // not integer or greater than 20
567 {
568 Real w = 0.0;
569 while( x < 10.0 ) // shift to
570 {
571 w += 1.0/x;
572 x += 1.0;
573 }
574 Real z = 1.0/(x * x);
575 y = std::log(x) - (0.5/x)
576 - z * evalPolynomial<6>( z, Const::bernouilliNumbersArrayDivBy2K) - w;
577 }
578 return (negative) ? y - nz : y;
579}

References STK::Const::bernouilliNumbersArrayDivBy2K.

Referenced by STK::JointGammaModel< Array, WColVector >::dloglikelihood::fImpl(), STK::ModelGamma_aj_bj< Data_, WColVector_ >::dloglikelihood::fImpl(), STK::hidden::invPsiMLog::fImpl(), and STK::hidden::invPsi::fImpl().

◆ SecantMethod()

template<class Function >
Real STK::Algo::SecantMethod ( IFunction< Function > const &  f,
Real const &  x0,
Real const &  x1,
Real  tol 
)

apply the secant method for finding the zero of a function.

Compute iteratively

\[ x_{n+1} = x_n - f(x_n)\frac{x_{n}-x_{n-1}}{f(x_{n}-f(x_{n-1})} \]

in order to find the zero of the function f.

Some details of the implementation :

  • If at some step the root is bracketed, the method switch immediately to the Brent's method.
  • If at some step the method diverge, a small line search is performed.
  • if no descent is found, the method return a NA number.
Parameters
fthe function
x0the first point of the algorithm
x1the second point of the algorithm
tolthe tolerance of the method
Returns
the zero of the function if the method converge, a NA value otherwise

Definition at line 79 of file STK_Algo_FindZero.h.

80{
81 Real a = x0, b = x1, fa = f(x0), fb = f(x1);
82 // set low and high values
83 if (std::abs(fa) < std::abs(fb))
84 {
85 std::swap(a, b);
86 std::swap(fa, fb);
87 }
88 // trivial case
89 if (std::abs(fb)<tol) return b;
90 // if root is bracketed use Brent method
91 if (fa * fb <0.) { return BrentMethod(f, a, b, tol);}
92 // start secant iterations
93 while( std::abs(b-a)>tol )
94 {
95 Real s = b - fb * (b - a) / (fb - fa);
96 Real delta = b-s;
97 // if we are out of bound, we try a desperate step
98 if (s<f.xmin())
99 { s = (b + std::max(f.xmin(), b - std::abs(b-a))/8 )/2.; delta = b-s;}
100 if (s>f.xmax())
101 { s = (b + std::min(f.xmax(), b + std::abs(b-a))/8 )/2.; delta = b-s;}
102 // compute f(s)
103 Real fs = f(s);
104 // check if we can switch to Brent method
105 if (fb * fs < 0) { return BrentMethod(f, b, s, tol); }
106 // handle divergence
107 if (std::abs(fs)>std::abs(fb))// divergence
108 {
109 bool dv = true;
110 for (int i=0; i<16; i++)
111 {
112 delta /=2.;
113 s = b - delta;
114 fs = f(s);
115 // check if we can switch to Brent method
116 if (fb * fs < 0) { return BrentMethod(f, b, s, tol); }
117 // ok !
118 if (std::abs(fs)<std::abs(fb)) { dv = false; break ;}
119 }
120 // we were not lucky...
121 if (dv) return Arithmetic<Real>::NA();
122 }
123 // update
124 a =b; fa = fb;
125 b =s; fb = fs;
126 // set low and high values
127 if (std::abs(fa) < std::abs(fb))
128 {
129 std::swap(a, b);
130 std::swap(fa, fb);
131 }
132 // trivial case
133 if (std::abs(fb) <tol) return b;
134 }
135 return b;
136}

References STK::Algo::BrentMethod(), and STK::Arithmetic< Type >::NA().

Referenced by STK::Algo::BrentMethod().

◆ stirlingSerie()

double STK::Funct::stirlingSerie ( Real const &  z)
inline

Compute the Stirling's series for the lgamma function.

Parameters
zgiven value for the stirling Series

Definition at line 98 of file STK_Funct_gamma.h.

99{
100 const Real z2 = z * z;
101 return (z <= 50) ? ( Const::stirlingCoefArray[0]
102 + ( Const::stirlingCoefArray[1]
103 + ( Const::stirlingCoefArray[2]
104 + Const::stirlingCoefArray[3]/z2
105 )/z2
106 )/z2
107 )/z
108 : ( Const::stirlingCoefArray[0]
109 + ( Const::stirlingCoefArray[1]
110 + Const::stirlingCoefArray[2]/z2
111 )/z2
112 )/z
113 ;
114}

Referenced by STK::Funct::gammaStirling(), STK::Funct::lgamma_raw(), STK::Funct::lgammaStirling(), STK::Funct::lgammaStirlingError(), and STK::Funct::lgammaStirlingError().

◆ sumAlternateSerie()

template<class Serie >
Real STK::sumAlternateSerie ( const ISerie< Serie > &  f,
int const &  n = 0 
)

Sum an alternating series using the Chebichev polynomials.

Compute

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

where the coefficients of the series are given or computed using the parameter ISerie.

The number 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
fthe ISerie giving the terms of the serie
nthe number of iterations

Definition at line 65 of file STK_Algo.h.

66{
67 // Compute the rate of convergence
68 Real rate = 3.0 + 2.0 * Const::_SQRT2_;
69 // compute the number of iterations if needed
70 int niter = n;
71 if (niter==0)
72 {
73 niter = (int )( std::log(Arithmetic<Real>::epsilon())
74 / (Const::_LN2_-std::log(rate))
75 ) + 1;
76 }
77 // initialization
78 rate = pow(rate, niter);
79 Real b = -2.0/(rate+1.0/rate), c = -1.0, sum = 0.0;
80 // first iteration (k = 0)
81 sum += (c = b-c) * f.first();
82 b *= (- 2.*niter*niter);
83 // iterations
84 for (int k = 1; k<niter; k++)
85 {
86 sum += (c = b-c) * f.next();
87 Real aux = k+1.;
88 b *= ((k + niter)/(k + aux)) * ((k - niter)/aux) * 2.;
89 }
90 return sum;
91}

References STK::sum().

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

◆ sumSerie()

template<class Serie >
Real STK::sumSerie ( const ISerie< Serie > &  f,
int 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 template parameter Serie.

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

Parameters
fthe functor giving the terms of the serie
iter_maxthe number of iterations

Definition at line 110 of file STK_Algo.h.

111{
112 Real e0[7]; // original serie
113 Real e1[6];
114 Real e2[5];
115 Real e3[4];
116 Real e4[3];
117 Real e5[2];
118 Real delta, sum; // e6
119
120 e0[0] = f.first(); //f[0];
121 e0[1] = e0[0] + f.next();
122 e0[2] = e0[1] + f.next();
123 e0[3] = e0[2] + f.next();
124 e0[4] = e0[3] + f.next();
125 e0[5] = e0[4] + f.next();
126 e0[6] = e0[5] + f.next();
127
128 // first epsilon e_{-1}[n] = 0
129 e1[0] = 1/(e0[1] - e0[0]);
130 e1[1] = 1/(e0[2] - e0[1]);
131 e1[2] = 1/(e0[3] - e0[2]);
132 e1[3] = 1/(e0[4] - e0[3]);
133 e1[4] = 1/(e0[5] - e0[4]);
134 e1[5] = 1/(e0[6] - e0[5]);
135
136 // second epsilon
137 e2[0] = e0[0] + 1/(e1[1] - e1[0]);
138 e2[1] = e0[1] + 1/(e1[2] - e1[1]);
139 e2[2] = e0[2] + 1/(e1[3] - e1[2]);
140 e2[3] = e0[3] + 1/(e1[4] - e1[3]);
141 e2[4] = e0[4] + 1/(e1[5] - e1[4]);
142
143 // third epsilon
144 e3[0] = e1[0] + 1/(e2[1] - e2[0]);
145 e3[1] = e1[1] + 1/(e2[2] - e2[1]);
146 e3[2] = e1[2] + 1/(e2[3] - e2[2]);
147 e3[3] = e1[3] + 1/(e2[4] - e2[3]);
148
149 // fourth epsilon
150 e4[0] = e2[0] + 1/(e3[1] - e3[0]);
151 e4[1] = e2[1] + 1/(e3[2] - e3[1]);
152 e4[2] = e2[2] + 1/(e3[3] - e3[2]);
153
154 // fifth epsilon
155 e5[0] = e3[0] + 1/(e4[1] - e4[0]);
156 e5[1] = e3[1] + 1/(e4[2] - e4[1]);
157
158 // sum e6[0]
159 sum = e4[0] + (delta = 1/(e5[1] - e5[0]));
160
161 // main loop
162 for (int n=1; n<=iter_max; n++)
163 { // roll s0
164 e0[4] = e0[5];
165 e0[5] = e0[6];
166 e0[6] += f.next();
167
168 // roll first epsilon
169 e1[3] = e1[4];
170 e1[4] = e1[5];
171 if ( (delta = (e0[6] - e0[5])) == 0) break;
172 e1[5] = 1./delta;
173
174 // roll second epsilon
175 e2[2] = e2[3];
176 e2[3] = e2[4];
177 if ( (delta = (e1[5] - e1[4])) == 0) break;
178 e2[4] = e0[4] + 1./delta;
179
180 // roll third epsilon
181 e3[1] = e3[2];
182 e3[2] = e3[3];
183 if ( (delta = (e2[4] - e2[3])) == 0) break;
184 e3[3] = e1[3] + 1./delta;
185
186 // roll fourth epsilon
187 e4[0] = e4[1];
188 e4[1] = e4[2];
189 if ( (delta = (e3[3] - e3[2])) == 0) break;
190 e4[2] = e2[2] + 1./delta;
191
192 // roll fifth epsilon
193 e5[0] = e5[1];
194 if ( (delta = (e4[2] - e4[1])) == 0) break;
195 e5[1] = e3[1] + 1./delta;
196
197 // roll sixth epsilon and compute sum
198 if ( (delta = (e5[1] - e5[0])) == 0) break;
199 if (!Arithmetic<Real>::isFinite(delta = 1./delta)) break;
200 // sum
201 Real sum1 = e4[0] + delta;
202 if (sum1 == sum) break;
203 sum = sum1;
204 }
205
206 return sum;
207}

References STK::sum().

Variable Documentation

◆ bernouilliNumbersArray

const Real STK::Const::bernouilliNumbersArray[21]
Initial value:
=
{
1.0,
1.0/ 6.0,
-1.0/ 30.0,
1.0/ 42.0,
-1.0/ 30.0,
5.0/ 66.0,
-691.0/ 2730.0,
7.0/ 6.0,
-3617.0/ 510.0,
43867.0/ 798.0,
-174611.0/ 330.0,
854513.0/ 138.0,
-236364091.0/ 2730.0,
8553103.0/ 6.0,
-23749461059.0/ 870.0,
8615841276005.0/ 14322.0,
-7709321041217.0/ 510.0,
2577687858367.0/ 6.0,
-26315271553053477373.0/1919190.0,
2929993913841559.0/ 6.0,
-261082718496449122051.0/ 13530.0
}

Array of the 40th first Bernouilli numbers Bernouilli(n) n=0, 2, 4, ... ,40.

Definition at line 696 of file STK_Const_Sequences.h.

697{
698 1.0, // 0
699 1.0/ 6.0, // 2
700 -1.0/ 30.0, // 4
701 1.0/ 42.0, // 6
702 -1.0/ 30.0, // 8
703 5.0/ 66.0, // 10
704 -691.0/ 2730.0, // 12
705 7.0/ 6.0, // 14
706 -3617.0/ 510.0, // 16
707 43867.0/ 798.0, // 18
708 -174611.0/ 330.0, // 20
709 854513.0/ 138.0, // 22
710 -236364091.0/ 2730.0, // 24
711 8553103.0/ 6.0, // 26
712 -23749461059.0/ 870.0, // 28
713 8615841276005.0/ 14322.0, // 30
714 -7709321041217.0/ 510.0, // 32
715 2577687858367.0/ 6.0, // 34
716 -26315271553053477373.0/1919190.0, // 36
717 2929993913841559.0/ 6.0, // 38
718 -261082718496449122051.0/ 13530.0 // 40
719};

◆ doubleFactorialArray

const Real STK::Const::doubleFactorialArray[39]

array for the double factorial.

The coefficients (2n)!! have been found on http://oeis.org/.

Definition at line 649 of file STK_Const_Sequences.h.

650{
651 1.0, // 0!!
652 1.0, // 1!!
653 2.0, // 2!!
654 3.0, // 3!!
655 8.0, // 4!!
656 15.0, // 5!!
657 48.0, // 6!!
658 105.0, // 7!!
659 384.0, // 8!!
660 945.0, // 9!!
661 3840.0, // 10!!
662 10395.0, // 11!!
663 46080.0, // 12!!
664 135135.0, // 13!!
665 645120.0, // 14!!
666 2027025.0, // 15!!
667 10321920.0, // 16!!
668 34459425.0, // 17!!
669 185794560.0, // 18!!
670 654729075.0, // 19!!
671 3715891200.0, // 20!!
672 13749310575.0, // 21!!
673 81749606400.0, // 22!!
674 316234143225.0, // 23!!
675 1961990553600.0, // 24!!
676 7905853580625.0, // 25!!
677 51011754393600.0, // 26!!
678 213458046676875.0, // 27!!
679 1428329123020800.0, // 28!!
680 6190283353629375.0, // 29!!
681 42849873690624000.0, // 30!!
682 191898783962510625.0, // 31!!
683 1371195958099968000.0, // 32!!
684 6332659870762850625.0, // 33!!
685 46620662575398912000.0, // 34!!
686 221643095476699771875.0, // 35!!
687 1678343852714360832000.0, // 36!!
688 8200794532637891559375.0, // 37!!
689 63777066403145711616000.0 // 38!!
690};

◆ doubleFactorialEvenArray

const Real STK::Const::doubleFactorialEvenArray[20]
Initial value:
=
{
1.0,
2.0,
8.0,
48.0,
384.0,
3840.0,
46080.0,
645120.0,
10321920.0,
185794560.0,
3715891200.0,
81749606400.0,
1961990553600.0,
51011754393600.0,
1428329123020800.0,
42849873690624000.0,
1371195958099968000.0,
46620662575398912000.0,
1678343852714360832000.0,
63777066403145711616000.0
}

array for the double factorial of even numbers.

The coefficients $ (2n)!! = 2\times 4\times 6\times...\times(2n)$. have been found on http://oeis.org/.

Definition at line 620 of file STK_Const_Sequences.h.

621{
622 1.0, // 0!!
623 2.0, // 2!!
624 8.0, // 4!!
625 48.0, // 6!!
626 384.0, // 8!!
627 3840.0, // 10!!
628 46080.0, // 12!!
629 645120.0, // 14!!
630 10321920.0, // 16!!
631 185794560.0, // 18!!
632 3715891200.0, // 20!!
633 81749606400.0, // 22!!
634 1961990553600.0, // 24!!
635 51011754393600.0, // 26!!
636 1428329123020800.0, // 28!!
637 42849873690624000.0, // 30!!
638 1371195958099968000.0, // 32!!
639 46620662575398912000.0, // 34!!
640 1678343852714360832000.0, // 36!!
641 63777066403145711616000.0 // 38!!
642};

◆ doubleFactorialOddArray

const Real STK::Const::doubleFactorialOddArray[19]
Initial value:
=
{
1.0,
3.0,
15.0,
105.0,
945.0,
10395.0,
135135.0,
2027025.0,
34459425.0,
654729075.0,
13749310575.0,
316234143225.0,
7905853580625.0,
213458046676875.0,
6190283353629375.0,
191898783962510625.0,
6332659870762850625.0,
221643095476699771875.0,
8200794532637891559375.0
}

array for the double factorial of odd numbers.

The coefficients $ (2n+1)!! = 1\times 3\times 5\times...\times(2n+1)$. have been found on http://oeis.org/.

Definition at line 591 of file STK_Const_Sequences.h.

592{
593 1.0, // 1!!
594 3.0, // 3!!
595 15.0, // 5!!
596 105.0, // 7!!
597 945.0, // 9!!
598 10395.0, // 11!!
599 135135.0, // 13!!
600 2027025.0, // 15!!
601 34459425.0, // 17!!
602 654729075.0, // 19!!
603 13749310575.0, // 21!!
604 316234143225.0, // 23!!
605 7905853580625.0, // 25!!
606 213458046676875.0, // 27!!
607 6190283353629375.0, // 29!!
608 191898783962510625.0, // 31!!
609 6332659870762850625.0, // 33!!
610 221643095476699771875.0, // 35!!
611 8200794532637891559375.0 // 37!!
612};

◆ factorialArray

const Real STK::Const::factorialArray[51]

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 55 of file STK_Const_Sequences.h.

56{
57 1.0, // 0!
58 1.0, // 1!
59 2.0, // 2!
60 6.0, // 3!
61 24.0, // 4!
62 120.0, // 5!
63 720.0, // 6!
64 5040.0, // 7!
65 40320.0, // 8!
66 362880.0, // 9!
67 3628800.0, // 10!
68 39916800.0, // 11!
69 479001600.0, // 12!
70 6227020800.0, // 13!
71 87178291200.0, // 14!
72 1307674368000.0, // 15!
73 20922789888000.0, // 16!
74 355687428096000.0, // 17!
75 6402373705728000.0, // 18!
76 121645100408832000.0, // 19!
77 2432902008176640000.0, // 20!
78 51090942171709440000.0, // 21!
79 1124000727777607680000.0, // 22!
80 25852016738884976640000.0, // 23!
81 620448401733239439360000.0, // 24!
82 15511210043330985984000000.0, // 25!
83 403291461126605635584000000.0, // 26!
84 10888869450418352160768000000.0, // 27!
85 304888344611713860501504000000.0, // 28!
86 8841761993739701954543616000000.0, // 29!
87 265252859812191058636308480000000.0, // 30!
88 8222838654177922817725562880000000.0, // 31!
89 263130836933693530167218012160000000.0, // 32!
90 8683317618811886495518194401280000000.0, // 33!
91 295232799039604140847618609643520000000.0, // 34!
92 10333147966386144929666651337523200000000.0, // 35!
93 371993326789901217467999448150835200000000.0, // 36!
94 13763753091226345046315979581580902400000000.0, // 37!
95 523022617466601111760007224100074291200000000.0, // 38!
96 20397882081197443358640281739902897356800000000.0, // 39!
97 815915283247897734345611269596115894272000000000.0, // 40!
98 33452526613163807108170062053440751665152000000000.0, // 41!
99 1405006117752879898543142606244511569936384000000000.0, // 42!
100 60415263063373835637355132068513997507264512000000000.0, // 43!
101 2658271574788448768043625811014615890319638528000000000.0, // 44!
102 119622220865480194561963161495657715064383733760000000000.0, // 45!
103 5502622159812088949850305428800254892961651752960000000000.0, // 46!
104 258623241511168180642964355153611979969197632389120000000000.0, // 47!
105 12413915592536072670862289047373375038521486354677760000000000.0, // 48!
106 608281864034267560872252163321295376887552831379210240000000000.0, // 49!
10730414093201713378043612608166064768844377641568960512000000000000.0 // 50!
108};

Referenced by STK::Funct::factorial_raw(), STK::Funct::factorial_raw(), and STK::Funct::gamma_raw().