Java Code Examples for cern.jet.random.engine.RandomEngine#raw()
The following examples show how to use
cern.jet.random.engine.RandomEngine#raw() .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 6 votes |
|
/**
* Returns a zipfian distributed random number with the given skew.
* <p>
* 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).
* @returns a zipfian distributed number in the closed interval <tt>[1,Integer.MAX_VALUE]</tt>.
*/
public static int nextZipfInt(double z, RandomEngine 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;
}
Example 2
| Source File: Distributions.java From database with GNU General Public License v2.0 | 6 votes |
|
/**
* Returns a random number from the standard Triangular distribution in (-1,1).
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>tra.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* <p>
*/
public static double nextTriangular(RandomEngine 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 */
}
Example 3
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 6 votes |
|
/**
* Returns a random number from the standard Triangular distribution in (-1,1).
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>tra.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* <p>
*/
public static double nextTriangular(RandomEngine 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 */
}
Example 4
| Source File: RandomSampler.java From jAudioGIT with GNU Lesser General Public License v2.1 | 5 votes |
|
/**
* Computes a sorted random set of <tt>count</tt> elements from the interval <tt>[low,low+N-1]</tt>.
* Since we are talking about a random set, no element will occur more than once.
*
* <p>Running time is <tt>O(N)</tt>, on average. Space requirements are zero.
*
* <p>Numbers are filled into the specified array starting at index <tt>fromIndex</tt> to the right.
* The array is returned sorted ascending in the range filled with numbers.
*
* @param n the total number of elements to choose (must be >= 0).
* @param N the interval to choose random numbers from is <tt>[low,low+N-1]</tt>.
* @param count the number of elements to be filled into <tt>values</tt> by this call (must be >= 0 and <=<tt>n</tt>). Normally, you will set <tt>count=n</tt>.
* @param low the interval to choose random numbers from is <tt>[low,low+N-1]</tt>. Hint: If <tt>low==0</tt>, then draws random numbers from the interval <tt>[0,N-1]</tt>.
* @param values the array into which the random numbers are to be filled; must have a length <tt>>= count+fromIndex</tt>.
* @param fromIndex the first index within <tt>values</tt> to be filled with numbers (inclusive).
* @param randomGenerator a random number generator.
*/
protected static void sampleMethodA(long n, long N, int count, long low, long[] values, int fromIndex, RandomEngine randomGenerator) {
double V, quot, Nreal, top;
long S;
long chosen = -1+low;
top = N-n;
Nreal = N;
while (n>=2 && count>0) {
V = randomGenerator.raw();
S = 0;
quot = top/Nreal;
while (quot > V) {
S++;
top--;
Nreal--;
quot = (quot*top)/Nreal;
}
chosen += S+1;
values[fromIndex++]=chosen;
count--;
Nreal--;
n--;
}
if (count>0) {
// special case n==1
S = (long) (Math.round(Nreal) * randomGenerator.raw());
chosen += S+1;
values[fromIndex]=chosen;
}
}
Example 5
| Source File: Distributions.java From database with GNU General Public License v2.0 | 5 votes |
|
/**
* Returns a discrete geometric distributed random number; <A HREF="http://www.statsoft.com/textbook/glosf.html#Geometric Distribution">Definition</A>.
* <p>
* <tt>p(k) = p * (1-p)^k</tt> for <tt> k >= 0</tt>.
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>geo.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* @param p must satisfy <tt>0 < p < 1</tt>.
* <p>
*/
public static int nextGeometric(double p, RandomEngine 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<p<1. *
* SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with *
* unsigned long integer *seed. *
* *
******************************************************************/
double u = randomGenerator.raw();
return (int)(Math.log(u)/Math.log(1.0-p));
}
Example 6
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 5 votes |
|
/**
* Returns a discrete geometric distributed random number; <A HREF="http://www.statsoft.com/textbook/glosf.html#Geometric Distribution">Definition</A>.
* <p>
* <tt>p(k) = p * (1-p)^k</tt> for <tt> k >= 0</tt>.
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>geo.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* @param p must satisfy <tt>0 < p < 1</tt>.
* <p>
*/
public static int nextGeometric(double p, RandomEngine 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<p<1. *
* SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with *
* unsigned long integer *seed. *
* *
******************************************************************/
double u = randomGenerator.raw();
return (int)(Math.log(u)/Math.log(1.0-p));
}
Example 7
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 5 votes |
|
/**
* Returns an erlang distributed random number with the given variance and mean.
*/
public static double nextErlang(double variance, double mean, RandomEngine 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;
}
Example 8
| Source File: RandomSampler.java From database with GNU General Public License v2.0 | 5 votes |
|
/**
* Computes a sorted random set of <tt>count</tt> elements from the interval <tt>[low,low+N-1]</tt>.
* Since we are talking about a random set, no element will occur more than once.
*
* <p>Running time is <tt>O(N)</tt>, on average. Space requirements are zero.
*
* <p>Numbers are filled into the specified array starting at index <tt>fromIndex</tt> to the right.
* The array is returned sorted ascending in the range filled with numbers.
*
* @param n the total number of elements to choose (must be >= 0).
* @param N the interval to choose random numbers from is <tt>[low,low+N-1]</tt>.
* @param count the number of elements to be filled into <tt>values</tt> by this call (must be >= 0 and <=<tt>n</tt>). Normally, you will set <tt>count=n</tt>.
* @param low the interval to choose random numbers from is <tt>[low,low+N-1]</tt>. Hint: If <tt>low==0</tt>, then draws random numbers from the interval <tt>[0,N-1]</tt>.
* @param values the array into which the random numbers are to be filled; must have a length <tt>>= count+fromIndex</tt>.
* @param fromIndex the first index within <tt>values</tt> to be filled with numbers (inclusive).
* @param randomGenerator a random number generator.
*/
protected static void sampleMethodA(long n, long N, int count, long low, long[] values, int fromIndex, RandomEngine randomGenerator) {
double V, quot, Nreal, top;
long S;
long chosen = -1+low;
top = N-n;
Nreal = N;
while (n>=2 && count>0) {
V = randomGenerator.raw();
S = 0;
quot = top/Nreal;
while (quot > V) {
S++;
top--;
Nreal--;
quot = (quot*top)/Nreal;
}
chosen += S+1;
values[fromIndex++]=chosen;
count--;
Nreal--;
n--;
}
if (count>0) {
// special case n==1
S = (long) (Math.round(Nreal) * randomGenerator.raw());
chosen += S+1;
values[fromIndex]=chosen;
}
}
Example 9
| Source File: Beta.java From jAudioGIT with GNU Lesser General Public License v2.1 | 4 votes |
|
/**
*
*/
protected double b00(double a, double b, RandomEngine randomGenerator) {
double U, V, X, Z;
if (a != a_last || b != b_last) {
a_last = a;
b_last = b;
a_ = a - 1.0;
b_ = b - 1.0;
c = (b * b_) / (a * a_); // b(1-b) / a(1-a)
t = (c == 1.0) ? 0.5 : (1.0 - Math.sqrt(c))/(1.0 - c); // t = t_opt
fa = Math.exp(a_ * Math.log(t));
fb = Math.exp(b_ * Math.log(1.0 - t)); // f(t) = fa * fb
p1 = t/a; // 0 < X < t
p2 = (1.0 - t)/b + p1; // t < X < 1
}
for (;;) {
if ((U = randomGenerator.raw() * p2) <= p1) { // X < t
Z = Math.exp(Math.log(U/p1) / a); X = t*Z;
// squeeze accept: L(x) = 1 + (1 - b)x
if ((V = randomGenerator.raw() * fb) <= 1.0 - b_*X) break;
// squeeze reject: U(x) = 1 + ((1 - t)^(b-1) - 1)/t * x
if (V <= 1.0 + (fb - 1.0)*Z) {
// quotient accept: q(x) = (1 - x)^(b-1) / fb
if (Math.log(V) <= b_ * Math.log(1.0 - X)) break;
}
}
else { // X > t
Z = Math.exp(Math.log((U-p1)/(p2-p1)) / b); X = 1.0 - (1.0 - t)*Z;
// squeeze accept: L(x) = 1 + (1 - a)(1 - x)
if ((V = randomGenerator.raw() * fa) <= 1.0 - a_*(1.0 - X)) break;
// squeeze reject: U(x) = 1 + (t^(a-1) - 1)/(1 - t) * (1 - x)
if (V <= 1.0 + (fa - 1.0)*Z) {
// quotient accept: q(x) = x^(a-1) / fa
if (Math.log(V) <= a_ * Math.log(X)) break;
}
}
}
return(X);
}
Example 10
| Source File: PoissonSlow.java From jAudioGIT with GNU Lesser General Public License v2.1 | 4 votes |
|
/**
* Returns a random number from the distribution; bypasses the internal state.
*/
private int nextInt(double theMean) {
/*
* Adapted from "Numerical Recipes in C".
*/
double xm = theMean;
double g = this.cached_g;
if (xm == -1.0 ) return 0; // not defined
if (xm < SWITCH_MEAN ) {
int poisson = -1;
double product = 1;
do {
poisson++;
product *= randomGenerator.raw();
} while ( product >= g );
// bug in CLHEP 1.4.0: was "} while ( product > g );"
return poisson;
}
else if (xm < MEAN_MAX ) {
double t;
double em;
double sq = this.cached_sq;
double alxm = this.cached_alxm;
RandomEngine rand = this.randomGenerator;
do {
double y;
do {
y = Math.tan(Math.PI*rand.raw());
em = sq*y + xm;
} while (em < 0.0);
em = (double) (int)(em); // faster than em = Math.floor(em); (em>=0.0)
t = 0.9*(1.0 + y*y)* Math.exp(em*alxm - logGamma(em + 1.0) - g);
} while (rand.raw() > t);
return (int) em;
}
else { // mean is too large
return (int) xm;
}
}
Example 11
| Source File: Zeta.java From database with GNU General Public License v2.0 | 4 votes |
|
/**
* Returns a zeta distributed random number.
*/
protected long generateZeta(double ro, double pk, RandomEngine randomGenerator) {
/******************************************************************
* *
* Zeta Distribution - Acceptance Rejection *
* *
******************************************************************
* *
* To sample from the Zeta distribution with parameters ro and pk *
* it suffices to sample variates x from the distribution with *
* density function f(x)=B*{[x+0.5]+pk}^(-(1+ro)) ( x > .5 ) *
* and then deliver k=[x+0.5]. *
* 1/B=Sum[(j+pk)^-(ro+1)] (j=1,2,...) converges for ro >= .5 . *
* It is not necessary to compute B, because variates x are *
* generated by acceptance rejection using density function *
* g(x)=ro*(c+0.5)^ro*(c+x)^-(ro+1). *
* * *
* Integer overflow is possible, when ro is small (ro <= .5) and *
* pk large. In this case a new sample is generated. If ro and pk *
* satisfy the inequality ro > .14 + pk*1.85e-8 + .02*ln(pk) *
* the percentage of overflow is less than 1%, so that the *
* result is reliable. *
* NOTE: The comment above is likely to be nomore valid since *
* the C-version operated on 32-bit integers, while this Java *
* version operates on 64-bit integers. However, the following is *
* still valid: * *
* * *
* If either ro > 100 or k > 10000 numerical problems in *
* computing the theoretical moments arise, therefore ro<=100 and *
* k<=10000 are recommended. *
* *
******************************************************************
* *
* FUNCTION: - zeta samples a random number from the *
* Zeta distribution with parameters ro > 0 and *
* pk >= 0. *
* REFERENCE: - J. Dagpunar (1988): Principles of Random *
* Variate Generation, Clarendon Press, Oxford. *
* *
******************************************************************/
double u,v,e,x;
long k;
if (ro != ro_prev || pk != pk_prev) { // Set-up
ro_prev = ro;
pk_prev = pk;
if (ro<pk) {
c = pk-0.5;
d = 0;
}
else {
c = ro-0.5;
d = (1.0+ro)*Math.log((1.0+pk)/(1.0+ro));
}
}
do {
do {
u=randomGenerator.raw();
v=randomGenerator.raw();
x = (c+0.5)*Math.exp(-Math.log(u)/ro) - c;
} while (x<=0.5 || x>=maxlongint);
k = (int) (x+0.5);
e = -Math.log(v);
} while ( e < (1.0+ro)*Math.log((k+pk)/(x+c)) - d );
return k;
}
Example 12
| Source File: Beta.java From jAudioGIT with GNU Lesser General Public License v2.1 | 4 votes |
|
/**
*
*/
protected double b01(double a, double b, RandomEngine randomGenerator) {
double U, V, X, Z;
if (a != a_last || b != b_last) {
a_last = a;
b_last = b;
a_ = a - 1.0;
b_ = b - 1.0;
t = a_/(a - b); // one step Newton * start value t
fb = Math.exp((b_ - 1.0) * Math.log(1.0 - t)); fa = a - (a + b_)*t;
t -= (t - (1.0 - fa) * (1.0 - t)*fb / b) / (1.0 - fa*fb);
fa = Math.exp(a_ * Math.log(t));
fb = Math.exp(b_ * Math.log(1.0 - t)); // f(t) = fa * fb
if (b_ <= 1.0) {
ml = (1.0 - fb) / t; // ml = -m1
mu = b_ * t; // mu = -m2 * t
}
else {
ml = b_;
mu = 1.0 - fb;
}
p1 = t/a; // 0 < X < t
p2 = fb * (1.0 - t)/b + p1; // t < X < 1
}
for (;;) {
if ((U = randomGenerator.raw() * p2) <= p1) { // X < t
Z = Math.exp(Math.log(U/p1) / a); X = t*Z;
// squeeze accept: L(x) = 1 + m1*x, ml = -m1
if ((V = randomGenerator.raw() ) <= 1.0 - ml*X) break;
// squeeze reject: U(x) = 1 + m2*x, mu = -m2 * t
if (V <= 1.0 - mu*Z) {
// quotient accept: q(x) = (1 - x)^(b-1)
if (Math.log(V) <= b_ * Math.log(1.0 - X)) break;
}
}
else { // X > t
Z = Math.exp(Math.log((U-p1)/(p2-p1)) / b); X = 1.0 - (1.0 - t)*Z;
// squeeze accept: L(x) = 1 + (1 - a)(1 - x)
if ((V = randomGenerator.raw() * fa) <= 1.0 - a_*(1.0 - X)) break;
// squeeze reject: U(x) = 1 + (t^(a-1) - 1)/(1 - t) * (1 - x)
if (V <= 1.0 + (fa - 1.0)*Z) {
// quotient accept: q(x) = (x)^(a-1) / fa
if (Math.log(V) <= a_ * Math.log(X)) break;
}
}
}
return(X);
}
Example 13
| Source File: Beta.java From database with GNU General Public License v2.0 | 4 votes |
|
/**
*
*/
protected double b01(double a, double b, RandomEngine randomGenerator) {
double U, V, X, Z;
if (a != a_last || b != b_last) {
a_last = a;
b_last = b;
a_ = a - 1.0;
b_ = b - 1.0;
t = a_/(a - b); // one step Newton * start value t
fb = Math.exp((b_ - 1.0) * Math.log(1.0 - t)); fa = a - (a + b_)*t;
t -= (t - (1.0 - fa) * (1.0 - t)*fb / b) / (1.0 - fa*fb);
fa = Math.exp(a_ * Math.log(t));
fb = Math.exp(b_ * Math.log(1.0 - t)); // f(t) = fa * fb
if (b_ <= 1.0) {
ml = (1.0 - fb) / t; // ml = -m1
mu = b_ * t; // mu = -m2 * t
}
else {
ml = b_;
mu = 1.0 - fb;
}
p1 = t/a; // 0 < X < t
p2 = fb * (1.0 - t)/b + p1; // t < X < 1
}
for (;;) {
if ((U = randomGenerator.raw() * p2) <= p1) { // X < t
Z = Math.exp(Math.log(U/p1) / a); X = t*Z;
// squeeze accept: L(x) = 1 + m1*x, ml = -m1
if ((V = randomGenerator.raw() ) <= 1.0 - ml*X) break;
// squeeze reject: U(x) = 1 + m2*x, mu = -m2 * t
if (V <= 1.0 - mu*Z) {
// quotient accept: q(x) = (1 - x)^(b-1)
if (Math.log(V) <= b_ * Math.log(1.0 - X)) break;
}
}
else { // X > t
Z = Math.exp(Math.log((U-p1)/(p2-p1)) / b); X = 1.0 - (1.0 - t)*Z;
// squeeze accept: L(x) = 1 + (1 - a)(1 - x)
if ((V = randomGenerator.raw() * fa) <= 1.0 - a_*(1.0 - X)) break;
// squeeze reject: U(x) = 1 + (t^(a-1) - 1)/(1 - t) * (1 - x)
if (V <= 1.0 + (fa - 1.0)*Z) {
// quotient accept: q(x) = (x)^(a-1) / fa
if (Math.log(V) <= a_ * Math.log(X)) break;
}
}
}
return(X);
}
Example 14
| Source File: Beta.java From database with GNU General Public License v2.0 | 4 votes |
|
/**
*
*/
protected double b00(double a, double b, RandomEngine randomGenerator) {
double U, V, X, Z;
if (a != a_last || b != b_last) {
a_last = a;
b_last = b;
a_ = a - 1.0;
b_ = b - 1.0;
c = (b * b_) / (a * a_); // b(1-b) / a(1-a)
t = (c == 1.0) ? 0.5 : (1.0 - Math.sqrt(c))/(1.0 - c); // t = t_opt
fa = Math.exp(a_ * Math.log(t));
fb = Math.exp(b_ * Math.log(1.0 - t)); // f(t) = fa * fb
p1 = t/a; // 0 < X < t
p2 = (1.0 - t)/b + p1; // t < X < 1
}
for (;;) {
if ((U = randomGenerator.raw() * p2) <= p1) { // X < t
Z = Math.exp(Math.log(U/p1) / a); X = t*Z;
// squeeze accept: L(x) = 1 + (1 - b)x
if ((V = randomGenerator.raw() * fb) <= 1.0 - b_*X) break;
// squeeze reject: U(x) = 1 + ((1 - t)^(b-1) - 1)/t * x
if (V <= 1.0 + (fb - 1.0)*Z) {
// quotient accept: q(x) = (1 - x)^(b-1) / fb
if (Math.log(V) <= b_ * Math.log(1.0 - X)) break;
}
}
else { // X > t
Z = Math.exp(Math.log((U-p1)/(p2-p1)) / b); X = 1.0 - (1.0 - t)*Z;
// squeeze accept: L(x) = 1 + (1 - a)(1 - x)
if ((V = randomGenerator.raw() * fa) <= 1.0 - a_*(1.0 - X)) break;
// squeeze reject: U(x) = 1 + (t^(a-1) - 1)/(1 - t) * (1 - x)
if (V <= 1.0 + (fa - 1.0)*Z) {
// quotient accept: q(x) = x^(a-1) / fa
if (Math.log(V) <= a_ * Math.log(X)) break;
}
}
}
return(X);
}
Example 15
| Source File: HyperGeometric.java From jAudioGIT with GNU Lesser General Public License v2.1 | 4 votes |
|
/**
* Returns a random number from the distribution.
*/
protected int hmdu(int N, int M, int n, RandomEngine randomGenerator) {
int I, K;
double p, nu, c, d, U;
if (N != N_last || M != M_last || n != n_last) { // set-up */
N_last = N;
M_last = M;
n_last = n;
Mp = (double) (M + 1);
np = (double) (n + 1); N_Mn = N - M - n;
p = Mp / (N + 2.0);
nu = np * p; /* mode, real */
if ((m = (int) nu) == nu && p == 0.5) { /* mode, integer */
mp = m--;
}
else {
mp = m + 1; /* mp = m + 1 */
}
/* mode probability, using the external function flogfak(k) = ln(k!) */
fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
+ Arithmetic.logFactorial(M) - Arithmetic.logFactorial(M - m) - Arithmetic.logFactorial(m)
- Arithmetic.logFactorial(N) + Arithmetic.logFactorial(N - n) + Arithmetic.logFactorial(n) );
/* safety bound - guarantees at least 17 significant decimal digits */
/* b = min(n, (long int)(nu + k*c')) */
b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));
if (b > n) b = n;
}
for (;;) {
if ((U = randomGenerator.raw() - fm) <= 0.0) return(m);
c = d = fm;
/* down- and upward search from the mode */
for (I = 1; I <= m; I++) {
K = mp - I; /* downward search */
c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
if ((U -= c) <= 0.0) return(K - 1);
K = m + I; /* upward search */
d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
if ((U -= d) <= 0.0) return(K);
}
/* upward search from K = 2m + 1 to K = b */
for (K = mp + m; K <= b; K++) {
d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
if ((U -= d) <= 0.0) return(K);
}
}
}
Example 16
| Source File: HyperGeometric.java From database with GNU General Public License v2.0 | 4 votes |
|
/**
* Returns a random number from the distribution.
*/
protected int hmdu(int N, int M, int n, RandomEngine randomGenerator) {
int I, K;
double p, nu, c, d, U;
if (N != N_last || M != M_last || n != n_last) { // set-up */
N_last = N;
M_last = M;
n_last = n;
Mp = (double) (M + 1);
np = (double) (n + 1); N_Mn = N - M - n;
p = Mp / (N + 2.0);
nu = np * p; /* mode, real */
if ((m = (int) nu) == nu && p == 0.5) { /* mode, integer */
mp = m--;
}
else {
mp = m + 1; /* mp = m + 1 */
}
/* mode probability, using the external function flogfak(k) = ln(k!) */
fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
+ Arithmetic.logFactorial(M) - Arithmetic.logFactorial(M - m) - Arithmetic.logFactorial(m)
- Arithmetic.logFactorial(N) + Arithmetic.logFactorial(N - n) + Arithmetic.logFactorial(n) );
/* safety bound - guarantees at least 17 significant decimal digits */
/* b = min(n, (long int)(nu + k*c')) */
b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));
if (b > n) b = n;
}
for (;;) {
if ((U = randomGenerator.raw() - fm) <= 0.0) return(m);
c = d = fm;
/* down- and upward search from the mode */
for (I = 1; I <= m; I++) {
K = mp - I; /* downward search */
c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
if ((U -= c) <= 0.0) return(K - 1);
K = m + I; /* upward search */
d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
if ((U -= d) <= 0.0) return(K);
}
/* upward search from K = 2m + 1 to K = b */
for (K = mp + m; K <= b; K++) {
d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
if ((U -= d) <= 0.0) return(K);
}
}
}
Example 17
| Source File: Distributions.java From database with GNU General Public License v2.0 | 4 votes |
|
/**
* Returns a random number from the Burr III, IV, V, VI, IX, XII distributions.
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>burr2.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* C-RAND's implementation, in turn, is based upon
* <p>
* L. Devroye (1986): Non-Uniform Random Variate Generation, Springer Verlag, New York.
* <p>
* @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, RandomEngine 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;
}
Example 18
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 3 votes |
|
/**
* Returns a Laplace (Double Exponential) distributed random number from the standard Laplace distribution L(0,1).
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>lapin.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
* <p>
* @returns a number in the open unit interval <code>(0.0,1.0)</code> (excluding 0.0 and 1.0).
*/
public static double nextLaplace(RandomEngine 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);
}
Example 19
| Source File: Distributions.java From database with GNU General Public License v2.0 | 2 votes |
|
/**
* Returns a random number from the standard Logistic distribution Log(0,1).
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>login.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
*/
public static double nextLogistic(RandomEngine randomGenerator) {
double u = randomGenerator.raw();
return(-Math.log(1.0 / u-1.0));
}
Example 20
| Source File: Distributions.java From jAudioGIT with GNU Lesser General Public License v2.1 | 2 votes |
|
/**
* Returns a random number from the standard Logistic distribution Log(0,1).
* <p>
* <b>Implementation:</b> Inversion method.
* This is a port of <tt>login.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
*/
public static double nextLogistic(RandomEngine randomGenerator) {
double u = randomGenerator.raw();
return(-Math.log(1.0 / u-1.0));
}
