Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
Boost165NormalDistribution.hpp
Go to the documentation of this file.
1
24/* boost random/normal_distribution.hpp header file
25 *
26 * Copyright Jens Maurer 2000-2001
27 * Copyright Steven Watanabe 2010-2011
28 * Distributed under the Boost Software License, Version 1.0. (See
29 * accompanying file LICENSE_1_0.txt or copy at
30 * http://www.boost.org/LICENSE_1_0.txt)
31 *
32 * See http://www.boost.org for most recent version including documentation.
33 *
34 * $Id$
35 *
36 * Revision history
37 * 2001-02-18 moved to individual header files
38 */
39
40#ifndef BOOST_165_RANDOM_NORMAL_DISTRIBUTION_HPP
41#define BOOST_165_RANDOM_NORMAL_DISTRIBUTION_HPP
42
43#include <boost/assert.hpp>
44#include <boost/config/no_tr1/cmath.hpp>
45#include <boost/limits.hpp>
46#include <boost/random/detail/config.hpp>
47//#include <boost/random/detail/int_float_pair.hpp> // replaced with Chaste version as not present in old boosts < 1.64
48#include <boost/random/detail/operators.hpp>
49//#include <boost/random/exponential_distribution.hpp> // replaced with Chaste version as different in old boosts < 1.64
50#include <boost/random/uniform_01.hpp>
51#include <boost/static_assert.hpp>
52#include <iosfwd>
53#include <istream>
54
55// Added by Chaste team
58
59namespace boost
60{
61namespace random
62{
63 namespace detail
64 {
65
66 // tables for the ziggurat algorithm
67 template <class RealType>
69 {
70 static const RealType table_x[129];
71 static const RealType table_y[129];
72 };
73
74 template <class RealType>
75 const RealType normal_table_v165<RealType>::table_x[129] = {
76 3.7130862467403632609, 3.4426198558966521214, 3.2230849845786185446, 3.0832288582142137009,
77 2.9786962526450169606, 2.8943440070186706210, 2.8231253505459664379, 2.7611693723841538514,
78 2.7061135731187223371, 2.6564064112581924999, 2.6109722484286132035, 2.5690336259216391328,
79 2.5300096723854666170, 2.4934545220919507609, 2.4590181774083500943, 2.4264206455302115930,
80 2.3954342780074673425, 2.3658713701139875435, 2.3375752413355307354, 2.3104136836950021558,
81 2.2842740596736568056, 2.2590595738653295251, 2.2346863955870569803, 2.2110814088747278106,
82 2.1881804320720206093, 2.1659267937448407377, 2.1442701823562613518, 2.1231657086697899595,
83 2.1025731351849988838, 2.0824562379877246441, 2.0627822745039633575, 2.0435215366506694976,
84 2.0246469733729338782, 2.0061338699589668403, 1.9879595741230607243, 1.9701032608497132242,
85 1.9525457295488889058, 1.9352692282919002011, 1.9182573008597320303, 1.9014946531003176140,
86 1.8849670357028692380, 1.8686611409895420085, 1.8525645117230870617, 1.8366654602533840447,
87 1.8209529965910050740, 1.8054167642140487420, 1.7900469825946189862, 1.7748343955807692457,
88 1.7597702248942318749, 1.7448461281083765085, 1.7300541605582435350, 1.7153867407081165482,
89 1.7008366185643009437, 1.6863968467734863258, 1.6720607540918522072, 1.6578219209482075462,
90 1.6436741568569826489, 1.6296114794646783962, 1.6156280950371329644, 1.6017183802152770587,
91 1.5878768648844007019, 1.5740982160167497219, 1.5603772223598406870, 1.5467087798535034608,
92 1.5330878776675560787, 1.5195095847593707806, 1.5059690368565502602, 1.4924614237746154081,
93 1.4789819769830978546, 1.4655259573357946276, 1.4520886428822164926, 1.4386653166774613138,
94 1.4252512545068615734, 1.4118417124397602509, 1.3984319141236063517, 1.3850170377251486449,
95 1.3715922024197322698, 1.3581524543224228739, 1.3446927517457130432, 1.3312079496576765017,
96 1.3176927832013429910, 1.3041418501204215390, 1.2905495919178731508, 1.2769102735516997175,
97 1.2632179614460282310, 1.2494664995643337480, 1.2356494832544811749, 1.2217602305309625678,
98 1.2077917504067576028, 1.1937367078237721994, 1.1795873846544607035, 1.1653356361550469083,
99 1.1509728421389760651, 1.1364898520030755352, 1.1218769225722540661, 1.1071236475235353980,
100 1.0922188768965537614, 1.0771506248819376573, 1.0619059636836193998, 1.0464709007525802629,
101 1.0308302360564555907, 1.0149673952392994716, 0.99886423348064351303, 0.98250080350276038481,
102 0.96585507938813059489, 0.94890262549791195381, 0.93161619660135381056, 0.91396525100880177644,
103 0.89591535256623852894, 0.87742742909771569142, 0.85845684317805086354, 0.83895221428120745572,
104 0.81885390668331772331, 0.79809206062627480454, 0.77658398787614838598, 0.75423066443451007146,
105 0.73091191062188128150, 0.70647961131360803456, 0.68074791864590421664, 0.65347863871504238702,
106 0.62435859730908822111, 0.59296294244197797913, 0.55869217837551797140, 0.52065603872514491759,
107 0.47743783725378787681, 0.42654798630330512490, 0.36287143102841830424, 0.27232086470466385065,
108 0
109 };
110
111 template <class RealType>
112 const RealType normal_table_v165<RealType>::table_y[129] = {
113 0, 0.0026696290839025035092, 0.0055489952208164705392, 0.0086244844129304709682,
114 0.011839478657982313715, 0.015167298010672042468, 0.018592102737165812650, 0.022103304616111592615,
115 0.025693291936149616572, 0.029356317440253829618, 0.033087886146505155566, 0.036884388786968774128,
116 0.040742868074790604632, 0.044660862200872429800, 0.048636295860284051878, 0.052667401903503169793,
117 0.056752663481538584188, 0.060890770348566375972, 0.065080585213631873753, 0.069321117394180252601,
118 0.073611501884754893389, 0.077950982514654714188, 0.082338898242957408243, 0.086774671895542968998,
119 0.091257800827634710201, 0.09578784912257815216, 0.10036444102954554013, 0.10498725541035453978,
120 0.10965602101581776100, 0.11437051244988827452, 0.11913054670871858767, 0.12393598020398174246,
121 0.12878670619710396109, 0.13368265258464764118, 0.13862377998585103702, 0.14361008009193299469,
122 0.14864157424369696566, 0.15371831220958657066, 0.15884037114093507813, 0.16400785468492774791,
123 0.16922089223892475176, 0.17447963833240232295, 0.17978427212496211424, 0.18513499701071343216,
124 0.19053204032091372112, 0.19597565311811041399, 0.20146611007620324118, 0.20700370944187380064,
125 0.21258877307373610060, 0.21822164655637059599, 0.22390269938713388747, 0.22963232523430270355,
126 0.23541094226572765600, 0.24123899354775131610, 0.24711694751469673582, 0.25304529850976585934,
127 0.25902456739871074263, 0.26505530225816194029, 0.27113807914102527343, 0.27727350292189771153,
128 0.28346220822601251779, 0.28970486044581049771, 0.29600215684985583659, 0.30235482778947976274,
129 0.30876363800925192282, 0.31522938806815752222, 0.32175291587920862031, 0.32833509837615239609,
130 0.33497685331697116147, 0.34167914123501368412, 0.34844296754987246935, 0.35526938485154714435,
131 0.36215949537303321162, 0.36911445366827513952, 0.37613546951445442947, 0.38322381105988364587,
132 0.39038080824138948916, 0.39760785649804255208, 0.40490642081148835099, 0.41227804010702462062,
133 0.41972433205403823467, 0.42724699830956239880, 0.43484783025466189638, 0.44252871528024661483,
134 0.45029164368692696086, 0.45813871627287196483, 0.46607215269457097924, 0.47409430069824960453,
135 0.48220764633483869062, 0.49041482528932163741, 0.49871863547658432422, 0.50712205108130458951,
136 0.51562823824987205196, 0.52424057267899279809, 0.53296265938998758838, 0.54179835503172412311,
137 0.55075179312105527738, 0.55982741271069481791, 0.56902999107472161225, 0.57836468112670231279,
138 0.58783705444182052571, 0.59745315095181228217, 0.60721953663260488551, 0.61714337082656248870,
139 0.62723248525781456578, 0.63749547734314487428, 0.64794182111855080873, 0.65858200005865368016,
140 0.66942766735770616891, 0.68049184100641433355, 0.69178914344603585279, 0.70333609902581741633,
141 0.71515150742047704368, 0.72725691835450587793, 0.73967724368333814856, 0.75244155918570380145,
142 0.76558417390923599480, 0.77914608594170316563, 0.79317701178385921053, 0.80773829469612111340,
143 0.82290721139526200050, 0.83878360531064722379, 0.85550060788506428418, 0.87324304892685358879,
144 0.89228165080230272301, 0.91304364799203805999, 0.93628268170837107547, 0.96359969315576759960,
145 1
146 };
147
148 template <class RealType = double>
150 {
151 template <class Engine>
152 RealType operator()(Engine& eng)
153 {
154 const double* const table_x = normal_table_v165<double>::table_x;
155 const double* const table_y = normal_table_v165<double>::table_y;
156 for (;;)
157 {
158 std::pair<RealType, int> vals = generate_int_float_pair_v165<RealType, 8>(eng);
159 int i = vals.second;
160 int sign = (i & 1) * 2 - 1;
161 i = i >> 1;
162 RealType x = vals.first * RealType(table_x[i]);
163 if (x < table_x[i + 1])
164 return x * sign;
165 if (i == 0)
166 return generate_tail_v165(eng) * sign;
167
168 RealType y01 = uniform_01<RealType>()(eng);
169 RealType y = RealType(table_y[i]) + y01 * RealType(table_y[i + 1] - table_y[i]);
170
171 // These store the value y - bound, or something proportional to that difference:
172 RealType y_above_ubound, y_above_lbound;
173
174 // There are three cases to consider:
175 // - convex regions (where x[i] > x[j] >= 1)
176 // - concave regions (where 1 <= x[i] < x[j])
177 // - region containing the inflection point (where x[i] > 1 > x[j])
178 // For convex (concave), exp^(-x^2/2) is bounded below (above) by the tangent at
179 // (x[i],y[i]) and is bounded above (below) by the diagonal line from (x[i+1],y[i+1]) to
180 // (x[i],y[i]).
181 //
182 // *If* the inflection point region satisfies slope(x[i+1]) < slope(diagonal), then we
183 // can treat the inflection region as a convex region: this condition is necessary and
184 // sufficient to ensure that the curve lies entirely below the diagonal (that, in turn,
185 // also implies that it will be above the tangent at x[i]).
186 //
187 // For the current table size (128), this is satisfied: slope(x[i+1]) = -0.60653 <
188 // slope(diag) = -0.60649, and so we have only two cases below instead of three.
189
190 if (table_x[i] >= 1)
191 { // convex (incl. inflection)
192 y_above_ubound = RealType(table_x[i] - table_x[i + 1]) * y01 - (RealType(table_x[i]) - x);
193 y_above_lbound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
194 }
195 else
196 { // concave
197 y_above_lbound = RealType(table_x[i] - table_x[i + 1]) * y01 - (RealType(table_x[i]) - x);
198 y_above_ubound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
199 }
200
201 if (y_above_ubound < 0 // if above the upper bound reject immediately
202 && (y_above_lbound < 0 // If below the lower bound accept immediately
203 || y < f(x) // Otherwise it's between the bounds and we need a full check
204 ))
205 {
206 return x * sign;
207 }
208 }
209 }
210 static RealType f(RealType x)
211 {
212 using std::exp;
213 return exp(-(x * x / 2));
214 }
215 // Generate from the tail using rejection sampling from the exponential(x_1) distribution,
216 // shifted by x_1. This looks a little different from the usual rejection sampling because it
217 // transforms the condition by taking the log of both sides, thus avoiding the costly exp() call
218 // on the RHS, then takes advantage of the fact that -log(unif01) is simply generating an
219 // exponential (by inverse cdf sampling) by replacing the log(unif01) on the LHS with a
220 // exponential(1) draw, y.
221 template <class Engine>
222 RealType generate_tail_v165(Engine& eng)
223 {
224 const RealType tail_start = RealType(normal_table_v165<double>::table_x[1]);
227 for (;;)
228 {
229 RealType x = exp_x(eng);
230 RealType y = exp_y(eng);
231 // If we were doing non-transformed rejection sampling, this condition would be:
232 // if (unif01 < exp(-.5*x*x)) return x + tail_start;
233 if (2 * y > x * x)
234 return x + tail_start;
235 }
236 }
237 };
238
239 } // namespace detail
240
258 template <class RealType = double>
260 {
261 public:
262 typedef RealType input_type;
263 typedef RealType result_type;
264
266 {
267 public:
269
276 explicit param_type(RealType mean_arg = RealType(0.0),
277 RealType sigma_arg = RealType(1.0))
278 : _mean(mean_arg),
279 _sigma(sigma_arg)
280 {
281 }
282
284 RealType mean() const { return _mean; }
285
287 RealType sigma() const { return _sigma; }
288
291 {
292 os << parm._mean << " " << parm._sigma;
293 return os;
294 }
295
298 {
299 is >> parm._mean >> std::ws >> parm._sigma;
300 return is;
301 }
302
305 {
306 return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
307 }
308
310 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
311
312 private:
313 RealType _mean;
314 RealType _sigma;
315 };
316
323 explicit normal_distribution_v165(const RealType& mean_arg = RealType(0.0),
324 const RealType& sigma_arg = RealType(1.0))
325 : _mean(mean_arg), _sigma(sigma_arg)
326 {
327 BOOST_ASSERT(_sigma >= RealType(0));
328 }
329
334 : _mean(parm.mean()), _sigma(parm.sigma())
335 {
336 }
337
339 RealType mean() const { return _mean; }
341 RealType sigma() const { return _sigma; }
342
345 {
346 return -std::numeric_limits<RealType>::infinity();
347 }
350 {
351 return std::numeric_limits<RealType>::infinity();
352 }
353
355 param_type param() const { return param_type(_mean, _sigma); }
357 void param(const param_type& parm)
358 {
359 _mean = parm.mean();
360 _sigma = parm.sigma();
361 }
362
367 void reset() {}
368
370 template <class Engine>
371 result_type operator()(Engine& eng)
372 {
374 return impl(eng) * _sigma + _mean;
375 }
376
378 template <class URNG>
379 result_type operator()(URNG& urng, const param_type& parm)
380 {
381 return normal_distribution_v165(parm)(urng);
382 }
383
386 {
387 os << nd._mean << " " << nd._sigma;
388 return os;
389 }
390
393 {
394 is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
395 return is;
396 }
397
403 {
404 return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
405 }
406
411 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(normal_distribution_v165)
412
413 private:
414 RealType _mean, _sigma;
415 };
416} // namespace random
417} // namespace boost
418
419#endif // BOOST_165_RANDOM_NORMAL_DISTRIBUTION_HPP
param_type(RealType mean_arg=RealType(0.0), RealType sigma_arg=RealType(1.0))
result_type operator()(URNG &urng, const param_type &parm)
normal_distribution_v165(const RealType &mean_arg=RealType(0.0), const RealType &sigma_arg=RealType(1.0))
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(normal_distribution_v165, lhs, rhs)
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution_v165, nd)
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution_v165, nd)