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 an <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">exponential distribution</a>. 23 * 24 * <p>Sampling uses:</p> 25 * 26 * <ul> 27 * <li>{@link UniformRandomProvider#nextLong()} 28 * <li>{@link UniformRandomProvider#nextDouble()} 29 * </ul> 30 * 31 * @since 1.0 32 */ 33 public class AhrensDieterExponentialSampler 34 extends SamplerBase 35 implements SharedStateContinuousSampler { 36 /** 37 * Table containing the constants 38 * \( q_i = sum_{j=1}^i (\ln 2)^j / j! = \ln 2 + (\ln 2)^2 / 2 + ... + (\ln 2)^i / i! \) 39 * until the largest representable fraction below 1 is exceeded. 40 * 41 * Note that 42 * \( 1 = 2 - 1 = \exp(\ln 2) - 1 = sum_{n=1}^\infinity (\ln 2)^n / n! \) 43 * thus \( q_i \rightarrow 1 as i \rightarrow +\infinity \), 44 * so the higher \( i \), the closer we get to 1 (the series is not alternating). 45 * 46 * By trying, n = 16 in Java is enough to reach 1. 47 */ 48 private static final double[] EXPONENTIAL_SA_QI = new double[16]; 49 /** The mean of this distribution. */ 50 private final double mean; 51 /** Underlying source of randomness. */ 52 private final UniformRandomProvider rng; 53 54 // 55 // Initialize tables. 56 // 57 static { 58 // 59 // Filling EXPONENTIAL_SA_QI table. 60 // Note that we don't want qi = 0 in the table. 61 // 62 final double ln2 = Math.log(2); 63 double qi = 0; 64 65 for (int i = 0; i < EXPONENTIAL_SA_QI.length; i++) { 66 qi += Math.pow(ln2, i + 1.0) / InternalUtils.factorial(i + 1); 67 EXPONENTIAL_SA_QI[i] = qi; 68 } 69 } 70 71 /** 72 * @param rng Generator of uniformly distributed random numbers. 73 * @param mean Mean of this distribution. 74 * @throws IllegalArgumentException if {@code mean <= 0} 75 */ 76 public AhrensDieterExponentialSampler(UniformRandomProvider rng, 77 double mean) { 78 super(null); 79 if (mean <= 0) { 80 throw new IllegalArgumentException("mean is not strictly positive: " + mean); 81 } 82 this.rng = rng; 83 this.mean = mean; 84 } 85 86 /** 87 * @param rng Generator of uniformly distributed random numbers. 88 * @param source Source to copy. 89 */ 90 private AhrensDieterExponentialSampler(UniformRandomProvider rng, 91 AhrensDieterExponentialSampler source) { 92 super(null); 93 this.rng = rng; 94 this.mean = source.mean; 95 } 96 97 /** {@inheritDoc} */ 98 @Override 99 public double sample() { 100 // Step 1: 101 double a = 0; 102 // Avoid u=0 which creates an infinite loop 103 double u = InternalUtils.makeNonZeroDouble(rng.nextLong()); 104 105 // Step 2 and 3: 106 while (u < 0.5) { 107 a += EXPONENTIAL_SA_QI[0]; 108 u *= 2; 109 } 110 111 // Step 4 (now u >= 0.5): 112 u += u - 1; 113 114 // Step 5: 115 if (u <= EXPONENTIAL_SA_QI[0]) { 116 return mean * (a + u); 117 } 118 119 // Step 6: 120 int i = 0; // Should be 1, be we iterate before it in while using 0. 121 double u2 = rng.nextDouble(); 122 double umin = u2; 123 124 // Step 7 and 8: 125 do { 126 ++i; 127 u2 = rng.nextDouble(); 128 129 if (u2 < umin) { 130 umin = u2; 131 } 132 133 // Step 8: 134 } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1. 135 136 return mean * (a + umin * EXPONENTIAL_SA_QI[0]); 137 } 138 139 /** {@inheritDoc} */ 140 @Override 141 public String toString() { 142 return "Ahrens-Dieter Exponential deviate [" + rng.toString() + "]"; 143 } 144 145 /** 146 * {@inheritDoc} 147 * 148 * @since 1.3 149 */ 150 @Override 151 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) { 152 return new AhrensDieterExponentialSampler(rng, this); 153 } 154 155 /** 156 * Create a new exponential distribution sampler. 157 * 158 * @param rng Generator of uniformly distributed random numbers. 159 * @param mean Mean of the distribution. 160 * @return the sampler 161 * @throws IllegalArgumentException if {@code mean <= 0} 162 * @since 1.3 163 */ 164 public static SharedStateContinuousSampler of(UniformRandomProvider rng, 165 double mean) { 166 return new AhrensDieterExponentialSampler(rng, mean); 167 } 168 }