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:</p> 025 * 026 * <ul> 027 * <li>{@link UniformRandomProvider#nextLong()} 028 * <li>{@link UniformRandomProvider#nextDouble()} 029 * </ul> 030 * 031 * @since 1.0 032 */ 033public class AhrensDieterExponentialSampler 034 extends SamplerBase 035 implements SharedStateContinuousSampler { 036 /** 037 * Table containing the constants 038 * \( q_i = sum_{j=1}^i (\ln 2)^j / j! = \ln 2 + (\ln 2)^2 / 2 + ... + (\ln 2)^i / i! \) 039 * until the largest representable fraction below 1 is exceeded. 040 * 041 * Note that 042 * \( 1 = 2 - 1 = \exp(\ln 2) - 1 = sum_{n=1}^\infinity (\ln 2)^n / n! \) 043 * thus \( q_i \rightarrow 1 as i \rightarrow +\infinity \), 044 * so the higher \( i \), the closer we get to 1 (the series is not alternating). 045 * 046 * By trying, n = 16 in Java is enough to reach 1. 047 */ 048 private static final double[] EXPONENTIAL_SA_QI = new double[16]; 049 /** The mean of this distribution. */ 050 private final double mean; 051 /** Underlying source of randomness. */ 052 private final UniformRandomProvider rng; 053 054 // 055 // Initialize tables. 056 // 057 static { 058 // 059 // Filling EXPONENTIAL_SA_QI table. 060 // Note that we don't want qi = 0 in the table. 061 // 062 final double ln2 = Math.log(2); 063 double qi = 0; 064 065 // Start with 0! 066 // This will not overflow a long as the length < 21 067 long factorial = 1; 068 for (int i = 0; i < EXPONENTIAL_SA_QI.length; i++) { 069 factorial *= i + 1; 070 qi += Math.pow(ln2, i + 1.0) / factorial; 071 EXPONENTIAL_SA_QI[i] = qi; 072 } 073 } 074 075 /** 076 * Create an instance. 077 * 078 * @param rng Generator of uniformly distributed random numbers. 079 * @param mean Mean of this distribution. 080 * @throws IllegalArgumentException if {@code mean <= 0} 081 */ 082 public AhrensDieterExponentialSampler(UniformRandomProvider rng, 083 double mean) { 084 // Validation before java.lang.Object constructor exits prevents partially initialized object 085 this(InternalUtils.requireStrictlyPositive(mean, "mean"), rng); 086 } 087 088 /** 089 * @param mean Mean. 090 * @param rng Generator of uniformly distributed random numbers. 091 */ 092 private AhrensDieterExponentialSampler(double mean, 093 UniformRandomProvider rng) { 094 super(null); 095 this.rng = rng; 096 this.mean = mean; 097 } 098 099 /** {@inheritDoc} */ 100 @Override 101 public double sample() { 102 // Step 1: 103 double a = 0; 104 // Avoid u=0 which creates an infinite loop 105 double u = InternalUtils.makeNonZeroDouble(rng.nextLong()); 106 107 // Step 2 and 3: 108 while (u < 0.5) { 109 a += EXPONENTIAL_SA_QI[0]; 110 u *= 2; 111 } 112 113 // Step 4 (now u >= 0.5): 114 u += u - 1; 115 116 // Step 5: 117 if (u <= EXPONENTIAL_SA_QI[0]) { 118 return mean * (a + u); 119 } 120 121 // Step 6: 122 int i = 0; // Should be 1, be we iterate before it in while using 0. 123 double u2 = rng.nextDouble(); 124 double umin = u2; 125 126 // Step 7 and 8: 127 do { 128 ++i; 129 u2 = rng.nextDouble(); 130 131 if (u2 < umin) { 132 umin = u2; 133 } 134 135 // Step 8: 136 } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1. 137 138 return mean * (a + umin * EXPONENTIAL_SA_QI[0]); 139 } 140 141 /** {@inheritDoc} */ 142 @Override 143 public String toString() { 144 return "Ahrens-Dieter Exponential deviate [" + rng.toString() + "]"; 145 } 146 147 /** 148 * {@inheritDoc} 149 * 150 * @since 1.3 151 */ 152 @Override 153 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) { 154 // Use private constructor without validation 155 return new AhrensDieterExponentialSampler(mean, rng); 156 } 157 158 /** 159 * Create a new exponential distribution sampler. 160 * 161 * @param rng Generator of uniformly distributed random numbers. 162 * @param mean Mean of the distribution. 163 * @return the sampler 164 * @throws IllegalArgumentException if {@code mean <= 0} 165 * @since 1.3 166 */ 167 public static SharedStateContinuousSampler of(UniformRandomProvider rng, 168 double mean) { 169 return new AhrensDieterExponentialSampler(rng, mean); 170 } 171}