1/* boost random/exponential_distribution.hpp header file
2 *
3 * Copyright Jens Maurer 2000-2001
4 * Copyright Steven Watanabe 2011
5 * Copyright Jason Rhinelander 2016
6 * Distributed under the Boost Software License, Version 1.0. (See
7 * accompanying file LICENSE_1_0.txt or copy at
8 * http://www.boost.org/LICENSE_1_0.txt)
9 *
10 * See http://www.boost.org for most recent version including documentation.
11 *
12 * $Id$
13 *
14 * Revision history
15 * 2001-02-18 moved to individual header files
16 */
17
18#ifndef BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
19#define BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
20
21#include <boost/config/no_tr1/cmath.hpp>
22#include <iosfwd>
23#include <boost/assert.hpp>
24#include <boost/limits.hpp>
25#include <boost/random/detail/config.hpp>
26#include <boost/random/detail/operators.hpp>
27#include <boost/random/detail/int_float_pair.hpp>
28#include <boost/random/uniform_01.hpp>
29
30namespace boost {
31namespace random {
32
33namespace detail {
34
35// tables for the ziggurat algorithm
36template<class RealType>
37struct exponential_table {
38 static const RealType table_x[257];
39 static const RealType table_y[257];
40};
41
42template<class RealType>
43const RealType exponential_table<RealType>::table_x[257] = {
44 8.6971174701310497140, 7.6971174701310497140, 6.9410336293772123602, 6.4783784938325698538,
45 6.1441646657724730491, 5.8821443157953997963, 5.6664101674540337371, 5.4828906275260628694,
46 5.3230905057543986131, 5.1814872813015010392, 5.0542884899813047117, 4.9387770859012514838,
47 4.8329397410251125881, 4.7352429966017412526, 4.6444918854200854873, 4.5597370617073515513,
48 4.4802117465284221949, 4.4052876934735729805, 4.3344436803172730116, 4.2672424802773661873,
49 4.2033137137351843802, 4.1423408656640511251, 4.0840513104082974638, 4.0282085446479365106,
50 3.9746060666737884793, 3.9230625001354895926, 3.8734176703995089983, 3.8255294185223367372,
51 3.7792709924116678992, 3.7345288940397975350, 3.6912010902374189454, 3.6491955157608538478,
52 3.6084288131289096339, 3.5688252656483374051, 3.5303158891293438633, 3.4928376547740601814,
53 3.4563328211327607625, 3.4207483572511205323, 3.3860354424603017887, 3.3521490309001100106,
54 3.3190474709707487166, 3.2866921715990692095, 3.2550473085704501813, 3.2240795652862645207,
55 3.1937579032122407483, 3.1640533580259734580, 3.1349388580844407393, 3.1063890623398246660,
56 3.0783802152540905188, 3.0508900166154554479, 3.0238975044556767713, 2.9973829495161306949,
57 2.9713277599210896472, 2.9457143948950456386, 2.9205262865127406647, 2.8957477686001416838,
58 2.8713640120155362592, 2.8473609656351888266, 2.8237253024500354905, 2.8004443702507381944,
59 2.7775061464397572041, 2.7548991965623453650, 2.7326126361947007411, 2.7106360958679293686,
60 2.6889596887418041593, 2.6675739807732670816, 2.6464699631518093905, 2.6256390267977886123,
61 2.6050729387408355373, 2.5847638202141406911, 2.5647041263169053687, 2.5448866271118700928,
62 2.5253043900378279427, 2.5059507635285939648, 2.4868193617402096807, 2.4679040502973649846,
63 2.4491989329782498908, 2.4306983392644199088, 2.4123968126888708336, 2.3942890999214583288,
64 2.3763701405361408194, 2.3586350574093374601, 2.3410791477030346875, 2.3236978743901964559,
65 2.3064868582835798692, 2.2894418705322694265, 2.2725588255531546952, 2.2558337743672190441,
66 2.2392628983129087111, 2.2228425031110364013, 2.2065690132576635755, 2.1904389667232199235,
67 2.1744490099377744673, 2.1585958930438856781, 2.1428764653998416425, 2.1272876713173679737,
68 2.1118265460190418108, 2.0964902118017147637, 2.0812758743932248696, 2.0661808194905755036,
69 2.0512024094685848641, 2.0363380802487695916, 2.0215853383189260770, 2.0069417578945183144,
70 1.9924049782135764992, 1.9779727009573602295, 1.9636426877895480401, 1.9494127580071845659,
71 1.9352807862970511135, 1.9212447005915276767, 1.9073024800183871196, 1.8934521529393077332,
72 1.8796917950722108462, 1.8660195276928275962, 1.8524335159111751661, 1.8389319670188793980,
73 1.8255131289035192212, 1.8121752885263901413, 1.7989167704602903934, 1.7857359354841254047,
74 1.7726311792313049959, 1.7596009308890742369, 1.7466436519460739352, 1.7337578349855711926,
75 1.7209420025219350428, 1.7081947058780575683, 1.6955145241015377061, 1.6829000629175537544,
76 1.6703499537164519163, 1.6578628525741725325, 1.6454374393037234057, 1.6330724165359912048,
77 1.6207665088282577216, 1.6085184617988580769, 1.5963270412864831349, 1.5841910325326886695,
78 1.5721092393862294810, 1.5600804835278879161, 1.5481036037145133070, 1.5361774550410318943,
79 1.5243009082192260050, 1.5124728488721167573, 1.5006921768428164936, 1.4889578055167456003,
80 1.4772686611561334579, 1.4656236822457450411, 1.4540218188487932264, 1.4424620319720121876,
81 1.4309432929388794104, 1.4194645827699828254, 1.4080248915695353509, 1.3966232179170417110,
82 1.3852585682631217189, 1.3739299563284902176, 1.3626364025050864742, 1.3513769332583349176,
83 1.3401505805295045843, 1.3289563811371163220, 1.3177933761763245480, 1.3066606104151739482,
84 1.2955571316866007210, 1.2844819902750125450, 1.2734342382962410994, 1.2624129290696153434,
85 1.2514171164808525098, 1.2404458543344064544, 1.2294981956938491599, 1.2185731922087903071,
86 1.2076698934267612830, 1.1967873460884031665, 1.1859245934042023557, 1.1750806743109117687,
87 1.1642546227056790397, 1.1534454666557748056, 1.1426522275816728928, 1.1318739194110786733,
88 1.1211095477013306083, 1.1103581087274114281, 1.0996185885325976575, 1.0888899619385472598,
89 1.0781711915113727024, 1.0674612264799681530, 1.0567590016025518414, 1.0460634359770445503,
90 1.0353734317905289496, 1.0246878730026178052, 1.0140056239570971074, 1.0033255279156973717,
91 0.99264640550727647009, 0.98196705308506317914, 0.97128624098390397896, 0.96060271166866709917,
92 0.94991517776407659940, 0.93922231995526297952, 0.92852278474721113999, 0.91781518207004493915,
93 0.90709808271569100600, 0.89637001558989069006, 0.88562946476175228052, 0.87487486629102585352,
94 0.86410460481100519511, 0.85331700984237406386, 0.84251035181036928333, 0.83168283773427388393,
95 0.82083260655441252290, 0.80995772405741906620, 0.79905617735548788109, 0.78812586886949324977,
96 0.77716460975913043936, 0.76617011273543541328, 0.75513998418198289808, 0.74407171550050873971,
97 0.73296267358436604916, 0.72181009030875689912, 0.71061105090965570413, 0.69936248110323266174,
98 0.68806113277374858613, 0.67670356802952337911, 0.66528614139267855405, 0.65380497984766565353,
99 0.64225596042453703448, 0.63063468493349100113, 0.61893645139487678178, 0.60715622162030085137,
100 0.59528858429150359384, 0.58332771274877027785, 0.57126731653258903915, 0.55910058551154127652,
101 0.54682012516331112550, 0.53441788123716615385, 0.52188505159213564105, 0.50921198244365495319,
102 0.49638804551867159754, 0.48340149165346224782, 0.47023927508216945338, 0.45688684093142071279,
103 0.44332786607355296305, 0.42954394022541129589, 0.41551416960035700100, 0.40121467889627836229,
104 0.38661797794112021568, 0.37169214532991786118, 0.35639976025839443721, 0.34069648106484979674,
105 0.32452911701691008547, 0.30783295467493287307, 0.29052795549123115167, 0.27251318547846547924,
106 0.25365836338591284433, 0.23379048305967553619, 0.21267151063096745264, 0.18995868962243277774,
107 0.16512762256418831796, 0.13730498094001380420, 0.10483850756582017915, 0.063852163815003480173,
108 0
109};
110
111template<class RealType>
112const RealType exponential_table<RealType>::table_y[257] = {
113 0, 0.00045413435384149675545, 0.00096726928232717452884, 0.0015362997803015723824,
114 0.0021459677437189061793, 0.0027887987935740759640, 0.0034602647778369039855, 0.0041572951208337952532,
115 0.0048776559835423925804, 0.0056196422072054831710, 0.0063819059373191794422, 0.0071633531836349841425,
116 0.0079630774380170392396, 0.0087803149858089752347, 0.0096144136425022094101, 0.010464810181029979488,
117 0.011331013597834597488, 0.012212592426255380661, 0.013109164931254991070, 0.014020391403181937334,
118 0.014945968011691148079, 0.015885621839973162490, 0.016839106826039946359, 0.017806200410911360563,
119 0.018786700744696029497, 0.019780424338009741737, 0.020787204072578117603, 0.021806887504283582125,
120 0.022839335406385238829, 0.023884420511558170348, 0.024942026419731782971, 0.026012046645134218076,
121 0.027094383780955798424, 0.028188948763978634421, 0.029295660224637394015, 0.030414443910466605492,
122 0.031545232172893605499, 0.032687963508959533317, 0.033842582150874329031, 0.035009037697397411067,
123 0.036187284781931419754, 0.037377282772959360128, 0.038578995503074859626, 0.039792391023374122670,
124 0.041017441380414820816, 0.042254122413316231413, 0.043502413568888183301, 0.044762297732943280694,
125 0.046033761076175166762, 0.047316792913181548703, 0.048611385573379494401, 0.049917534282706374944,
126 0.051235237055126279830, 0.052564494593071689595, 0.053905310196046085104, 0.055257689676697038322,
127 0.056621641283742874438, 0.057997175631200659098, 0.059384305633420264487, 0.060783046445479636051,
128 0.062193415408540996150, 0.063615431999807331076, 0.065049117786753755036, 0.066494496385339779043,
129 0.067951593421936607770, 0.069420436498728751675, 0.070901055162371828426, 0.072393480875708743023,
130 0.073897746992364746308, 0.075413888734058408453, 0.076941943170480510100, 0.078481949201606426042,
131 0.080033947542319910023, 0.081597980709237420930, 0.083174093009632380354, 0.084762330532368125386,
132 0.086362741140756912277, 0.087975374467270219300, 0.089600281910032864534, 0.091237516631040162057,
133 0.092887133556043546523, 0.094549189376055853718, 0.096223742550432800103, 0.097910853311492199618,
134 0.099610583670637128826, 0.10132299742595363588, 0.10304816017125771553, 0.10478613930657016928,
135 0.10653700405000166218, 0.10830082545103379867, 0.11007767640518539026, 0.11186763167005629731,
136 0.11367076788274431301, 0.11548716357863353664, 0.11731689921155557057, 0.11916005717532768467,
137 0.12101672182667483729, 0.12288697950954513498, 0.12477091858083096578, 0.12666862943751066518,
138 0.12858020454522817870, 0.13050573846833078225, 0.13244532790138752023, 0.13439907170221363078,
139 0.13636707092642885841, 0.13834942886358021406, 0.14034625107486244210, 0.14235764543247220043,
140 0.14438372216063476473, 0.14642459387834493787, 0.14848037564386679222, 0.15055118500103990354,
141 0.15263714202744286154, 0.15473836938446807312, 0.15685499236936522013, 0.15898713896931420572,
142 0.16113493991759203183, 0.16329852875190180795, 0.16547804187493600915, 0.16767361861725019322,
143 0.16988540130252766513, 0.17211353531532005700, 0.17435816917135348788, 0.17661945459049489581,
144 0.17889754657247831241, 0.18119260347549629488, 0.18350478709776746150, 0.18583426276219711495,
145 0.18818119940425430485, 0.19054576966319540013, 0.19292814997677133873, 0.19532852067956322315,
146 0.19774706610509886464, 0.20018397469191127727, 0.20263943909370901930, 0.20511365629383770880,
147 0.20760682772422204205, 0.21011915938898825914, 0.21265086199297827522, 0.21520215107537867786,
148 0.21777324714870053264, 0.22036437584335949720, 0.22297576805812018050, 0.22560766011668406495,
149 0.22826029393071670664, 0.23093391716962742173, 0.23362878343743333945, 0.23634515245705964715,
150 0.23908329026244917002, 0.24184346939887722761, 0.24462596913189210901, 0.24743107566532763894,
151 0.25025908236886230967, 0.25311029001562948171, 0.25598500703041538015, 0.25888354974901621678,
152 0.26180624268936295243, 0.26475341883506220209, 0.26772541993204481808, 0.27072259679906003167,
153 0.27374530965280298302, 0.27679392844851734458, 0.27986883323697289920, 0.28297041453878076010,
154 0.28609907373707684673, 0.28925522348967773308, 0.29243928816189258772, 0.29565170428126120948,
155 0.29889292101558177099, 0.30216340067569352897, 0.30546361924459023541, 0.30879406693456016794,
156 0.31215524877417956945, 0.31554768522712893632, 0.31897191284495723773, 0.32242848495608914289,
157 0.32591797239355619822, 0.32944096426413633091, 0.33299806876180896713, 0.33658991402867758144,
158 0.34021714906678004560, 0.34388044470450243010, 0.34758049462163698567, 0.35131801643748334681,
159 0.35509375286678745925, 0.35890847294874976196, 0.36276297335481777335, 0.36665807978151414890,
160 0.37059464843514599421, 0.37457356761590215193, 0.37859575940958081092, 0.38266218149600982112,
161 0.38677382908413768115, 0.39093173698479710717, 0.39513698183329015336, 0.39939068447523107877,
162 0.40369401253053026739, 0.40804818315203238238, 0.41245446599716116772, 0.41691418643300289465,
163 0.42142872899761659635, 0.42599954114303435739, 0.43062813728845883923, 0.43531610321563659758,
164 0.44006510084235387501, 0.44487687341454851593, 0.44975325116275498919, 0.45469615747461548049,
165 0.45970761564213768669, 0.46478975625042618067, 0.46994482528395999841, 0.47517519303737738299,
166 0.48048336393045423016, 0.48587198734188493564, 0.49134386959403255500, 0.49690198724154955294,
167 0.50254950184134769289, 0.50828977641064283495, 0.51412639381474855788, 0.52006317736823356823,
168 0.52610421398361972602, 0.53225388026304326945, 0.53851687200286186590, 0.54489823767243963663,
169 0.55140341654064131685, 0.55803828226258748140, 0.56480919291240022434, 0.57172304866482579008,
170 0.57878735860284503057, 0.58601031847726802755, 0.59340090169173341521, 0.60096896636523224742,
171 0.60872538207962206507, 0.61668218091520762326, 0.62485273870366592605, 0.63325199421436607968,
172 0.64189671642726607018, 0.65080583341457104881, 0.66000084107899974178, 0.66950631673192477684,
173 0.67935057226476538741, 0.68956649611707798890, 0.70019265508278816709, 0.71127476080507597882,
174 0.72286765959357200702, 0.73503809243142351530, 0.74786862198519510742, 0.76146338884989624862,
175 0.77595685204011559675, 0.79152763697249565519, 0.80842165152300838005, 0.82699329664305033399,
176 0.84778550062398962096, 0.87170433238120363669, 0.90046992992574643800, 0.93814368086217467916,
177 1
178};
179
180template<class RealType = double>
181struct unit_exponential_distribution
182{
183 template<class Engine>
184 RealType operator()(Engine& eng) {
185 const double * const table_x = exponential_table<double>::table_x;
186 const double * const table_y = exponential_table<double>::table_y;
187 RealType shift(0);
188 for(;;) {
189 std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
190 int i = vals.second;
191 RealType x = vals.first * RealType(table_x[i]);
192 if(x < RealType(table_x[i + 1])) return shift + x;
193 // For i=0 we need to generate from the tail, but because this is an exponential
194 // distribution, the tail looks exactly like the body, so we can simply repeat with a
195 // shift:
196 if (i == 0) shift += RealType(table_x[1]);
197 else {
198 RealType y01 = uniform_01<RealType>()(eng);
199 RealType y = RealType(table_y[i]) + y01 * RealType(table_y[i+1] - table_y[i]);
200
201 // All we care about is whether these are < or > 0; these values are equal to
202 // (lbound) or proportional to (ubound) `y` minus the lower/upper bound.
203 RealType y_above_ubound = RealType(table_x[i] - table_x[i+1]) * y01 - (RealType(table_x[i]) - x),
204 y_above_lbound = y - (RealType(table_y[i+1]) + (RealType(table_x[i+1]) - x) * RealType(table_y[i+1]));
205
206 if (y_above_ubound < 0 // if above the upper bound reject immediately
207 &&
208 (
209 y_above_lbound < 0 // If below the lower bound accept immediately
210 ||
211 y < f(x) // Otherwise it's between the bounds and we need a full check
212 )
213 ) {
214 return x + shift;
215 }
216 }
217 }
218 }
219 static RealType f(RealType x) {
220 using std::exp;
221 return exp(-x);
222 }
223};
224
225} // namespace detail
226
227
228/**
229 * The exponential distribution is a model of \random_distribution with
230 * a single parameter lambda.
231 *
232 * It has \f$\displaystyle p(x) = \lambda e^{-\lambda x}\f$
233 *
234 * The implementation uses the "ziggurat" algorithm, as described in
235 *
236 * @blockquote
237 * "The Ziggurat Method for Generating Random Variables",
238 * George Marsaglia and Wai Wan Tsang, Journal of Statistical Software
239 * Volume 5, Number 8 (2000), 1-7.
240 * @endblockquote
241 */
242template<class RealType = double>
243class exponential_distribution
244{
245public:
246 typedef RealType input_type;
247 typedef RealType result_type;
248
249 class param_type
250 {
251 public:
252
253 typedef exponential_distribution distribution_type;
254
255 /**
256 * Constructs parameters with a given lambda.
257 *
258 * Requires: lambda > 0
259 */
260 param_type(RealType lambda_arg = RealType(1.0))
261 : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }
262
263 /** Returns the lambda parameter of the distribution. */
264 RealType lambda() const { return _lambda; }
265
266 /** Writes the parameters to a @c std::ostream. */
267 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
268 {
269 os << parm._lambda;
270 return os;
271 }
272
273 /** Reads the parameters from a @c std::istream. */
274 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
275 {
276 is >> parm._lambda;
277 return is;
278 }
279
280 /** Returns true if the two sets of parameters are equal. */
281 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
282 { return lhs._lambda == rhs._lambda; }
283
284 /** Returns true if the two sets of parameters are different. */
285 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
286
287 private:
288 RealType _lambda;
289 };
290
291 /**
292 * Constructs an exponential_distribution with a given lambda.
293 *
294 * Requires: lambda > 0
295 */
296 explicit exponential_distribution(RealType lambda_arg = RealType(1.0))
297 : _lambda(lambda_arg) { BOOST_ASSERT(_lambda > RealType(0)); }
298
299 /**
300 * Constructs an exponential_distribution from its parameters
301 */
302 explicit exponential_distribution(const param_type& parm)
303 : _lambda(parm.lambda()) {}
304
305 // compiler-generated copy ctor and assignment operator are fine
306
307 /** Returns the lambda parameter of the distribution. */
308 RealType lambda() const { return _lambda; }
309
310 /** Returns the smallest value that the distribution can produce. */
311 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
312 { return RealType(0); }
313 /** Returns the largest value that the distribution can produce. */
314 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
315 { return (std::numeric_limits<RealType>::infinity)(); }
316
317 /** Returns the parameters of the distribution. */
318 param_type param() const { return param_type(_lambda); }
319 /** Sets the parameters of the distribution. */
320 void param(const param_type& parm) { _lambda = parm.lambda(); }
321
322 /**
323 * Effects: Subsequent uses of the distribution do not depend
324 * on values produced by any engine prior to invoking reset.
325 */
326 void reset() { }
327
328 /**
329 * Returns a random variate distributed according to the
330 * exponential distribution.
331 */
332 template<class Engine>
333 result_type operator()(Engine& eng) const
334 {
335 detail::unit_exponential_distribution<RealType> impl;
336 return impl(eng) / _lambda;
337 }
338
339 /**
340 * Returns a random variate distributed according to the exponential
341 * distribution with parameters specified by param.
342 */
343 template<class Engine>
344 result_type operator()(Engine& eng, const param_type& parm) const
345 {
346 return exponential_distribution(parm)(eng);
347 }
348
349 /** Writes the distribution to a std::ostream. */
350 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, exponential_distribution, ed)
351 {
352 os << ed._lambda;
353 return os;
354 }
355
356 /** Reads the distribution from a std::istream. */
357 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, exponential_distribution, ed)
358 {
359 is >> ed._lambda;
360 return is;
361 }
362
363 /**
364 * Returns true iff the two distributions will produce identical
365 * sequences of values given equal generators.
366 */
367 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(exponential_distribution, lhs, rhs)
368 { return lhs._lambda == rhs._lambda; }
369
370 /**
371 * Returns true iff the two distributions will produce different
372 * sequences of values given equal generators.
373 */
374 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(exponential_distribution)
375
376private:
377 result_type _lambda;
378};
379
380} // namespace random
381
382using random::exponential_distribution;
383
384} // namespace boost
385
386#endif // BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
387