RNG-30: Discrete distributions. Code copied and adapted from classes in "o.a.c.math4.distribution" package.
Classes "InternalGamma" and "InternalUtils" provide stripped down versions of code present in Commons Math. Ideally, they should be replaced with a dependency upon a truly low-level math component. Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/74bbdd27 Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/74bbdd27 Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/74bbdd27 Branch: refs/heads/master Commit: 74bbdd27cd38c55b9b9980b0cdfcc889cbe85c64 Parents: 76e1700 Author: Gilles <er...@apache.org> Authored: Sat Nov 12 16:39:55 2016 +0100 Committer: Gilles <er...@apache.org> Committed: Sat Nov 12 16:39:55 2016 +0100 ---------------------------------------------------------------------- .../rng/sampling/AbstractBaseSampler.java | 14 + .../distribution/DiscreteUniformSampler.java | 76 ++++ .../sampling/distribution/InternalGamma.java | 394 +++++++++++++++++++ .../sampling/distribution/InternalUtils.java | 145 +++++++ .../InverseMethodParetoSampler.java | 55 +++ .../sampling/distribution/PoissonSampler.java | 164 ++++++++ .../RejectionInversionZipfSampler.java | 244 ++++++++++++ 7 files changed, 1092 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/AbstractBaseSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/AbstractBaseSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/AbstractBaseSampler.java index 00789da..f97687c 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/AbstractBaseSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/AbstractBaseSampler.java @@ -40,6 +40,20 @@ public class AbstractBaseSampler { return rng.nextDouble(); } + /** + * @return a random {@code int} value. + */ + protected int nextInt() { + return rng.nextInt(); + } + + /** + * @return a random {@code int} value in the interval {@code [0, max)}. + */ + protected int nextInt(int max) { + return rng.nextInt(max); + } + /** {@inheritDoc} */ @Override public String toString() { http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DiscreteUniformSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DiscreteUniformSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DiscreteUniformSampler.java new file mode 100644 index 0000000..f0a3e8d --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DiscreteUniformSampler.java @@ -0,0 +1,76 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ + +package org.apache.commons.rng.sampling.distribution; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.AbstractDiscreteSampler; + +/** + * Discrete uniform distribution sampler. + */ +public class DiscreteUniformSampler extends AbstractDiscreteSampler { + /** Lower bound. */ + private final int lower; + /** Upper bound. */ + private final int upper; + + /** + * @param rng Generator of uniformly distributed random numbers. + * @param lower Lower bound (inclusive) of the distribution. + * @param upper Upper bound (inclusive) of the distribution. + * @throws IllegalArgumentException if {@code lower > upper}. + */ + public DiscreteUniformSampler(UniformRandomProvider rng, + int lower, + int upper) { + super(rng); + if (lower > upper) { + throw new IllegalArgumentException(lower + " > " + upper); + } + + this.lower = lower; + this.upper = upper; + } + + /** {@inheritDoc} */ + @Override + public int sample() { + final int max = (upper - lower) + 1; + if (max <= 0) { + // The range is too wide to fit in a positive int (larger + // than 2^31); as it covers more than half the integer range, + // we use a simple rejection method. + while (true) { + final int r = nextInt(); + if (r >= lower && + r <= upper) { + return r; + } + } + } else { + // We can shift the range and directly generate a positive int. + return lower + nextInt(max); + } + } + + /** {@inheritDoc} */ + @Override + public String toString() { + return "[Uniform distribution " + super.toString() + "]"; + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalGamma.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalGamma.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalGamma.java new file mode 100644 index 0000000..2feb273 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalGamma.java @@ -0,0 +1,394 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ +package org.apache.commons.rng.sampling.distribution; + +/** + * <h3> + * Adapted and stripped down copy of class + * {@code "org.apache.commons.math4.special.Gamma"}. + * TODO: Include it in a "core" component upon which high-level functionality + * such as sampling can depend. + * </h3> + * + * <p> + * This is a utility class that provides computation methods related to the + * Γ (Gamma) family of functions. + * </p> + * <p> + * Implementation of {@link #invGamma1pm1(double)} and + * {@link #logGamma1p(double)} is based on the algorithms described in + * <ul> + * <li><a href="http://dx.doi.org/10.1145/22721.23109";>Didonato and Morris + * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and + * their Inverse</em>, TOMS 12(4), 377-393,</li> + * <li><a href="http://dx.doi.org/10.1145/131766.131776";>Didonato and Morris + * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the + * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li> + * </ul> + * and implemented in the + * <a href="http://www.dtic.mil/docs/citations/ADA476840";>NSWC Library of Mathematical Functions</a>, + * available + * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html";>here</a>. + * This library is "approved for public release", and the + * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf";>Copyright guidance</a> + * indicates that unless otherwise stated in the code, all FORTRAN functions in + * this library are license free. Since no such notice appears in the code these + * functions can safely be ported to Commons-Math. + * </p> + * + */ +class InternalGamma { // Class is package-private on purpose; do not make it public. + /** + * Constant \( g = \frac{607}{128} \) in the {@link #lanczos(double) Lanczos approximation}. + */ + public static final double LANCZOS_G = 607.0 / 128.0; + + /** Lanczos coefficients */ + private static final double[] LANCZOS = { + 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5, + }; + + /** Avoid repeated computation of log of 2 PI in logGamma */ + private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI); + + /* + * Constants for the computation of double invGamma1pm1(double). + * Copied from DGAM1 in the NSWC library. + */ + + /** The constant {@code A0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08; + + /** The constant {@code A1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08; + + /** The constant {@code B1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00; + + /** The constant {@code B2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01; + + /** The constant {@code B3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03; + + /** The constant {@code B4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05; + + /** The constant {@code B5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05; + + /** The constant {@code B6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06; + + /** The constant {@code B7} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07; + + /** The constant {@code B8} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09; + + /** The constant {@code P0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08; + + /** The constant {@code P1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08; + + /** The constant {@code P2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09; + + /** The constant {@code P3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10; + + /** The constant {@code P4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11; + + /** The constant {@code P5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12; + + /** The constant {@code P6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14; + + /** The constant {@code Q1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00; + + /** The constant {@code Q2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01; + + /** The constant {@code Q3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02; + + /** The constant {@code Q4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03; + + /** The constant {@code C} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00; + + /** The constant {@code C0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00; + + /** The constant {@code C1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00; + + /** The constant {@code C2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01; + + /** The constant {@code C3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00; + + /** The constant {@code C4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01; + + /** The constant {@code C5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02; + + /** The constant {@code C6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02; + + /** The constant {@code C7} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02; + + /** The constant {@code C8} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03; + + /** The constant {@code C9} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03; + + /** The constant {@code C10} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04; + + /** The constant {@code C11} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05; + + /** The constant {@code C12} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05; + + /** The constant {@code C13} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06; + + /** + * Class contains only static methods. + */ + private InternalGamma() {} + + /** + * Computes the function \( \ln \Gamma(x) \) for \( x > 0 \). + * + * <p> + * For \( x \leq 8 \), the implementation is based on the double precision + * implementation in the <em>NSWC Library of Mathematics Subroutines</em>, + * {@code DGAMLN}. For \( x \geq 8 \), the implementation is based on + * </p> + * + * <ul> + * <li><a href="http://mathworld.wolfram.com/GammaFunction.html";>Gamma + * Function</a>, equation (28).</li> + * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html";> + * Lanczos Approximation</a>, equations (1) through (5).</li> + * <li><a href="http://my.fit.edu/~gabdo/gamma.txt";>Paul Godfrey, A note on + * the computation of the convergent Lanczos complex Gamma + * approximation</a></li> + * </ul> + * + * @param x Argument. + * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}. + */ + public static double logGamma(double x) { + double ret; + + if (Double.isNaN(x) || (x <= 0.0)) { + ret = Double.NaN; + } else if (x < 0.5) { + return logGamma1p(x) - Math.log(x); + } else if (x <= 2.5) { + return logGamma1p((x - 0.5) - 0.5); + } else if (x <= 8.0) { + final int n = (int) Math.floor(x - 1.5); + double prod = 1.0; + for (int i = 1; i <= n; i++) { + prod *= x - i; + } + return logGamma1p(x - (n + 1)) + Math.log(prod); + } else { + double sum = lanczos(x); + double tmp = x + LANCZOS_G + .5; + ret = ((x + .5) * Math.log(tmp)) - tmp + + HALF_LOG_2_PI + Math.log(sum / x); + } + + return ret; + } + + /** + * Computes the Lanczos approximation used to compute the gamma function. + * + * <p> + * The Lanczos approximation is related to the Gamma function by the + * following equation + * \[ + * \Gamma(x) = \sqrt{2\pi} \, \frac{(g + x + \frac{1}{2})^{x + \frac{1}{2}} \, e^{-(g + x + \frac{1}{2})} \, \mathrm{lanczos}(x)} + * {x} + * \] + * where \(g\) is the Lanczos constant. + * </p> + * + * @param x Argument. + * @return The Lanczos approximation. + * + * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html";>Lanczos Approximation</a> + * equations (1) through (5), and Paul Godfrey's + * <a href="http://my.fit.edu/~gabdo/gamma.txt";>Note on the computation + * of the convergent Lanczos complex Gamma approximation</a> + */ + public static double lanczos(final double x) { + double sum = 0.0; + for (int i = LANCZOS.length - 1; i > 0; --i) { + sum += LANCZOS[i] / (x + i); + } + return sum + LANCZOS[0]; + } + + /** + * Computes the function \( \frac{1}{\Gamma(1 + x)} - 1 \) for \( -0.5 \leq x \leq 1.5 \). + * <p> + * This implementation is based on the double precision implementation in + * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGAM1}. + * </p> + * + * @param x Argument. + * @return \( \frac{1}{\Gamma(1 + x)} - 1 \) + * @throws IllegalArgumentException if {@code x < -0.5} + * @throws IllegalArgumentException if {@code x > 1.5} + */ + public static double invGamma1pm1(final double x) { + + if (x < -0.5) { + throw new IllegalArgumentException(); + } + if (x > 1.5) { + throw new IllegalArgumentException(); + } + + final double ret; + final double t = x <= 0.5 ? x : (x - 0.5) - 0.5; + if (t < 0.0) { + final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1; + double b = INV_GAMMA1P_M1_B8; + b = INV_GAMMA1P_M1_B7 + t * b; + b = INV_GAMMA1P_M1_B6 + t * b; + b = INV_GAMMA1P_M1_B5 + t * b; + b = INV_GAMMA1P_M1_B4 + t * b; + b = INV_GAMMA1P_M1_B3 + t * b; + b = INV_GAMMA1P_M1_B2 + t * b; + b = INV_GAMMA1P_M1_B1 + t * b; + b = 1.0 + t * b; + + double c = INV_GAMMA1P_M1_C13 + t * (a / b); + c = INV_GAMMA1P_M1_C12 + t * c; + c = INV_GAMMA1P_M1_C11 + t * c; + c = INV_GAMMA1P_M1_C10 + t * c; + c = INV_GAMMA1P_M1_C9 + t * c; + c = INV_GAMMA1P_M1_C8 + t * c; + c = INV_GAMMA1P_M1_C7 + t * c; + c = INV_GAMMA1P_M1_C6 + t * c; + c = INV_GAMMA1P_M1_C5 + t * c; + c = INV_GAMMA1P_M1_C4 + t * c; + c = INV_GAMMA1P_M1_C3 + t * c; + c = INV_GAMMA1P_M1_C2 + t * c; + c = INV_GAMMA1P_M1_C1 + t * c; + c = INV_GAMMA1P_M1_C + t * c; + if (x > 0.5) { + ret = t * c / x; + } else { + ret = x * ((c + 0.5) + 0.5); + } + } else { + double p = INV_GAMMA1P_M1_P6; + p = INV_GAMMA1P_M1_P5 + t * p; + p = INV_GAMMA1P_M1_P4 + t * p; + p = INV_GAMMA1P_M1_P3 + t * p; + p = INV_GAMMA1P_M1_P2 + t * p; + p = INV_GAMMA1P_M1_P1 + t * p; + p = INV_GAMMA1P_M1_P0 + t * p; + + double q = INV_GAMMA1P_M1_Q4; + q = INV_GAMMA1P_M1_Q3 + t * q; + q = INV_GAMMA1P_M1_Q2 + t * q; + q = INV_GAMMA1P_M1_Q1 + t * q; + q = 1.0 + t * q; + + double c = INV_GAMMA1P_M1_C13 + (p / q) * t; + c = INV_GAMMA1P_M1_C12 + t * c; + c = INV_GAMMA1P_M1_C11 + t * c; + c = INV_GAMMA1P_M1_C10 + t * c; + c = INV_GAMMA1P_M1_C9 + t * c; + c = INV_GAMMA1P_M1_C8 + t * c; + c = INV_GAMMA1P_M1_C7 + t * c; + c = INV_GAMMA1P_M1_C6 + t * c; + c = INV_GAMMA1P_M1_C5 + t * c; + c = INV_GAMMA1P_M1_C4 + t * c; + c = INV_GAMMA1P_M1_C3 + t * c; + c = INV_GAMMA1P_M1_C2 + t * c; + c = INV_GAMMA1P_M1_C1 + t * c; + c = INV_GAMMA1P_M1_C0 + t * c; + + if (x > 0.5) { + ret = (t / x) * ((c - 0.5) - 0.5); + } else { + ret = x * c; + } + } + + return ret; + } + + /** + * Computes the function \( \ln \Gamma(1 + x) \) for \( -0.5 \leq x \leq 1.5 \). + * <p> + * This implementation is based on the double precision implementation in + * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}. + * </p> + * + * @param x Argument. + * @return \( \ln \Gamma(1 + x) \) + * @throws IllegalArgumentException if {@code x < -0.5}. + * @throws IllegalArgumentException if {@code x > 1.5}. + */ + public static double logGamma1p(final double x) { + + if (x < -0.5) { + throw new IllegalArgumentException(); + } + if (x > 1.5) { + throw new IllegalArgumentException(); + } + + return -Math.log1p(invGamma1pm1(x)); + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalUtils.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalUtils.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalUtils.java new file mode 100644 index 0000000..4fcf591 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InternalUtils.java @@ -0,0 +1,145 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ + +package org.apache.commons.rng.sampling.distribution; + +/** + * Functions used by some of the samplers. + * This class is not part of the public API, as it would be + * better to group these utilities in a dedicated components. + */ +class InternalUtils { // Class is package-private on purpose; do not make it public. + /** All long-representable factorials. */ + private static final long[] FACTORIALS = new long[] { + 1L, 1L, 2L, + 6L, 24L, 120L, + 720L, 5040L, 40320L, + 362880L, 3628800L, 39916800L, + 479001600L, 6227020800L, 87178291200L, + 1307674368000L, 20922789888000L, 355687428096000L, + 6402373705728000L, 121645100408832000L, 2432902008176640000L }; + + /** Utility class. */ + private InternalUtils() {} + + /** + * @param n Argument. + * @return {@code n!} + * @throws IllegalArgumentException if the result is too large to be represented + * by a {@code long} (i.e. if {@code n > 20}). + */ + public static long factorial(int n) { + if (n < 0) { + throw new IllegalArgumentException(n + " < " + 0); + } + if (n > 20) { + throw new IllegalArgumentException(n + " > " + 20); + } + return FACTORIALS[n]; + } + + /** + * Class for computing the natural logarithm of the factorial of {@code n}. + * It allows to allocate a cache of precomputed values. + * In case of cache miss, computation is preformed by a call to + * {@link Gamma#logGamma(double)}. + */ + public static final class FactorialLog { + /** + * Precomputed values of the function: + * {@code LOG_FACTORIALS[i] = log(i!)}. + */ + private final double[] LOG_FACTORIALS; + + /** + * Creates an instance, reusing the already computed values if available. + * + * @param numValues Number of values of the function to compute. + * @param cache Existing cache. + * @throw IllegalArgumentException if {@code numValues < 0}. + */ + private FactorialLog(int numValues, + double[] cache) { + if (numValues < 0) { + throw new IllegalArgumentException(numValues + " < " + 0); + } + + LOG_FACTORIALS = new double[numValues]; + + final int beginCopy = 2; + final int endCopy = cache == null || cache.length <= beginCopy ? + beginCopy : cache.length <= numValues ? + cache.length : numValues; + + // Copy available values. + for (int i = beginCopy; i < endCopy; i++) { + LOG_FACTORIALS[i] = cache[i]; + } + + // Precompute. + for (int i = endCopy; i < numValues; i++) { + LOG_FACTORIALS[i] = LOG_FACTORIALS[i - 1] + Math.log(i); + } + } + + /** + * Creates an instance with no precomputed values. + * @return an instance with no precomputed values. + */ + public static FactorialLog create() { + return new FactorialLog(0, null); + } + + /** + * Creates an instance with the specified cache size. + * + * @param cacheSize Number of precomputed values of the function. + * @return a new instance where {@code cacheSize} values have been + * precomputed. + * @throws IllegalArgumentException if {@code n < 0}. + */ + public FactorialLog withCache(final int cacheSize) { + return new FactorialLog(cacheSize, LOG_FACTORIALS); + } + + /** + * Computes {@code log(n!)}. + * + * @param n Argument. + * @return {@code log(n!)}. + * @throws IllegalArgumentException if {@code n < 0}. + */ + public double value(final int n) { + if (n < 0) { + throw new IllegalArgumentException(n + " < " + 0); + } + + // Use cache of precomputed values. + if (n < LOG_FACTORIALS.length) { + return LOG_FACTORIALS[n]; + } + + // Use cache of precomputed factorial values. + if (n < FACTORIALS.length) { + return Math.log(FACTORIALS[n]); + } + + // Delegate. + return InternalGamma.logGamma(n + 1); + } + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InverseMethodParetoSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InverseMethodParetoSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InverseMethodParetoSampler.java new file mode 100644 index 0000000..91a8bf9 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/InverseMethodParetoSampler.java @@ -0,0 +1,55 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ +package org.apache.commons.rng.sampling.distribution; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.AbstractContinuousSampler; + +/** + * Sampling from a <a href="https://en.wikipedia.org/wiki/Pareto_distribution";>Pareto distribution</a>. + */ +public class InverseMethodParetoSampler extends AbstractContinuousSampler { + /** Scale. */ + private final double scale; + /** Shape. */ + private final double shape; + + /** + * @param rng Generator of uniformly distributed random numbers. + * @param scale Scale of the distribution. + * @param shape Shape of the distribution. + */ + public InverseMethodParetoSampler(UniformRandomProvider rng, + double scale, + double shape) { + super(rng); + this.scale = scale; + this.shape = shape; + } + + /** {@inheritDoc} */ + @Override + public double sample() { + return scale / Math.pow(nextUniform(), 1 / shape); + } + + /** {@inheritDoc} */ + @Override + public String toString() { + return "[Inverse method for Pareto distribution " + super.toString() + "]"; + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/PoissonSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/PoissonSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/PoissonSampler.java new file mode 100644 index 0000000..6662068 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/PoissonSampler.java @@ -0,0 +1,164 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ +package org.apache.commons.rng.sampling.distribution; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.AbstractBaseSampler; +import org.apache.commons.rng.sampling.DiscreteSampler; +import org.apache.commons.rng.sampling.ContinuousSampler; + +/** + * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html";>Poisson distribution</a>. + */ +public class PoissonSampler + extends AbstractBaseSampler + implements DiscreteSampler { + /** Mean of the distribution. */ + private final double mean; + /** Exponential. */ + private final ContinuousSampler exponential; + /** Gaussian. */ + private final ContinuousSampler gaussian; + /** {@code log(n!)}. */ + private final InternalUtils.FactorialLog factorialLog; + + /** + * @param rng Generator of uniformly distributed random numbers. + * @param mean Mean. + * @throws IllegalArgumentException if {@code mean <= 0}. + */ + public PoissonSampler(UniformRandomProvider rng, + double mean) { + super(rng); + if (mean <= 0) { + throw new IllegalArgumentException(mean + " <= " + 0); + } + + this.mean = mean; + + gaussian = new BoxMullerGaussianSampler(rng, 0, 1); + exponential = new AhrensDieterExponentialSampler(rng, 1); + factorialLog = InternalUtils.FactorialLog.create(); + } + + /** {@inheritDoc} */ + @Override + public int sample() { + return (int) Math.min(nextPoisson(mean), Integer.MAX_VALUE); + } + + /** + * @param meanPoisson Mean. + * @return the next sample. + */ + private long nextPoisson(double meanPoisson) { + final double pivot = 40; + if (meanPoisson < pivot) { + double p = Math.exp(-meanPoisson); + long n = 0; + double r = 1; + double rnd = 1; + + while (n < 1000 * meanPoisson) { + rnd = nextUniform(); + r *= rnd; + if (r >= p) { + n++; + } else { + return n; + } + } + return n; + } else { + final double lambda = Math.floor(meanPoisson); + final double lambdaFractional = meanPoisson - lambda; + final double logLambda = Math.log(lambda); + final double logLambdaFactorial = factorialLog((int) lambda); + final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional); + final double delta = Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1)); + final double halfDelta = delta / 2; + final double twolpd = 2 * lambda + delta; + final double a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(1 / (8 * lambda)); + final double a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd); + final double aSum = a1 + a2 + 1; + final double p1 = a1 / aSum; + final double p2 = a2 / aSum; + final double c1 = 1 / (8 * lambda); + + double x = 0; + double y = 0; + double v = 0; + int a = 0; + double t = 0; + double qr = 0; + double qa = 0; + while (true) { + final double u = nextUniform(); + if (u <= p1) { + final double n = gaussian.sample(); + x = n * Math.sqrt(lambda + halfDelta) - 0.5d; + if (x > delta || x < -lambda) { + continue; + } + y = x < 0 ? Math.floor(x) : Math.ceil(x); + final double e = exponential.sample(); + v = -e - (n * n / 2) + c1; + } else { + if (u > p1 + p2) { + y = lambda; + break; + } else { + x = delta + (twolpd / delta) * exponential.sample(); + y = Math.ceil(x); + v = -exponential.sample() - delta * (x + 1) / twolpd; + } + } + a = x < 0 ? 1 : 0; + t = y * (y + 1) / (2 * lambda); + if (v < -t && a == 0) { + y = lambda + y; + break; + } + qr = t * ((2 * y + 1) / (6 * lambda) - 1); + qa = qr - (t * t) / (3 * (lambda + a * (y + 1))); + if (v < qa) { + y = lambda + y; + break; + } + if (v > qr) { + continue; + } + if (v < y * logLambda - factorialLog((int) (y + lambda)) + logLambdaFactorial) { + y = lambda + y; + break; + } + } + return y2 + (long) y; + } + } + + /** + * Compute the natural logarithm of the factorial of {@code n}. + * + * @param n Argument. + * @return {@code log(n!)} + * @throws IllegalArgumentException if {@code n < 0}. + */ + private double factorialLog(int n) { + return factorialLog.value(n); + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/74bbdd27/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/RejectionInversionZipfSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/RejectionInversionZipfSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/RejectionInversionZipfSampler.java new file mode 100644 index 0000000..de3ebb4 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/RejectionInversionZipfSampler.java @@ -0,0 +1,244 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You 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. + */ + +package org.apache.commons.rng.sampling.distribution; + +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.AbstractBaseSampler; +import org.apache.commons.rng.sampling.DiscreteSampler; + +/** + * Implementation of the <a href="https://en.wikipedia.org/wiki/Zipf's_law">Zipf distribution</a>. + * + * @since 1.0 + */ +public class RejectionInversionZipfSampler + extends AbstractBaseSampler + implements DiscreteSampler { + /** Number of elements. */ + private final int numberOfElements; + /** Exponent parameter of the distribution. */ + private final double exponent; + /** {@code hIntegral(1.5) - 1}. */ + private final double hIntegralX1; + /** {@code hIntegral(numberOfElements + 0.5)}. */ + private final double hIntegralNumberOfElements; + /** {@code 2 - hIntegralInverse(hIntegral(2.5) - h(2)}. */ + private final double s; + + /** + * @param rng Generator of uniformly distributed random numbers. + * @param numberOfElements Number of elements. + * @param exponent Exponent. + * @throws IllegalArgumentException if {@code numberOfElements <= 0} + * or {@code exponent <= 0}. + */ + public RejectionInversionZipfSampler(UniformRandomProvider rng, + int numberOfElements, + double exponent) { + super(rng); + if (numberOfElements <= 0) { + throw new IllegalArgumentException(numberOfElements + " <= " + 0); + } + if (exponent <= 0) { + throw new IllegalArgumentException(exponent + " <= " + 0); + } + + this.numberOfElements = numberOfElements; + this.exponent = exponent; + this.hIntegralX1 = hIntegral(1.5) - 1; + this.hIntegralNumberOfElements = hIntegral(numberOfElements + 0.5); + this.s = 2 - hIntegralInverse(hIntegral(2.5) - h(2)); + } + + /** + * Rejection inversion sampling method for a discrete, bounded Zipf + * distribution that is based on the method described in + * <p> + * Wolfgang Hörmann and Gerhard Derflinger + * "Rejection-inversion to generate variates from monotone discrete + * distributions", ACM Transactions on Modeling and Computer Simulation + * (TOMACS) 6.3 (1996): 169-184. + * </p> + * + * <p> + * The paper describes an algorithm for exponents larger than 1 + * (Algorithm ZRI). + * The original method uses {@code H(x) = (v + x)^(1 - q) / (1 - q)} + * as the integral of the hat function. + * This function is undefined for q = 1, which is the reason for + * the limitation of the exponent. + * If instead the integral function + * {@code H(x) = ((v + x)^(1 - q) - 1) / (1 - q)} is used, + * for which a meaningful limit exists for q = 1, the method works + * for all positive exponents. + * </p> + * + * <p> + * The following implementation uses v = 0 and generates integral + * number in the range [1, numberOfElements]. + * This is different to the original method where v is defined to + * be positive and numbers are taken from [0, i_max]. + * This explains why the implementation looks slightly different. + * </p> + */ + @Override + public int sample() { + while(true) { + final double u = hIntegralNumberOfElements + nextUniform() * (hIntegralX1 - hIntegralNumberOfElements); + // u is uniformly distributed in (hIntegralX1, hIntegralNumberOfElements] + + double x = hIntegralInverse(u); + int k = (int)(x + 0.5); + + // Limit k to the range [1, numberOfElements] + // (k could be outside due to numerical inaccuracies) + if (k < 1) { + k = 1; + } + else if (k > numberOfElements) { + k = numberOfElements; + } + + // Here, the distribution of k is given by: + // + // P(k = 1) = C * (hIntegral(1.5) - hIntegralX1) = C + // P(k = m) = C * (hIntegral(m + 1/2) - hIntegral(m - 1/2)) for m >= 2 + // + // where C = 1 / (hIntegralNumberOfElements - hIntegralX1) + + if (k - x <= s || u >= hIntegral(k + 0.5) - h(k)) { + + // Case k = 1: + // + // The right inequality is always true, because replacing k by 1 gives + // u >= hIntegral(1.5) - h(1) = hIntegralX1 and u is taken from + // (hIntegralX1, hIntegralNumberOfElements]. + // + // Therefore, the acceptance rate for k = 1 is P(accepted | k = 1) = 1 + // and the probability that 1 is returned as random value is + // P(k = 1 and accepted) = P(accepted | k = 1) * P(k = 1) = C = C / 1^exponent + // + // Case k >= 2: + // + // The left inequality (k - x <= s) is just a short cut + // to avoid the more expensive evaluation of the right inequality + // (u >= hIntegral(k + 0.5) - h(k)) in many cases. + // + // If the left inequality is true, the right inequality is also true: + // Theorem 2 in the paper is valid for all positive exponents, because + // the requirements h'(x) = -exponent/x^(exponent + 1) < 0 and + // (-1/hInverse'(x))'' = (1+1/exponent) * x^(1/exponent-1) >= 0 + // are both fulfilled. + // Therefore, f(x) = x - hIntegralInverse(hIntegral(x + 0.5) - h(x)) + // is a non-decreasing function. If k - x <= s holds, + // k - x <= s + f(k) - f(2) is obviously also true which is equivalent to + // -x <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)), + // -hIntegralInverse(u) <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)), + // and finally u >= hIntegral(k + 0.5) - h(k). + // + // Hence, the right inequality determines the acceptance rate: + // P(accepted | k = m) = h(m) / (hIntegrated(m+1/2) - hIntegrated(m-1/2)) + // The probability that m is returned is given by + // P(k = m and accepted) = P(accepted | k = m) * P(k = m) = C * h(m) = C / m^exponent. + // + // In both cases the probabilities are proportional to the probability mass function + // of the Zipf distribution. + + return k; + } + } + } + + /** + * {@code H(x)} is defined as + * <ul> + * <li>{@code (x^(1 - exponent) - 1) / (1 - exponent)}, if {@code exponent != 1}</li> + * <li>{@code log(x)}, if {@code exponent == 1}</li> + * </ul> + * H(x) is an integral function of h(x), the derivative of H(x) is h(x). + * + * @param x free parameter + * @return {@code H(x)}. + */ + private double hIntegral(final double x) { + final double logX = Math.log(x); + return helper2((1d-exponent)*logX)*logX; + } + + /** + * {@code h(x) = 1 / x^exponent} + * + * @param x Free parameter. + * @return {@code h(x)}. + */ + private double h(final double x) { + return Math.exp(-exponent * Math.log(x)); + } + + /** + * The inverse function of {@code H(x)}. + * + * @param x Free parameter + * @return y for which {@code H(y) = x}. + */ + private double hIntegralInverse(final double x) { + double t = x * (1 - exponent); + if (t < -1) { + // Limit value to the range [-1, +inf). + // t could be smaller than -1 in some rare cases due to numerical errors. + t = -1; + } + return Math.exp(helper1(t) * x); + } + + /** + * Helper function that calculates {@code log(1 + x) / x}. + * <p> + * A Taylor series expansion is used, if x is close to 0. + * </p> + * + * @param x A value larger than or equal to -1. + * @return {@code log(1 + x) / x}. + */ + private static double helper1(final double x) { + if (Math.abs(x) > 1e-8) { + return Math.log1p(x) / x; + } + else { + return 1 - x * (0.5 - x * (0.33333333333333333 - 0.25 * x)); + } + } + + /** + * Helper function to calculate {@code (exp(x) - 1) / x}. + * <p> + * A Taylor series expansion is used, if x is close to 0. + * </p> + * + * @param x Free parameter. + * @return {@code (exp(x) - 1) / x} if x is non-zero, or 1 if x = 0. + */ + private static double helper2(final double x) { + if (Math.abs(x) > 1e-8) { + return Math.expm1(x) / x; + } + else { + return 1 + x * 0.5 * (1 + x * 0.33333333333333333 * (1 + 0.25 * x)); + } + } +}
