math.cern.ProbabilityFuncs Maven / Gradle / Ivy
Show all versions of math-base Show documentation
/*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/*
* Copyright ? 1999 CERN - European Organization for Nuclear Research.
* Permission to use, copy, modify, distribute and sell this software and
* its documentation for any purpose is hereby granted without fee, provided
* that the above copyright notice appear in all copies and that both that
* copyright notice and this permission notice appear in supporting
* documentation. CERN makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without expressed or
* implied warranty.
*/
package math.cern;
import static math.MathConsts.*;
/**
* Custom tailored numerical integration of certain probability distributions.
*
* Implementation:
*
* Some code taken and adapted from the Java 2D Graph
* Package 2.4, which in turn is a port from the Cephes 2.2
* Math Library (C). Most Cephes code (missing from the 2D Graph Package)
* directly ported.
*
* @author [email protected]
* @author [email protected]
*/
public final class ProbabilityFuncs {
/*************************************************
* COEFFICIENTS FOR METHOD normalInverse() *
*************************************************/
/* approximation for 0 <= |y - 0.5| <= 3/8 */
private static final double P0[] = { -5.99633501014107895267E1,
9.80010754185999661536E1, -5.66762857469070293439E1,
1.39312609387279679503E1, -1.23916583867381258016E0, };
private static final double Q0[] = { 1.95448858338141759834E0,
4.67627912898881538453E0, 8.63602421390890590575E1,
-2.25462687854119370527E2, 2.00260212380060660359E2,
-8.20372256168333339912E1, 1.59056225126211695515E1,
-1.18331621121330003142E0, };
/*
* Approximation for interval z = sqrt(-2 log y ) between 2 and 8 i.e., y
* between exp(-2) = .135 and exp(-32) = 1.27e-14.
*/
private static final double P1[] = { 4.05544892305962419923E0,
3.15251094599893866154E1, 5.71628192246421288162E1,
4.40805073893200834700E1, 1.46849561928858024014E1,
2.18663306850790267539E0, -1.40256079171354495875E-1,
-3.50424626827848203418E-2, -8.57456785154685413611E-4, };
private static final double Q1[] = { 1.57799883256466749731E1,
4.53907635128879210584E1, 4.13172038254672030440E1,
1.50425385692907503408E1, 2.50464946208309415979E0,
-1.42182922854787788574E-1, -3.80806407691578277194E-2,
-9.33259480895457427372E-4, };
/*
* Approximation for interval z = sqrt(-2 log y) between 8 and 64 i.e., y
* between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
*/
private static final double P2[] = { 3.23774891776946035970E0,
6.91522889068984211695E0, 3.93881025292474443415E0,
1.33303460815807542389E0, 2.01485389549179081538E-1,
1.23716634817820021358E-2, 3.01581553508235416007E-4,
2.65806974686737550832E-6, 6.23974539184983293730E-9, };
private static final double Q2[] = { 6.02427039364742014255E0,
3.67983563856160859403E0, 1.37702099489081330271E0,
2.16236993594496635890E-1, 1.34204006088543189037E-2,
3.28014464682127739104E-4, 2.89247864745380683936E-6,
6.79019408009981274425E-9, };
/*
* Coefficients for erf()
*/
private static final double T_[] = { 9.60497373987051638749E0,
9.00260197203842689217E1, 2.23200534594684319226E3,
7.00332514112805075473E3, 5.55923013010394962768E4 };
private static final double U_[] = {
// 1.00000000000000000000E0,
3.35617141647503099647E1, 5.21357949780152679795E2,
4.59432382970980127987E3, 2.26290000613890934246E4,
4.92673942608635921086E4 };
/*
* Coefficients for erfc()
*/
private static double P_[] = { 2.46196981473530512524E-10,
5.64189564831068821977E-1, 7.46321056442269912687E0,
4.86371970985681366614E1, 1.96520832956077098242E2,
5.26445194995477358631E2, 9.34528527171957607540E2,
1.02755188689515710272E3, 5.57535335369399327526E2 };
private static double Q_[] = {
// 1.0
1.32281951154744992508E1, 8.67072140885989742329E1,
3.54937778887819891062E2, 9.75708501743205489753E2,
1.82390916687909736289E3, 2.24633760818710981792E3,
1.65666309194161350182E3, 5.57535340817727675546E2 };
private static double R_[] = { 5.64189583547755073984E-1,
1.27536670759978104416E0, 5.01905042251180477414E0,
6.16021097993053585195E0, 7.40974269950448939160E0,
2.97886665372100240670E0 };
private static double S_[] = {
// 1.00000000000000000000E0,
2.26052863220117276590E0, 9.39603524938001434673E0,
1.20489539808096656605E1, 1.70814450747565897222E1,
9.60896809063285878198E0, 3.36907645100081516050E0 };
/**
* Returns the area from zero to x under the beta density function.
*
*
* x
* - -
* | (a+b) | | a-1 b-1
* P(x) = ---------- | t (1-t) dt
* - - | |
* | (a) | (b) -
* 0
*
*
* This function is identical to the incomplete beta integral function
* Gamma.incompleteBeta(a, b, x).
*
* The complemented function is
*
* 1 - P(1-x) = Gamma.incompleteBeta( b, a, x );
*
* @param a parameter a
* @param b parameter b
* @param x the value
* @return the value for x
*/
public static double beta(final double a, final double b, final double x) {
return BetaFun.incompleteBeta(a, b, x);
}
/**
* Returns the area under the right hand tail (from x to infinity)
* of the beta density function.
*
* This function is identical to the incomplete beta integral function
* Gamma.incompleteBeta(b, a, x).
*
* @param a parameter a
* @param b parameter b
* @param x the value
* @return the value for x
*/
public static double betaComplemented(final double a, final double b,
final double x) {
return BetaFun.incompleteBeta(b, a, x);
}
/**
* Returns the sum of the terms 0 through k of the
* Binomial probability density.
*
*
* k
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=0
*
*
* The terms are not summed directly; instead the incomplete beta integral
* is employed, according to the formula
*
* y = binomial(k, n, p) = Gamma.incompleteBeta(n - k, k + 1, 1 - p).
*
* All arguments must be positive,
*
* @param k
* end term.
* @param n
* the number of trials.
* @param p
* the probability of success (must be in (0.0,1.0)).
* @return the sum
*/
public static double binomial(final int k, final int n, final double p) {
if ((p < 0.0) || (p > 1.0)) {
throw new IllegalArgumentException("(p < 0.0) || (p > 1.0)");
}
if ((k < 0) || (n < k)) {
throw new IllegalArgumentException("(k < 0) || (n < k)");
}
if (k == n) {
return 1.0;
}
if (k == 0) {
return Math.pow(1.0 - p, n - k);
}
return BetaFun.incompleteBeta(n - k, k + 1, 1.0 - p);
}
/**
* Returns the sum of the terms k+1 through n of the
* Binomial probability density.
*
*
* n
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=k+1
*
*
* The terms are not summed directly; instead the incomplete beta integral
* is employed, according to the formula
*
* y = binomialComplemented(k, n, p) = Gamma.incompleteBeta(k + 1, n - k, p ).
*
* All arguments must be positive,
*
* @param k
* end term.
* @param n
* the number of trials.
* @param p
* the probability of success (must be in (0.0,1.0)).
* @return the sum
*/
public static double binomialComplemented(final int k, final int n,
final double p) {
if ((p < 0.0) || (p > 1.0)) {
throw new IllegalArgumentException("(p < 0.0) || (p > 1.0)");
}
if ((k < 0) || (n < k)) {
throw new IllegalArgumentException("(k < 0) || (n < k)");
}
if (k == n) {
return 0.0;
}
if (k == 0) {
return 1.0 - Math.pow(1.0 - p, n - k);
}
return BetaFun.incompleteBeta(k + 1, n - k, p);
}
/**
* Returns the area under the left hand tail (from 0 to x) of the
* Chi square probability density function with v degrees of
* freedom.
*
*
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
*
*
* where x is the Chi-square variable.
*
* The incomplete gamma integral is used, according to the formula
*
* y = chiSquare(v, x) = incompleteGamma(v / 2.0, x / 2.0).
*
* The arguments must both be positive.
*
* @param v
* degrees of freedom.
* @param x
* integration end point.
* @return the chi square value
*/
public static double chiSquare(final double v, final double x) {
if (x < 0.0 || v < 1.0) {
return 0.0;
}
return GammaFun.incompleteGamma(v / 2.0, x / 2.0);
}
/**
* Returns the area under the right hand tail (from x to infinity)
* of the Chi square probability density function with v degrees of
* freedom.
*
*
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
*
*
* where x is the Chi-square variable.
*
* The incomplete gamma integral is used, according to the formula
*
* y = chiSquareComplemented(v, x) = incompleteGammaComplement(v / 2.0, x / 2.0)
* .
*
* The arguments must both be positive.
*
* @param v
* degrees of freedom.
* @param x the value
* @return the chi square complemented value
*/
public static double chiSquareComplemented(final double v, final double x) {
if (x < 0.0 || v < 1.0) {
return 0.0;
}
return GammaFun.incompleteGammaComplement(v / 2.0, x / 2.0);
}
/**
* Returns the error function of the normal distribution. The integral is
*
*
* x
* -
* 2 | | 2
* erf(x) = -------- | exp( - t ) dt.
* sqrt(pi) | |
* -
* 0
*
*
* @param x
* the argument to the function.
* @return error function value for x
*/
public static double errorFunction(final double x) {
return erf_(x);
}
/**
* Returns the complementary Error function of the normal distribution.
*
*
* 1 - erf(x) =
*
* inf.
* -
* 2 | | 2
* erfc(x) = -------- | exp( - t ) dt
* sqrt(pi) | |
* -
* x
*
*
* @param x
* the argument to the function.
* @return error function complement for x
*/
public static double errorFunctionComplemented(final double x) {
return erfc_(x);
}
/**
* Returns the inverse errorFunction (inverse of erf).
*
* This implementation is described in the paper: Approximating
* the erfinv function by Mike Giles, Oxford-Man Institute of
* Quantitative Finance, which was published in GPU Computing Gems, volume
* 2, 2010. The source code is available
* here.
*
*
* @param x
* the argument to the function.
* @return v such that x = errorFunction(v)
*/
public static double errorFunctionInverse(final double x) {
if (x < -1.0) {
throw new IllegalArgumentException("x < -1.0");
}
if (x > +1.0) {
throw new IllegalArgumentException("x > +1.0");
}
// beware that the logarithm argument must be
// computed as (1.0 - x) * (1.0 + x),
// it must NOT be simplified as 1.0 - x * x as this
// would induce rounding errors near the boundaries +/-1
double w = -Math.log((1.0 - x) * (1.0 + x));
double p;
if (w < 6.25) {
w -= 3.125;
p = -3.6444120640178196996e-21;
p = -1.685059138182016589e-19 + p * w;
p = 1.2858480715256400167e-18 + p * w;
p = 1.115787767802518096e-17 + p * w;
p = -1.333171662854620906e-16 + p * w;
p = 2.0972767875968561637e-17 + p * w;
p = 6.6376381343583238325e-15 + p * w;
p = -4.0545662729752068639e-14 + p * w;
p = -8.1519341976054721522e-14 + p * w;
p = 2.6335093153082322977e-12 + p * w;
p = -1.2975133253453532498e-11 + p * w;
p = -5.4154120542946279317e-11 + p * w;
p = 1.051212273321532285e-09 + p * w;
p = -4.1126339803469836976e-09 + p * w;
p = -2.9070369957882005086e-08 + p * w;
p = 4.2347877827932403518e-07 + p * w;
p = -1.3654692000834678645e-06 + p * w;
p = -1.3882523362786468719e-05 + p * w;
p = 0.0001867342080340571352 + p * w;
p = -0.00074070253416626697512 + p * w;
p = -0.0060336708714301490533 + p * w;
p = 0.24015818242558961693 + p * w;
p = 1.6536545626831027356 + p * w;
} else if (w < 16.0) {
w = Math.sqrt(w) - 3.25;
p = 2.2137376921775787049e-09;
p = 9.0756561938885390979e-08 + p * w;
p = -2.7517406297064545428e-07 + p * w;
p = 1.8239629214389227755e-08 + p * w;
p = 1.5027403968909827627e-06 + p * w;
p = -4.013867526981545969e-06 + p * w;
p = 2.9234449089955446044e-06 + p * w;
p = 1.2475304481671778723e-05 + p * w;
p = -4.7318229009055733981e-05 + p * w;
p = 6.8284851459573175448e-05 + p * w;
p = 2.4031110387097893999e-05 + p * w;
p = -0.0003550375203628474796 + p * w;
p = 0.00095328937973738049703 + p * w;
p = -0.0016882755560235047313 + p * w;
p = 0.0024914420961078508066 + p * w;
p = -0.0037512085075692412107 + p * w;
p = 0.005370914553590063617 + p * w;
p = 1.0052589676941592334 + p * w;
p = 3.0838856104922207635 + p * w;
} else if (!Double.isInfinite(w)) {
w = Math.sqrt(w) - 5.0;
p = -2.7109920616438573243e-11;
p = -2.5556418169965252055e-10 + p * w;
p = 1.5076572693500548083e-09 + p * w;
p = -3.7894654401267369937e-09 + p * w;
p = 7.6157012080783393804e-09 + p * w;
p = -1.4960026627149240478e-08 + p * w;
p = 2.9147953450901080826e-08 + p * w;
p = -6.7711997758452339498e-08 + p * w;
p = 2.2900482228026654717e-07 + p * w;
p = -9.9298272942317002539e-07 + p * w;
p = 4.5260625972231537039e-06 + p * w;
p = -1.9681778105531670567e-05 + p * w;
p = 7.5995277030017761139e-05 + p * w;
p = -0.00021503011930044477347 + p * w;
p = -0.00013871931833623122026 + p * w;
p = 1.0103004648645343977 + p * w;
p = 4.8499064014085844221 + p * w;
} else {
// this branch does not appear in the original code, it
// was added because the previous branch does not handle
// x = +/-1 correctly. In this case, w is positive infinity
// and as the first coefficient (-2.71e-11) is negative.
// Once the first multiplication is done, p becomes negative
// infinity and remains so throughout the polynomial evaluation.
// So the branch above incorrectly returns negative infinity
// instead of the correct positive infinity.
p = Double.POSITIVE_INFINITY;
}
return p * x;
}
/**
* Returns the integral from zero to x of the gamma probability
* density function.
*
*
* x
* alpha -
* beta | | alpha-1 -beta t
* y = --------- | t e dt
* - | |
* | (alpha) -
* 0
*
*
* The incomplete gamma integral is used, according to the relation
*
* y = Gamma.incompleteGamma(alpha, beta * x).
*
* See http://en.wikipedia.org/wiki/Gamma_distribution#
* Probability_density_function
*
* @param alpha
* the shape parameter of the gamma distribution.
* @param beta
* the rate parameter of the gamma distribution.
* @param x
* integration end point.
* @return gamma function probability for x
*/
public static double gamma(double alpha, double beta, double x) {
if (x < 0.0) {
return 0.0;
}
return GammaFun.incompleteGamma(alpha, beta * x);
}
/**
* Returns the integral from x to infinity of the gamma probability
* density function:
*
*
* inf.
* b -
* a | | b-1 -at
* y = ----- | t e dt
* - | |
* | (b) -
* x
*
*
* The incomplete gamma integral is used, according to the relation
*
* y = Gamma.incompleteGammaComplement(b, a*x).
*
* @param a
* the paramater a (alpha) of the gamma distribution.
* @param b
* the paramater b (beta, lambda) of the gamma distribution.
* @param x
* integration end point.
* @return gamma complemented for x
*/
public static double gammaComplemented(final double a, final double b,
final double x) {
if (x < 0.0) {
return 0.0;
}
return GammaFun.incompleteGammaComplement(b, a * x);
}
/**
* Returns the sum of the terms 0 through k of the
* Negative Binomial Distribution.
*
*
* k
* -- ( n+j-1 ) n j
* > ( ) p (1-p)
* -- ( j )
* j=0
*
*
* In a sequence of Bernoulli trials, this is the probability that
* k or fewer failures precede the n-th success.
*
* The terms are not computed individually; instead the incomplete beta
* integral is employed, according to the formula
*
* y = negativeBinomial(k, n, p) = Gamma.incompleteBeta(n, k + 1, p)
* . All arguments must be positive,
*
* @param k
* end term.
* @param n
* the number of trials.
* @param p
* the probability of success (must be in (0.0,1.0)).
* @return the negative binomial
*/
public static double negativeBinomial(final int k, final int n,
final double p) {
if ((p < 0.0) || (p > 1.0)) {
throw new IllegalArgumentException("(p < 0.0) || (p > 1.0)");
}
if (k < 0) {
return 0.0;
}
return BetaFun.incompleteBeta(n, k + 1, p);
}
/**
* Returns the sum of the terms k+1 to infinity of the Negative
* Binomial distribution.
*
*
* inf
* -- ( n+j-1 ) n j
* > ( ) p (1-p)
* -- ( j )
* j=k+1
*
*
* The terms are not computed individually; instead the incomplete beta
* integral is employed, according to the formula
*
* y = negativeBinomialComplemented(k, n, p) = Gamma.incompleteBeta(k + 1,
* n, 1 - p).
*
* All arguments must be positive,
*
* @param k
* end term.
* @param n
* the number of trials.
* @param p
* the probability of success (must be in (0.0,1.0)).
* @return the negative binomial complement
*/
public static double negativeBinomialComplemented(final int k, final int n,
final double p) {
if ((p < 0.0) || (p > 1.0)) {
throw new IllegalArgumentException("(p < 0.0) || (p > 1.0)");
}
if (k < 0) {
return 0.0;
}
return BetaFun.incompleteBeta(k + 1, n, 1.0 - p);
}
/**
* Returns the area under the Normal (Gaussian) probability density
* function, integrated from minus infinity to x (assumes mean is
* zero, variance is one).
*
*
* x
* -
* 1 | | 2
* normal(x) = --------- | exp( - t /2 ) dt
* sqrt(2pi) | |
* -
* -inf.
*
* = ( 1 + erf(z) ) / 2
* = erfc(z) / 2
*
*
* where z = x/sqrt(2). Computation is via the functions
* errorFunction and errorFunctionComplement.
*
* @param a the value
* @return normal probability for the value
*/
public static double normal(final double a) {
final double x = a * SQRT_TWO_HALF;
final double z = Math.abs(x);
if (z < SQRT_TWO_HALF) {
return 0.5 + 0.5 * erf_(x);
} else {
double y = 0.5 * erfc_(z);
if (x > 0) {
y = 1.0 - y;
}
return y;
}
}
/**
* Returns the area under the Normal (Gaussian) probability density
* function, integrated from minus infinity to x.
*
*
* x
* -
* 1 | | 2
* normal(x) = --------- | exp( - (t-mean) / 2v ) dt
* sqrt(2pi*v)| |
* -
* -inf.
*
*
*
* where v = variance. Computation is via the function
* errorFunction.
*
* @param mean
* the mean of the normal distribution.
* @param variance
* the variance of the normal distribution.
* @param x
* the integration limit.
* @return the normal probability for the integration limit
*/
public static double normal(final double mean, final double variance,
final double x) {
if (x > 0) {
return 0.5 + 0.5 * erf_((x - mean) / Math.sqrt(2.0 * variance));
} else {
return 0.5 - 0.5 * erf_(-(x - mean) / Math.sqrt(2.0 * variance));
}
}
/**
* Returns the error function of the normal distribution; formerly named
* erf. The integral is
*
*
* x
* -
* 2 | | 2
* erf(x) = -------- | exp( - t ) dt.
* sqrt(pi) | |
* -
* 0
*
*
* Implementation: For
* 0 <= |x| < 1, erf_(x) = x * P4(x**2)/Q5(x**2); otherwise
* erf_(x) = 1 - erfc_(x).
*
* Code adapted from the Java 2D
* Graph Package 2.4, which in turn is a port from the Cephes
* 2.2 Math Library (C).
*
* @param x
* the argument to the function.
*/
private static double erf_(final double x) {
if (Math.abs(x) > 1.0) {
return 1.0 - erfc_(x);
}
final double z = x * x;
return x * Polynomial.polevl(z, T_, 4) / Polynomial.p1evl(z, U_, 5);
}
/**
* Returns the complementary Error function of the normal distribution;
* formerly named erfc.
*
*
* 1 - erf(x) =
*
* inf.
* -
* 2 | | 2
* erfc(x) = -------- | exp( - t ) dt
* sqrt(pi) | |
* -
* x
*
*
* Implementation: For small x, erfc_(x) = 1 - erf_(x);
* otherwise rational approximations are computed.
*
* Code adapted from the Java 2D
* Graph Package 2.4, which in turn is a port from the Cephes
* 2.2 Math Library (C).
*
* @param x
* the argument to the function.
*/
private static double erfc_(final double x) {
final double t = (x < 0.0) ? -x : x;
if (t < 1.0) {
return 1.0 - erf_(x);
}
double z = -x * x;
if (z < -MAX_LOG) {
if (x < 0.0) {
return 2.0;
} else {
return 0.0;
}
}
z = Math.exp(z);
double p;
double q;
if (t < 8.0) {
p = Polynomial.polevl(t, P_, 8);
q = Polynomial.p1evl(t, Q_, 8);
} else {
p = Polynomial.polevl(t, R_, 5);
q = Polynomial.p1evl(t, S_, 6);
}
double y = (z * p) / q;
if (x < 0.0) {
y = 2.0 - y;
}
if (y == 0.0) {
if (x < 0.0) {
return 2.0;
} else {
return 0.0;
}
}
return y;
}
/**
* Returns the value, x, for which the area under the Normal
* (Gaussian) probability density function (integrated from minus infinity
* to x) is equal to the argument y (assumes mean is zero,
* variance is one).
*
* For small arguments 0 < y < exp(-2), the program computes
* z = sqrt(-2.0 * log(y)); then the approximation is
* x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). There are two
* rational functions P/Q, one for 0 < y < exp(-32) and the other
* for y up to exp(-2). For larger arguments,
* w = y - 0.5, and
* x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
*
* @param y the value
* @return normal inverse for the value
*/
public static double normalInverse(final double y) {
if (y <= 0.0) {
throw new IllegalArgumentException("y <= 0.0");
}
if (y >= 1.0) {
throw new IllegalArgumentException("y >= 1.0");
}
boolean flag = true;
double y_ = y;
if (y_ > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */
y_ = 1.0 - y_;
flag = false;
}
double x;
if (y_ > 0.13533528323661269189) {
y_ = y_ - 0.5;
final double y2 = y_ * y_;
x = y_
+ y_
* (y2 * Polynomial.polevl(y2, P0, 4) / Polynomial.p1evl(y2,
Q0, 8));
x = x * SQRT_TWO_PI;
return x;
}
x = Math.sqrt(-2.0 * Math.log(y_));
final double x0 = x - Math.log(x) / x;
final double z = 1.0 / x;
double x1;
if (x < 8.0) { /* y > exp(-32) = 1.2664165549e-14 */
x1 = z * Polynomial.polevl(z, P1, 8) / Polynomial.p1evl(z, Q1, 8);
} else {
x1 = z * Polynomial.polevl(z, P2, 8) / Polynomial.p1evl(z, Q2, 8);
}
x = x0 - x1;
if (flag) {
x = -x;
}
return x;
}
/**
* Returns the sum of the first k terms of the Poisson
* distribution.
*
*
* k j
* -- -m m
* > e --
* -- j!
* j=0
*
*
* The terms are not summed directly; instead the incomplete gamma integral
* is employed, according to the relation
*
* y = poisson(k, m) = Gamma.incompleteGammaComplement(k + 1, m).
*
* The arguments must both be positive.
*
* @param k
* number of terms.
* @param mean
* the mean of the Poisson distribution.
* @return sum for the Poisson distribution
*/
public static double poisson(final int k, final double mean) {
if (mean < 0) {
throw new IllegalArgumentException("mean < 0");
}
if (k < 0) {
return 0.0;
}
return GammaFun.incompleteGammaComplement(k + 1, mean);
}
/**
* Returns the sum of the terms k+1 to Infinity of the
* Poisson distribution.
*
*
* inf. j
* -- -m m
* > e --
* -- j!
* j=k+1
*
*
* The terms are not summed directly; instead the incomplete gamma integral
* is employed, according to the formula
*
* y = poissonComplemented(k, m) = Gamma.incompleteGamma(k + 1, m)
* . The arguments must both be positive.
*
* @param k
* start term.
* @param mean
* the mean of the Poisson distribution.
* @return sum of the Poisson distribution
*/
public static double poissonComplemented(final int k, final double mean) {
if (mean < 0) {
throw new IllegalArgumentException("mean < 0");
}
if (k < -1) {
return 0.0;
}
return GammaFun.incompleteGamma(k + 1, mean);
}
/**
* Returns the integral from minus infinity to t of the Student-t
* distribution with k > 0 degrees of freedom.
*
*
* t
* -
* | |
* - | 2 -(k+1)/2
* | ( (k+1)/2 ) | ( x )
* ---------------------- | ( 1 + --- ) dx
* - | ( k )
* sqrt( k pi ) | ( k/2 ) |
* | |
* -
* -inf.
*
*
* Relation to incomplete beta integral:
*
* 1 - studentT(k, t) = 0.5 * Gamma.incompleteBeta(k/2, 1/2, z)
* where z = k/(k + t**2).
*
* Since the function is symmetric around t=0, the area under the
* right tail of the density is found by calling the function with
* -t instead of t.
*
* @param k
* degrees of freedom.
* @param t
* integration end point.
* @return the integral for the StudentT distribution
*/
public static double studentT(final double k, final double t) {
if (k <= 0.0) {
throw new IllegalArgumentException("k <= 0");
}
if (t == 0.0) {
return 0.5;
}
double cdf = 0.5 * BetaFun.incompleteBeta(0.5 * k, 0.5, k / (k + t * t));
if (t >= 0.0) {
// fixes bug reported by
// [email protected]
cdf = 1.0 - cdf;
}
return cdf;
}
/**
* Returns the value, t, for which the area under the Student-t
* probability density function (integrated from minus infinity to
* t) is equal to 1 - alpha/2. The value returned
* corresponds to usual Student t-distribution lookup table for
* talpha[size].
*
* The function uses the studentT function to determine the return value
* iteratively.
*
* @param alpha
* probability
* @param size
* size of data set
* @return the StudentT inverse
*/
public static double studentTInverse(final double alpha, final int size) {
final double cumProb = 1 - alpha / 2; // Cumulative probability
double x1 = normalInverse(cumProb);
// Return inverse of normal for large size
if (size > 200) {
return x1;
}
// Find a pair of x1,x2 that brackets zero
double f1 = studentT(size, x1) - cumProb;
double x2 = x1;
double f2 = f1;
do {
if (f1 > 0) {
x2 = x2 / 2;
} else {
x2 = x2 + x1;
}
f2 = studentT(size, x2) - cumProb;
} while (f1 * f2 > 0);
// Find better approximation
// Pegasus-method
do {
// Calculate slope of secant and t value for which it is 0
final double s12 = (f2 - f1) / (x2 - x1);
final double x3 = x2 - f2 / s12;
// Calculate function value at x3
final double f3 = studentT(size, x3) - cumProb;
if (Math.abs(f3) < 1e-8) { // This criteria needs to be very tight!
// We found a perfect value -> return
return x3;
}
if (f3 * f2 < 0) {
x1 = x2;
f1 = f2;
x2 = x3;
f2 = f3;
} else {
final double g = f2 / (f2 + f3);
f1 = g * f1;
x2 = x3;
f2 = f3;
}
} while (Math.abs(x2 - x1) > 0.001);
if (Math.abs(f2) <= Math.abs(f1)) {
return x2;
} else {
return x1;
}
}
private ProbabilityFuncs() {
}
}