cern.jet.random.tdouble.Distributions Maven / Gradle / Ivy
Show all versions of parallelcolt Show documentation
/*
Copyright (C) 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 cern.jet.random.tdouble;
import cern.jet.random.tdouble.engine.DoubleRandomEngine;
/**
* Contains methods for conveniently generating pseudo-random numbers from
* special distributions such as the Burr, Cauchy, Erlang, Geometric, Lambda,
* Laplace, Logistic, Weibull, etc.
*
* About this class:
*
All distributions are obtained by using a uniform pseudo-random
* number generator. followed by a transformation to the desired distribution.
*
* Example usage:
*
*
* cern.jet.random.engine.RandomEngine generator;
* generator = new cern.jet.random.engine.MersenneTwister(new java.util.Date());
* //generator = new edu.cornell.lassp.houle.RngPack.Ranecu(new java.util.Date());
* //generator = new edu.cornell.lassp.houle.RngPack.Ranmar(new java.util.Date());
* //generator = new edu.cornell.lassp.houle.RngPack.Ranlux(new java.util.Date());
* //generator = AbstractDistribution.makeDefaultGenerator();
* for (int i=1000000; --i >=0; ) {
* int cauchy = Distributions.nextCauchy(generator);
* ...
* }
*
*
* @see cern.jet.random.tdouble.engine.DoubleMersenneTwister
* @see java.util.Random
* @see java.lang.Math
* @author [email protected]
* @version 1.0, 09/24/99
*/
public class Distributions {
/**
* Makes this class non instantiable, but still let's others inherit from
* it.
*/
protected Distributions() {
throw new RuntimeException("Non instantiable");
}
/**
* Returns the probability distribution function of the discrete geometric
* distribution.
*
* p(k) = p * (1-p)^k for k >= 0.
*
*
* @param k
* the argument to the probability distribution function.
* @param p
* the parameter of the probability distribution function.
*/
public static double geometricPdf(int k, double p) {
if (k < 0)
throw new IllegalArgumentException();
return p * Math.pow(1 - p, k);
}
/**
* Returns a random number from the Burr II, VII, VIII, X Distributions.
*
* Implementation: Inversion method. This is a port of
* burr1.c from the C-RAND /
* WIN-RAND library. C-RAND's implementation, in turn, is based upon
*
* L. Devroye (1986): Non-Uniform Random Variate Generation, Springer
* Verlag, New York.
*
*
* @param r
* must be > 0.
* @param nr
* the number of the burr distribution (e.g. 2,7,8,10).
*/
public static double nextBurr1(double r, int nr, DoubleRandomEngine randomGenerator) {
/***********************************************************************
* * Burr II, VII, VIII, X Distributions - Inversion * *
* ***************************************************************** *
* FUNCTION : - burr1 samples a random number from one of the * Burr II,
* VII, VIII, X distributions with * parameter r > 0 , where the no. of
* the * distribution is indicated by a pointer * variable. * REFERENCE
* : - L. Devroye (1986): Non-Uniform Random Variate * Generation,
* Springer Verlag, New York. * SUBPROGRAM : - drand(seed) ...
* (0,1)-uniform generator with * unsigned long integer *seed. * *
**********************************************************************/
double y;
y = Math.exp(Math.log(randomGenerator.raw()) / r); /* y=u^(1/r) */
switch (nr) {
// BURR II
case 2:
return (-Math.log(1 / y - 1));
// BURR VII
case 7:
return (Math.log(2 * y / (2 - 2 * y)) / 2);
// BURR VIII
case 8:
return (Math.log(Math.tan(y * Math.PI / 2.0)));
// BURR X
case 10:
return (Math.sqrt(-Math.log(1 - y)));
}
return 0;
}
/**
* Returns a random number from the Burr III, IV, V, VI, IX, XII
* distributions.
*
* Implementation: Inversion method. This is a port of
* burr2.c from the C-RAND /
* WIN-RAND library. C-RAND's implementation, in turn, is based upon
*
* L. Devroye (1986): Non-Uniform Random Variate Generation, Springer
* Verlag, New York.
*
*
* @param r
* must be > 0.
* @param k
* must be > 0.
* @param nr
* the number of the burr distribution (e.g. 3,4,5,6,9,12).
*/
public static double nextBurr2(double r, double k, int nr, DoubleRandomEngine randomGenerator) {
/***********************************************************************
* * Burr III, IV, V, VI, IX, XII Distribution - Inversion * *
* ***************************************************************** *
* FUNCTION : - burr2 samples a random number from one of the * Burr
* III, IV, V, VI, IX, XII distributions with * parameters r > 0 and k >
* 0, where the no. of * the distribution is indicated by a pointer *
* variable. * REFERENCE : - L. Devroye (1986): Non-Uniform Random
* Variate * Generation, Springer Verlag, New York. * SUBPROGRAM : -
* drand(seed) ... (0,1)-Uniform generator with * unsigned long integer
* *seed. * *
**********************************************************************/
double y, u;
u = randomGenerator.raw(); // U(0/1)
y = Math.exp(-Math.log(u) / r) - 1.0; // u^(-1/r) - 1
switch (nr) {
case 3: // BURR III
return (Math.exp(-Math.log(y) / k)); // y^(-1/k)
case 4: // BURR IV
y = Math.exp(k * Math.log(y)) + 1.0; // y^k + 1
y = k / y;
return (y);
case 5: // BURR V
y = Math.atan(-Math.log(y / k)); // arctan[log(y/k)]
return (y);
case 6: // BURR VI
y = -Math.log(y / k) / r;
y = Math.log(y + Math.sqrt(y * y + 1.0));
return (y);
case 9: // BURR IX
y = 1.0 + 2.0 * u / (k * (1.0 - u));
y = Math.exp(Math.log(y) / r) - 1.0; // y^(1/r) -1
return Math.log(y);
case 12: // BURR XII
return Math.exp(Math.log(y) / k); // y^(1/k)
}
return 0;
}
/**
* Returns a cauchy distributed random number from the standard Cauchy
* distribution C(0,1). math definition and
* animated definition.
*
* p(x) = 1/ (mean*pi * (1+(x/mean)^2)).
*
* Implementation: This is a port of cin.c from the C-RAND /
* WIN-RAND library.
*
*
* @return a number in the open unit interval (0.0,1.0)
* (excluding 0.0 and 1.0).
*/
public static double nextCauchy(DoubleRandomEngine randomGenerator) {
return Math.tan(Math.PI * randomGenerator.raw());
}
/**
* Returns an erlang distributed random number with the given variance and
* mean.
*/
public static double nextErlang(double variance, double mean, DoubleRandomEngine randomGenerator) {
int k = (int) ((mean * mean) / variance + 0.5);
k = (k > 0) ? k : 1;
double a = k / mean;
double prod = 1.0;
for (int i = 0; i < k; i++)
prod *= randomGenerator.raw();
return -Math.log(prod) / a;
}
/**
* Returns a discrete geometric distributed random number; Definition.
*
* p(k) = p * (1-p)^k for k >= 0.
*
* Implementation: Inversion method. This is a port of geo.c
* from the C-RAND /
* WIN-RAND library.
*
* @param p
* must satisfy 0 < p < 1.
*
*/
public static int nextGeometric(double p, DoubleRandomEngine randomGenerator) {
/***********************************************************************
* * Geometric Distribution - Inversion * *
* ***************************************************************** *
* On generating random numbers of a discrete distribution by *
* Inversion normally sequential search is necessary, but in the * case
* of the Geometric distribution a direct transformation is * possible
* because of the special parallel to the continuous * Exponential
* distribution Exp(t): * X - Exp(t): G(x)=1-exp(-tx) * Geo(p):
* pk=G(k+1)-G(k)=exp(-tk)*(1-exp(-t)) * p=1-exp(-t) * A random number
* of the Geometric distribution Geo(p) is * obtained by k=(long int)x,
* where x is from Exp(t) with * parameter t=-log(1-p). * *
* ***************************************************************** *
* FUNCTION: - geo samples a random number from the Geometric *
* distribution with parameter 0
* Implementation: Inversion method. This is a port of
* lamin.c from the C-RAND /
* WIN-RAND library. C-RAND's implementation, in turn, is based upon
*
* J.S. Ramberg, B:W. Schmeiser (1974): An approximate method for generating
* asymmetric variables, Communications ACM 17, 78-82.
*
*/
public static double nextLambda(double l3, double l4, DoubleRandomEngine randomGenerator) {
double l_sign;
if ((l3 < 0) || (l4 < 0))
l_sign = -1.0; // sign(l)
else
l_sign = 1.0;
double u = randomGenerator.raw(); // U(0/1)
double x = l_sign * (Math.exp(Math.log(u) * l3) - Math.exp(Math.log(1.0 - u) * l4));
return x;
}
/**
* Returns a Laplace (Double Exponential) distributed random number from the
* standard Laplace distribution L(0,1).
*
* Implementation: Inversion method. This is a port of
* lapin.c from the C-RAND /
* WIN-RAND library.
*
*
* @return a number in the open unit interval (0.0,1.0)
* (excluding 0.0 and 1.0).
*/
public static double nextLaplace(DoubleRandomEngine randomGenerator) {
double u = randomGenerator.raw();
u = u + u - 1.0;
if (u > 0)
return -Math.log(1.0 - u);
else
return Math.log(1.0 + u);
}
/**
* Returns a random number from the standard Logistic distribution Log(0,1).
*
* Implementation: Inversion method. This is a port of
* login.c from the C-RAND /
* WIN-RAND library.
*/
public static double nextLogistic(DoubleRandomEngine randomGenerator) {
double u = randomGenerator.raw();
return (-Math.log(1.0 / u - 1.0));
}
/**
* Returns a power-law distributed random number with the given exponent and
* lower cutoff.
*
* @param alpha
* the exponent
* @param cut
* the lower cutoff
*/
public static double nextPowLaw(double alpha, double cut, DoubleRandomEngine randomGenerator) {
return cut * Math.pow(randomGenerator.raw(), 1.0 / (alpha + 1.0));
}
/**
* Returns a random number from the standard Triangular distribution in
* (-1,1).
*
* Implementation: Inversion method. This is a port of tra.c
* from the C-RAND /
* WIN-RAND library.
*
*/
public static double nextTriangular(DoubleRandomEngine randomGenerator) {
/***********************************************************************
* * Triangular Distribution - Inversion: x = +-(1-sqrt(u)) * *
* ***************************************************************** *
* FUNCTION : - tra samples a random number from the * standard
* Triangular distribution in (-1,1) * SUBPROGRAM : - drand(seed) ...
* (0,1)-Uniform generator with * unsigned long integer *seed. * *
**********************************************************************/
double u;
u = randomGenerator.raw();
if (u <= 0.5)
return (Math.sqrt(2.0 * u) - 1.0); /* -1 <= x <= 0 */
else
return (1.0 - Math.sqrt(2.0 * (1.0 - u))); /* 0 <= x <= 1 */
}
/**
* Returns a weibull distributed random number. Polar method. See
* Simulation, Modelling & Analysis by Law & Kelton, pp259
*/
public static double nextWeibull(double alpha, double beta, DoubleRandomEngine randomGenerator) {
// Polar method.
// See Simulation, Modelling & Analysis by Law & Kelton, pp259
return Math.pow(beta * (-Math.log(1.0 - randomGenerator.raw())), 1.0 / alpha);
}
/**
* Returns a zipfian distributed random number with the given skew.
*
* Algorithm from page 551 of: Devroye, Luc (1986) `Non-uniform random
* variate generation', Springer-Verlag: Berlin. ISBN 3-540-96305-7 (also
* 0-387-96305-7)
*
* @param z
* the skew of the distribution (must be >1.0).
* @return a zipfian distributed number in the closed interval
* [1,Integer.MAX_VALUE].
*/
public static int nextZipfInt(double z, DoubleRandomEngine randomGenerator) {
/*
* Algorithm from page 551 of: Devroye, Luc (1986) `Non-uniform random
* variate generation', Springer-Verlag: Berlin. ISBN 3-540-96305-7
* (also 0-387-96305-7)
*/
final double b = Math.pow(2.0, z - 1.0);
final double constant = -1.0 / (z - 1.0);
int result = 0;
for (;;) {
double u = randomGenerator.raw();
double v = randomGenerator.raw();
result = (int) (Math.floor(Math.pow(u, constant)));
double t = Math.pow(1.0 + 1.0 / result, z - 1.0);
if (v * result * (t - 1.0) / (b - 1.0) <= t / b)
break;
}
return result;
}
}