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 a <a href="https://en.wikipedia.org/wiki/Geometric_distribution">geometric
023 * distribution</a>.
024 *
025 * <p>This distribution samples the number of failures before the first success taking values in the
026 * set {@code [0, 1, 2, ...]}.</p>
027 *
028 * <p>The sample is computed using a related exponential distribution. If \( X \) is an
029 * exponentially distributed random variable with parameter \( \lambda \), then
030 * \( Y = \left \lfloor X \right \rfloor \) is a geometrically distributed random variable with
031 * parameter \( p = 1 − e^\lambda \), with \( p \) the probability of success.</p>
032 *
033 * <p>This sampler outperforms using the {@link InverseTransformDiscreteSampler} with an appropriate
034 * geometric inverse cumulative probability function.</p>
035 *
036 * <p>Usage note: As the probability of success (\( p \)) tends towards zero the mean of the
037 * distribution (\( \frac{1-p}{p} \)) tends towards infinity and due to the use of {@code int}
038 * for the sample this can result in truncation of the distribution.</p>
039 *
040 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
041 *
042 * @see <a
043 * href="https://en.wikipedia.org/wiki/Geometric_distribution#Related_distributions">Geometric
044 * distribution - related distributions</a>
045 * @since 1.3
046 */
047public final class GeometricSampler {
048    /**
049     * Sample from the geometric distribution when the probability of success is 1.
050     */
051    private static final class GeometricP1Sampler
052        implements SharedStateDiscreteSampler {
053        /** The single instance. */
054        static final GeometricP1Sampler INSTANCE = new GeometricP1Sampler();
055
056        @Override
057        public int sample() {
058            // When probability of success is 1 the sample is always zero
059            return 0;
060        }
061
062        @Override
063        public String toString() {
064            return "Geometric(p=1) deviate";
065        }
066
067        @Override
068        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
069            // No requirement for a new instance
070            return this;
071        }
072    }
073
074    /**
075     * Sample from the geometric distribution by using a related exponential distribution.
076     */
077    private static final class GeometricExponentialSampler
078        implements SharedStateDiscreteSampler {
079        /** Underlying source of randomness. Used only for the {@link #toString()} method. */
080        private final UniformRandomProvider rng;
081        /** The related exponential sampler for the geometric distribution. */
082        private final SharedStateContinuousSampler exponentialSampler;
083
084        /**
085         * @param rng Generator of uniformly distributed random numbers
086         * @param probabilityOfSuccess The probability of success (must be in the range
087         * {@code [0 < probabilityOfSuccess < 1]})
088         */
089        GeometricExponentialSampler(UniformRandomProvider rng, double probabilityOfSuccess) {
090            this.rng = rng;
091            // Use a related exponential distribution:
092            // λ = −ln(1 − probabilityOfSuccess)
093            // exponential mean = 1 / λ
094            // --
095            // Note on validation:
096            // If probabilityOfSuccess == Math.nextDown(1.0) the exponential mean is >0 (valid).
097            // If probabilityOfSuccess == Double.MIN_VALUE the exponential mean is +Infinity
098            // and the sample will always be Integer.MAX_VALUE (the distribution is truncated). It
099            // is noted in the class Javadoc that the use of a small p leads to truncation so
100            // no checks are made for this case.
101            final double exponentialMean = 1.0 / (-Math.log1p(-probabilityOfSuccess));
102            exponentialSampler = ZigguratSampler.Exponential.of(rng, exponentialMean);
103        }
104
105        /**
106         * @param rng Generator of uniformly distributed random numbers
107         * @param source Source to copy.
108         */
109        GeometricExponentialSampler(UniformRandomProvider rng, GeometricExponentialSampler source) {
110            this.rng = rng;
111            exponentialSampler = source.exponentialSampler.withUniformRandomProvider(rng);
112        }
113
114        @Override
115        public int sample() {
116            // Return the floor of the exponential sample
117            return (int) Math.floor(exponentialSampler.sample());
118        }
119
120        @Override
121        public String toString() {
122            return "Geometric deviate [" + rng.toString() + "]";
123        }
124
125        @Override
126        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
127            return new GeometricExponentialSampler(rng, this);
128        }
129    }
130
131    /** Class contains only static methods. */
132    private GeometricSampler() {}
133
134    /**
135     * Creates a new geometric distribution sampler. The samples will be provided in the set
136     * {@code k=[0, 1, 2, ...]} where {@code k} indicates the number of failures before the first
137     * success.
138     *
139     * @param rng Generator of uniformly distributed random numbers.
140     * @param probabilityOfSuccess The probability of success.
141     * @return the sampler
142     * @throws IllegalArgumentException if {@code probabilityOfSuccess} is not in the range
143     * {@code [0 < probabilityOfSuccess <= 1]})
144     */
145    public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
146                                                double probabilityOfSuccess) {
147        if (probabilityOfSuccess <= 0 || probabilityOfSuccess > 1) {
148            throw new IllegalArgumentException(
149                "Probability of success (p) must be in the range [0 < p <= 1]: " +
150                    probabilityOfSuccess);
151        }
152        return probabilityOfSuccess == 1 ?
153            GeometricP1Sampler.INSTANCE :
154            new GeometricExponentialSampler(rng, probabilityOfSuccess);
155    }
156}