View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.rng.sampling.distribution;
18  
19  import org.apache.commons.rng.UniformRandomProvider;
20  
21  /**
22   * Sampling from a <a href="https://en.wikipedia.org/wiki/Geometric_distribution">geometric
23   * distribution</a>.
24   *
25   * <p>This distribution samples the number of failures before the first success taking values in the
26   * set {@code [0, 1, 2, ...]}.</p>
27   *
28   * <p>The sample is computed using a related exponential distribution. If \( X \) is an
29   * exponentially distributed random variable with parameter \( \lambda \), then
30   * \( Y = \left \lfloor X \right \rfloor \) is a geometrically distributed random variable with
31   * parameter \( p = 1 − e^\lambda \), with \( p \) the probability of success.</p>
32   *
33   * <p>This sampler outperforms using the {@link InverseTransformDiscreteSampler} with an appropriate
34   * geometric inverse cumulative probability function.</p>
35   *
36   * <p>Usage note: As the probability of success (\( p \)) tends towards zero the mean of the
37   * distribution (\( \frac{1-p}{p} \)) tends towards infinity and due to the use of {@code int}
38   * for the sample this can result in truncation of the distribution.</p>
39   *
40   * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
41   *
42   * @see <a
43   * href="https://en.wikipedia.org/wiki/Geometric_distribution#Related_distributions">Geometric
44   * distribution - related distributions</a>
45   * @since 1.3
46   */
47  public final class GeometricSampler {
48      /**
49       * Sample from the geometric distribution when the probability of success is 1.
50       */
51      private static class GeometricP1Sampler
52          implements SharedStateDiscreteSampler {
53          /** The single instance. */
54          static final GeometricP1Sampler INSTANCE = new GeometricP1Sampler();
55  
56          @Override
57          public int sample() {
58              // When probability of success is 1 the sample is always zero
59              return 0;
60          }
61  
62          @Override
63          public String toString() {
64              return "Geometric(p=1) deviate";
65          }
66  
67          @Override
68          public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
69              // No requirement for a new instance
70              return this;
71          }
72      }
73  
74      /**
75       * Sample from the geometric distribution by using a related exponential distribution.
76       */
77      private static class GeometricExponentialSampler
78          implements SharedStateDiscreteSampler {
79          /** Underlying source of randomness. Used only for the {@link #toString()} method. */
80          private final UniformRandomProvider rng;
81          /** The related exponential sampler for the geometric distribution. */
82          private final SharedStateContinuousSampler exponentialSampler;
83  
84          /**
85           * @param rng Generator of uniformly distributed random numbers
86           * @param probabilityOfSuccess The probability of success (must be in the range
87           * {@code [0 < probabilityOfSuccess < 1]})
88           */
89          GeometricExponentialSampler(UniformRandomProvider rng, double probabilityOfSuccess) {
90              this.rng = rng;
91              // Use a related exponential distribution:
92              // λ = −ln(1 − probabilityOfSuccess)
93              // exponential mean = 1 / λ
94              // --
95              // Note on validation:
96              // If probabilityOfSuccess == Math.nextDown(1.0) the exponential mean is >0 (valid).
97              // If probabilityOfSuccess == Double.MIN_VALUE the exponential mean is +Infinity
98              // and the sample will always be Integer.MAX_VALUE (the distribution is truncated). It
99              // 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 }