001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.sampling.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Sampling from an <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">exponential distribution</a>.
023 *
024 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
025 *
026 * @since 1.0
027 */
028public class AhrensDieterExponentialSampler
029    extends SamplerBase
030    implements SharedStateContinuousSampler {
031    /**
032     * Table containing the constants
033     * \( q_i = sum_{j=1}^i (\ln 2)^j / j! = \ln 2 + (\ln 2)^2 / 2 + ... + (\ln 2)^i / i! \)
034     * until the largest representable fraction below 1 is exceeded.
035     *
036     * Note that
037     * \( 1 = 2 - 1 = \exp(\ln 2) - 1 = sum_{n=1}^\infinity (\ln 2)^n / n! \)
038     * thus \( q_i \rightarrow 1 as i \rightarrow +\infinity \),
039     * so the higher \( i \), the closer we get to 1 (the series is not alternating).
040     *
041     * By trying, n = 16 in Java is enough to reach 1.
042     */
043    private static final double[] EXPONENTIAL_SA_QI = new double[16];
044    /** The mean of this distribution. */
045    private final double mean;
046    /** Underlying source of randomness. */
047    private final UniformRandomProvider rng;
048
049    /**
050     * Initialize tables.
051     */
052    static {
053        /**
054         * Filling EXPONENTIAL_SA_QI table.
055         * Note that we don't want qi = 0 in the table.
056         */
057        final double ln2 = Math.log(2);
058        double qi = 0;
059
060        for (int i = 0; i < EXPONENTIAL_SA_QI.length; i++) {
061            qi += Math.pow(ln2, i + 1.0) / InternalUtils.factorial(i + 1);
062            EXPONENTIAL_SA_QI[i] = qi;
063        }
064    }
065
066    /**
067     * @param rng Generator of uniformly distributed random numbers.
068     * @param mean Mean of this distribution.
069     * @throws IllegalArgumentException if {@code mean <= 0}
070     */
071    public AhrensDieterExponentialSampler(UniformRandomProvider rng,
072                                          double mean) {
073        super(null);
074        if (mean <= 0) {
075            throw new IllegalArgumentException("mean is not strictly positive: " + mean);
076        }
077        this.rng = rng;
078        this.mean = mean;
079    }
080
081    /**
082     * @param rng Generator of uniformly distributed random numbers.
083     * @param source Source to copy.
084     */
085    private AhrensDieterExponentialSampler(UniformRandomProvider rng,
086                                           AhrensDieterExponentialSampler source) {
087        super(null);
088        this.rng = rng;
089        this.mean = source.mean;
090    }
091
092    /** {@inheritDoc} */
093    @Override
094    public double sample() {
095        // Step 1:
096        double a = 0;
097        double u = rng.nextDouble();
098
099        // Step 2 and 3:
100        while (u < 0.5) {
101            a += EXPONENTIAL_SA_QI[0];
102            u *= 2;
103        }
104
105        // Step 4 (now u >= 0.5):
106        u += u - 1;
107
108        // Step 5:
109        if (u <= EXPONENTIAL_SA_QI[0]) {
110            return mean * (a + u);
111        }
112
113        // Step 6:
114        int i = 0; // Should be 1, be we iterate before it in while using 0.
115        double u2 = rng.nextDouble();
116        double umin = u2;
117
118        // Step 7 and 8:
119        do {
120            ++i;
121            u2 = rng.nextDouble();
122
123            if (u2 < umin) {
124                umin = u2;
125            }
126
127            // Step 8:
128        } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1.
129
130        return mean * (a + umin * EXPONENTIAL_SA_QI[0]);
131    }
132
133    /** {@inheritDoc} */
134    @Override
135    public String toString() {
136        return "Ahrens-Dieter Exponential deviate [" + rng.toString() + "]";
137    }
138
139    /**
140     * {@inheritDoc}
141     *
142     * @since 1.3
143     */
144    @Override
145    public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
146        return new AhrensDieterExponentialSampler(rng, this);
147    }
148
149    /**
150     * Create a new exponential distribution sampler.
151     *
152     * @param rng Generator of uniformly distributed random numbers.
153     * @param mean Mean of the distribution.
154     * @return the sampler
155     * @throws IllegalArgumentException if {@code mean <= 0}
156     * @since 1.3
157     */
158    public static SharedStateContinuousSampler of(UniformRandomProvider rng,
159                                                  double mean) {
160        return new AhrensDieterExponentialSampler(rng, mean);
161    }
162}