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="http://en.wikipedia.org/wiki/Beta_distribution">beta distribution</a>.
23   *
24   * <p>Uses Cheng's algorithms for beta distribution sampling:</p>
25   *
26   * <blockquote>
27   * <pre>
28   * R. C. H. Cheng,
29   * "Generating beta variates with nonintegral shape parameters",
30   * Communications of the ACM, 21, 317-322, 1978.
31   * </pre>
32   * </blockquote>
33   *
34   * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
35   *
36   * @since 1.0
37   */
38  public class ChengBetaSampler
39      extends SamplerBase
40      implements SharedStateContinuousSampler {
41      /** Natural logarithm of 4. */
42      private static final double LN_4 = Math.log(4.0);
43  
44      /** The appropriate beta sampler for the parameters. */
45      private final SharedStateContinuousSampler delegate;
46  
47      /**
48       * Base class to implement Cheng's algorithms for the beta distribution.
49       */
50      private abstract static class BaseChengBetaSampler
51              implements SharedStateContinuousSampler {
52          /**
53           * Flag set to true if {@code a} is the beta distribution alpha shape parameter.
54           * Otherwise {@code a} is the beta distribution beta shape parameter.
55           *
56           * <p>From the original Cheng paper this is equal to the result of: {@code a == a0}.</p>
57           */
58          protected final boolean aIsAlphaShape;
59          /**
60           * First shape parameter.
61           * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
62           */
63          protected final double a;
64          /**
65           * Second shape parameter.
66           * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
67           */
68          protected final double b;
69          /** Underlying source of randomness. */
70          protected final UniformRandomProvider rng;
71          /**
72           * The algorithm alpha factor. This is not the beta distribution alpha shape parameter.
73           * It is the sum of the two shape parameters ({@code a + b}.
74           */
75          protected final double alpha;
76          /** The logarithm of the alpha factor. */
77          protected final double logAlpha;
78  
79          /**
80           * @param rng Generator of uniformly distributed random numbers.
81           * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
82           * @param a Distribution first shape parameter.
83           * @param b Distribution second shape parameter.
84           */
85          BaseChengBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
86              this.rng = rng;
87              this.aIsAlphaShape = aIsAlphaShape;
88              this.a = a;
89              this.b = b;
90              alpha = a + b;
91              logAlpha = Math.log(alpha);
92          }
93  
94          /**
95           * @param rng Generator of uniformly distributed random numbers.
96           * @param source Source to copy.
97           */
98          private BaseChengBetaSampler(UniformRandomProvider rng,
99                                       BaseChengBetaSampler source) {
100             this.rng = rng;
101             aIsAlphaShape = source.aIsAlphaShape;
102             a = source.a;
103             b = source.b;
104             // Compute algorithm factors.
105             alpha = source.alpha;
106             logAlpha = source.logAlpha;
107         }
108 
109         /** {@inheritDoc} */
110         @Override
111         public String toString() {
112             return "Cheng Beta deviate [" + rng.toString() + "]";
113         }
114 
115         /**
116          * Compute the sample result X.
117          *
118          * <blockquote>
119          * If a == a0 deliver X = W/(b + W); otherwise deliver X = b/(b + W).
120          * </blockquote>
121          *
122          * <p>The finalisation step is shared between the BB and BC algorithm (as step 5 of the
123          * BB algorithm and step 6 of the BC algorithm).</p>
124          *
125          * @param w Algorithm value W.
126          * @return the sample value
127          */
128         protected double computeX(double w) {
129             // Avoid (infinity / infinity) producing NaN
130             final double tmp = Math.min(w, Double.MAX_VALUE);
131             return aIsAlphaShape ? tmp / (b + tmp) : b / (b + tmp);
132         }
133     }
134 
135     /**
136      * Computes one sample using Cheng's BB algorithm, when beta distribution {@code alpha} and
137      * {@code beta} shape parameters are both larger than 1.
138      */
139     private static class ChengBBBetaSampler extends BaseChengBetaSampler {
140         /** 1 + natural logarithm of 5. */
141         private static final double LN_5_P1 = 1 + Math.log(5.0);
142 
143         /** The algorithm beta factor. This is not the beta distribution beta shape parameter. */
144         private final double beta;
145         /** The algorithm gamma factor. */
146         private final double gamma;
147 
148         /**
149          * @param rng Generator of uniformly distributed random numbers.
150          * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
151          * @param a min(alpha, beta) shape parameter.
152          * @param b max(alpha, beta) shape parameter.
153          */
154         ChengBBBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
155             super(rng, aIsAlphaShape, a, b);
156             beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
157             gamma = a + 1 / beta;
158         }
159 
160         /**
161          * @param rng Generator of uniformly distributed random numbers.
162          * @param source Source to copy.
163          */
164         private ChengBBBetaSampler(UniformRandomProvider rng,
165                                    ChengBBBetaSampler source) {
166             super(rng, source);
167             // Compute algorithm factors.
168             beta = source.beta;
169             gamma = source.gamma;
170         }
171 
172         @Override
173         public double sample() {
174             double r;
175             double w;
176             double t;
177             do {
178                 // Step 1:
179                 final double u1 = rng.nextDouble();
180                 final double u2 = rng.nextDouble();
181                 final double v = beta * (Math.log(u1) - Math.log1p(-u1));
182                 w = a * Math.exp(v);
183                 final double z = u1 * u1 * u2;
184                 r = gamma * v - LN_4;
185                 final double s = a + r - w;
186                 // Step 2:
187                 if (s + LN_5_P1 >= 5 * z) {
188                     break;
189                 }
190 
191                 // Step 3:
192                 t = Math.log(z);
193                 if (s >= t) {
194                     break;
195                 }
196                 // Step 4:
197             } while (r + alpha * (logAlpha - Math.log(b + w)) < t);
198 
199             // Step 5:
200             return computeX(w);
201         }
202 
203         @Override
204         public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
205             return new ChengBBBetaSampler(rng, this);
206         }
207     }
208 
209     /**
210      * Computes one sample using Cheng's BC algorithm, when at least one of beta distribution
211      * {@code alpha} or {@code beta} shape parameters is smaller than 1.
212      */
213     private static class ChengBCBetaSampler extends BaseChengBetaSampler {
214         /** 1/2. */
215         private static final double ONE_HALF = 1d / 2;
216         /** 1/4. */
217         private static final double ONE_QUARTER = 1d / 4;
218 
219         /** The algorithm beta factor. This is not the beta distribution beta shape parameter. */
220         private final double beta;
221         /** The algorithm delta factor. */
222         private final double delta;
223         /** The algorithm k1 factor. */
224         private final double k1;
225         /** The algorithm k2 factor. */
226         private final double k2;
227 
228         /**
229          * @param rng Generator of uniformly distributed random numbers.
230          * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
231          * @param a max(alpha, beta) shape parameter.
232          * @param b min(alpha, beta) shape parameter.
233          */
234         ChengBCBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
235             super(rng, aIsAlphaShape, a, b);
236             // Compute algorithm factors.
237             beta = 1 / b;
238             delta = 1 + a - b;
239             // The following factors are divided by 4:
240             // k1 = delta * (1 + 3b) / (18a/b - 14)
241             // k2 = 1 + (2 + 1/delta) * b
242             k1 = delta * (1.0 / 72.0 + 3.0 / 72.0 * b) / (a * beta - 7.0 / 9.0);
243             k2 = 0.25 + (0.5 + 0.25 / delta) * b;
244         }
245 
246         /**
247          * @param rng Generator of uniformly distributed random numbers.
248          * @param source Source to copy.
249          */
250         private ChengBCBetaSampler(UniformRandomProvider rng,
251                                    ChengBCBetaSampler source) {
252             super(rng, source);
253             beta = source.beta;
254             delta = source.delta;
255             k1 = source.k1;
256             k2 = source.k2;
257         }
258 
259         @Override
260         public double sample() {
261             double w;
262             while (true) {
263                 // Step 1:
264                 final double u1 = rng.nextDouble();
265                 final double u2 = rng.nextDouble();
266                 // Compute Y and Z
267                 final double y = u1 * u2;
268                 final double z = u1 * y;
269                 if (u1 < ONE_HALF) {
270                     // Step 2:
271                     if (ONE_QUARTER * u2 + z - y >= k1) {
272                         continue;
273                     }
274                 } else {
275                     // Step 3:
276                     if (z <= ONE_QUARTER) {
277                         final double v = beta * (Math.log(u1) - Math.log1p(-u1));
278                         w = a * Math.exp(v);
279                         break;
280                     }
281 
282                     // Step 4:
283                     if (z >= k2) {
284                         continue;
285                     }
286                 }
287 
288                 // Step 5:
289                 final double v = beta * (Math.log(u1) - Math.log1p(-u1));
290                 w = a * Math.exp(v);
291                 if (alpha * (logAlpha - Math.log(b + w) + v) - LN_4 >= Math.log(z)) {
292                     break;
293                 }
294             }
295 
296             // Step 6:
297             return computeX(w);
298         }
299 
300         @Override
301         public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
302             return new ChengBCBetaSampler(rng, this);
303         }
304     }
305 
306     /**
307      * Creates a sampler instance.
308      *
309      * @param rng Generator of uniformly distributed random numbers.
310      * @param alpha Distribution first shape parameter.
311      * @param beta Distribution second shape parameter.
312      * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}
313      */
314     public ChengBetaSampler(UniformRandomProvider rng,
315                             double alpha,
316                             double beta) {
317         super(null);
318         delegate = of(rng, alpha, beta);
319     }
320 
321     /** {@inheritDoc} */
322     @Override
323     public double sample() {
324         return delegate.sample();
325     }
326 
327     /** {@inheritDoc} */
328     @Override
329     public String toString() {
330         return delegate.toString();
331     }
332 
333     /**
334      * {@inheritDoc}
335      *
336      * @since 1.3
337      */
338     @Override
339     public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
340         return delegate.withUniformRandomProvider(rng);
341     }
342 
343     /**
344      * Creates a new beta distribution sampler.
345      *
346      * @param rng Generator of uniformly distributed random numbers.
347      * @param alpha Distribution first shape parameter.
348      * @param beta Distribution second shape parameter.
349      * @return the sampler
350      * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}
351      * @since 1.3
352      */
353     public static SharedStateContinuousSampler of(UniformRandomProvider rng,
354                                                   double alpha,
355                                                   double beta) {
356         if (alpha <= 0) {
357             throw new IllegalArgumentException("alpha is not strictly positive: " + alpha);
358         }
359         if (beta <= 0) {
360             throw new IllegalArgumentException("beta is not strictly positive: " + beta);
361         }
362 
363         // Choose the algorithm.
364         final double a = Math.min(alpha, beta);
365         final double b = Math.max(alpha, beta);
366         final boolean aIsAlphaShape = a == alpha;
367 
368         return a > 1 ?
369             // BB algorithm
370             new ChengBBBetaSampler(rng, aIsAlphaShape, a, b) :
371             // The BC algorithm is deliberately invoked with reversed parameters
372             // as the argument order is: max(alpha,beta), min(alpha,beta).
373             // Also invert the 'a is alpha' flag.
374             new ChengBCBetaSampler(rng, !aIsAlphaShape, b, a);
375     }
376 }