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;
}
}
}
}
Ads help maintain this website.