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 exponential distribution. 023 * 024 * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a> 025 */ 026public class AhrensDieterExponentialSampler 027 extends SamplerBase 028 implements ContinuousSampler { 029 /** 030 * Table containing the constants 031 * \( q_i = sum_{j=1}^i (\ln 2)^j / j! = \ln 2 + (\ln 2)^2 / 2 + ... + (\ln 2)^i / i! \) 032 * until the largest representable fraction below 1 is exceeded. 033 * 034 * Note that 035 * \( 1 = 2 - 1 = \exp(\ln 2) - 1 = sum_{n=1}^\infinity (\ln 2)^n / n! \) 036 * thus \( q_i \rightarrow 1 as i \rightarrow +\infinity \), 037 * so the higher \( i \), the closer we get to 1 (the series is not alternating). 038 * 039 * By trying, n = 16 in Java is enough to reach 1. 040 */ 041 private static final double[] EXPONENTIAL_SA_QI = new double[16]; 042 /** The mean of this distribution. */ 043 private final double mean; 044 045 /** 046 * Initialize tables. 047 */ 048 static { 049 /** 050 * Filling EXPONENTIAL_SA_QI table. 051 * Note that we don't want qi = 0 in the table. 052 */ 053 final double ln2 = Math.log(2); 054 double qi = 0; 055 056 for (int i = 0; i < EXPONENTIAL_SA_QI.length; i++) { 057 qi += Math.pow(ln2, i + 1) / InternalUtils.factorial(i + 1); 058 EXPONENTIAL_SA_QI[i] = qi; 059 } 060 } 061 062 /** 063 * @param rng Generator of uniformly distributed random numbers. 064 * @param mean Mean of this distribution. 065 */ 066 public AhrensDieterExponentialSampler(UniformRandomProvider rng, 067 double mean) { 068 super(rng); 069 this.mean = mean; 070 } 071 072 /** {@inheritDoc} */ 073 @Override 074 public double sample() { 075 // Step 1: 076 double a = 0; 077 double u = nextDouble(); 078 079 // Step 2 and 3: 080 while (u < 0.5) { 081 a += EXPONENTIAL_SA_QI[0]; 082 u *= 2; 083 } 084 085 // Step 4 (now u >= 0.5): 086 u += u - 1; 087 088 // Step 5: 089 if (u <= EXPONENTIAL_SA_QI[0]) { 090 return mean * (a + u); 091 } 092 093 // Step 6: 094 int i = 0; // Should be 1, be we iterate before it in while using 0. 095 double u2 = nextDouble(); 096 double umin = u2; 097 098 // Step 7 and 8: 099 do { 100 ++i; 101 u2 = nextDouble(); 102 103 if (u2 < umin) { 104 umin = u2; 105 } 106 107 // Step 8: 108 } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1. 109 110 return mean * (a + umin * EXPONENTIAL_SA_QI[0]); 111 } 112 113 /** {@inheritDoc} */ 114 @Override 115 public String toString() { 116 return "Ahrens-Dieter Exponential deviate [" + super.toString() + "]"; 117 } 118}