umontreal.iro.lecuyer.randvar.GammaAcceptanceRejectionGen Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: GammaAcceptanceRejectionGen
* Description: gamma random variate generators using a method that combines
acceptance-rejection with acceptance-complement
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
package umontreal.iro.lecuyer.randvar;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
/**
* This class implements gamma random variate generators using
* a method that combines acceptance-rejection
* with acceptance-complement.
* It uses acceptance-rejection for α < 1 and acceptance-complement
* for
* α >= 1.
* For each gamma variate, the first uniform required is taken from the
* main stream and all additional uniforms (after the first rejection)
* are obtained from the auxiliary stream.
*
*/
public class GammaAcceptanceRejectionGen extends GammaGen {
private RandomStream auxStream;
private static final double q1 = 0.0416666664;
private static final double q2 = 0.0208333723;
private static final double q3 = 0.0079849875;
private static final double q4 = 0.0015746717;
private static final double q5 = -0.0003349403;
private static final double q6 = 0.0003340332;
private static final double q7 = 0.0006053049;
private static final double q8 = -0.0004701849;
private static final double q9 = 0.0001710320;
private static final double a1 = 0.333333333;
private static final double a2 = -0.249999949;
private static final double a3 = 0.199999867;
private static final double a4 = -0.166677482;
private static final double a5 = 0.142873973;
private static final double a6 = -0.124385581;
private static final double a7 = 0.110368310;
private static final double a8 = -0.112750886;
private static final double a9 = 0.104089866;
private static final double e1 = 1.000000000;
private static final double e2 = 0.499999994;
private static final double e3 = 0.166666848;
private static final double e4 = 0.041664508;
private static final double e5 = 0.008345522;
private static final double e6 = 0.001353826;
private static final double e7 = 0.000247453;
private int gen;
// UNURAN parameters for the distribution
private double beta;
private double gamma;
// Generator parameters
// Acceptance rejection combined with acceptance complement
private static final int gs = 0;
private static final int gd = 1;
private double b;
private double ss;
private double s;
private double d;
private double r;
private double q0;
private double c;
private double si;
/**
* Creates a gamma random variate generator with parameters α =
* alpha and λ = lambda, using main stream s and
* auxiliary stream aux.
* The auxiliary stream is used when a random number of uniforms
* is required for a rejection-type generation method.
*
*/
public GammaAcceptanceRejectionGen (RandomStream s, RandomStream aux,
double alpha, double lambda) {
super (s, null);
auxStream = aux;
setParams (alpha, lambda);
beta = 1.0/lambda;
gamma = 0.0;
init ();
}
/**
* Creates a gamma random variate generator with parameters α =
* alpha and λ = lambda, using stream s.
*
*/
public GammaAcceptanceRejectionGen (RandomStream s,
double alpha, double lambda) {
this (s, s, alpha, lambda);
}
/**
* Creates a new generator object for the gamma
* distribution dist, using main stream s and
* auxiliary stream aux.
* The auxiliary stream is used when a random number of uniforms
* is required for a rejection-type generation method.
*
*/
public GammaAcceptanceRejectionGen (RandomStream s, RandomStream aux,
GammaDist dist) {
super (s, dist);
auxStream = aux;
if (dist != null)
setParams (dist.getAlpha(), dist.getLambda());
beta = 1.0/dist.getLambda();
gamma = 0.0;
init ();
}
/**
* Creates a new generator object for the gamma
* distribution dist and stream s for both the main and
* auxiliary stream.
*
*/
public GammaAcceptanceRejectionGen (RandomStream s, GammaDist dist) {
this (s, s, dist);
}
/**
* Returns the auxiliary stream associated with this object.
*
*/
public RandomStream getAuxStream() {
return auxStream;
}
/**
* Generates a new gamma variate with parameters
* α = alpha and λ = lambda, using
* main stream s and auxiliary stream aux.
*
*/
public static double nextDouble (RandomStream s, RandomStream aux,
double alpha, double lambda) {
int gen = gs;
double s_ = 0, ss = 0, d = 0, q0 = 0, r = 0, c = 0, si = 0, b = 0;
// Code taken from UNURAN
if (alpha < 1.0) {
gen = gs;
b = 1.0 + 0.36788794412*alpha; // Step 1
}
else {
gen = gd;
// Step 1. Preparations
ss = alpha - 0.5;
s_ = Math.sqrt (ss);
d = 5.656854249 - 12.0*s_;
// Step 4. Set-up for hat case
r = 1.0 / alpha;
q0 = ((((((((q9*r + q8)*r + q7)*r + q6)*r + q5)*r + q4)*
r + q3)*r + q2)*r + q1)*r;
if (alpha > 3.686) {
if (alpha > 13.022) {
b = 1.77;
si = 0.75;
c = 0.1515/s_;
}
else {
b = 1.654 + 0.0076 * ss;
si = 1.68/s_ + 0.275;
c = 0.062/s_ + 0.024;
}
}
else {
b = 0.463 + s_ - 0.178*ss;
si = 1.235;
c = 0.195/s_ - 0.079 + 0.016*s_;
}
}
return acceptanceRejection
(s, aux, alpha, 1.0/lambda, 0, gen, b, s_, ss, d, r, q0, c, si);
}
public double nextDouble() {
return acceptanceRejection
(stream, auxStream, alpha, beta, gamma, gen, b, s, ss, d, r, q0, c, si);
}
/**
* Same as nextDouble (s, s, alpha, lambda).
*
*/
public static double nextDouble (RandomStream s, double alpha,
double lambda) {
return nextDouble (s, s, alpha, lambda);
}
private static double acceptanceRejection
(RandomStream stream, RandomStream auxStream,
double alpha, double beta, double gamma, int gen,
double b, double s, double ss,
double d, double r, double q0, double c, double si) {
// Code taken from UNURAN
double X, p, U, E;
double q,sign_U,t,v,w,x;
switch (gen) {
case gs:
while (true) {
p = b*stream.nextDouble();
stream = auxStream;
if (p <= 1.0) { // Step 2. Case gds <= 1
X = Math.exp (Math.log (p)/alpha);
if (Math.log (stream.nextDouble()) <= -X)
break;
}
else { // Step 3. Case gds > 1
X = -Math.log ((b - p) / alpha);
if ( Math.log (stream.nextDouble()) <= ((alpha - 1.0)*Math.log (X)))
break;
}
}
break;
case gd:
do {
// Step 2. Normal deviate
t = NormalDist.inverseF01 (stream.nextDouble());
stream = auxStream;
x = s + 0.5*t;
X = x*x;
if (t >= 0.)
break; // Immediate acceptance
// Step 3. Uniform random number
U = stream.nextDouble();
if (d*U <= t*t*t)
break; // Squeeze acceptance
// Step 5. Calculation of q
if (x > 0.) {
// Step 6.
v = t/(s + s);
if (Math.abs (v) > 0.25)
q = q0 - s*t + 0.25*t*t + (ss + ss)*Math.log (1. + v);
else
q = q0 + 0.5*t*t*((((((((a9*v + a8)*v + a7)*v + a6)*
v + a5)*v + a4)*v + a3)*v + a2)*v + a1)*v;
// Step 7. Quotient acceptance
if (Math.log (1. - U) <= q)
break;
}
// Step 8. Double exponential deviate t
while (true) {
do {
E = -Math.log (stream.nextDouble());
U = stream.nextDouble();
U = U + U - 1.;
sign_U = (U > 0) ? 1. : -1.;
t = b + (E*si)*sign_U;
} while (t <= -0.71874483771719); // Step 9. Rejection of t
// Step 10. New q(t)
v = t/(s + s);
if (Math.abs (v) > 0.25)
q = q0 - s*t + 0.25*t*t + (ss + ss)*Math.log (1. + v);
else
q = q0 + 0.5*t*t * ((((((((a9*v + a8)*v + a7)*v + a6)*
v + a5)*v + a4)*v + a3)*v + a2)*v + a1)*v;
// Step 11.
if (q <= 0.)
continue;
if (q > 0.5)
w = Math.exp (q) - 1.;
else
w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
q + e1) * q;
// Step 12. Hat acceptance
if ( c * U * sign_U <= w*Math.exp (E - 0.5*t*t)) {
x = s + 0.5 * t;
X = x * x;
break;
}
}
} while (false);
break;
default: throw new IllegalStateException();
}
return gamma + beta*X;
}
private void init() {
// Code taken from UNURAN
if (alpha < 1.0) {
gen = gs;
b = 1.0 + 0.36788794412*alpha; // Step 1
}
else {
gen = gd;
// Step 1. Preparations
ss = alpha - 0.5;
s = Math.sqrt (ss);
d = 5.656854249 - 12.0*s;
// Step 4. Set-up for hat case
r = 1.0 / alpha;
q0 = ((((((((q9*r + q8)*r + q7)*r + q6)*r + q5)*r + q4)*
r + q3)*r + q2)*r + q1)*r;
if (alpha > 3.686) {
if (alpha > 13.022) {
b = 1.77;
si = 0.75;
c = 0.1515/s;
}
else {
b = 1.654 + 0.0076 * ss;
si = 1.68/s + 0.275;
c = 0.062/s + 0.024;
}
}
else {
b = 0.463 + s - 0.178*ss;
si = 1.235;
c = 0.195/s - 0.079 + 0.016*s;
}
}
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy