Line data Source code
1 : /* Copyright (c) 2018-2021, The Tor Project, Inc. */
2 : /* See LICENSE for licensing information */
3 :
4 : /**
5 : * \file prob_distr.c
6 : *
7 : * \brief
8 : * Implements various probability distributions.
9 : * Almost all code is courtesy of Riastradh.
10 : *
11 : * \details
12 : * Here are some details that might help you understand this file:
13 : *
14 : * - Throughout this file, `eps' means the largest relative error of a
15 : * correctly rounded floating-point operation, which in binary64
16 : * floating-point arithmetic is 2^-53. Here the relative error of a
17 : * true value x from a computed value y is |x - y|/|x|. This
18 : * definition of epsilon is conventional for numerical analysts when
19 : * writing error analyses. (If your libm doesn't provide correctly
20 : * rounded exp and log, their relative error is usually below 2*2^-53
21 : * and probably closer to 1.1*2^-53 instead.)
22 : *
23 : * The C constant DBL_EPSILON is actually twice this, and should
24 : * perhaps rather be named ulp(1) -- that is, it is the distance from
25 : * 1 to the next greater floating-point number, which is usually of
26 : * more interest to programmers and hardware engineers.
27 : *
28 : * Since this file is concerned mainly with error bounds rather than
29 : * with low-level bit-hacking of floating-point numbers, we adopt the
30 : * numerical analysts' definition in the comments, though we do use
31 : * DBL_EPSILON in a handful of places where it is convenient to use
32 : * some function of eps = DBL_EPSILON/2 in a case analysis.
33 : *
34 : * - In various functions (e.g. sample_log_logistic()) we jump through hoops so
35 : * that we can use reals closer to 0 than closer to 1, since we achieve much
36 : * greater accuracy for floating point numbers near 0. In particular, we can
37 : * represent differences as small as 10^-300 for numbers near 0, but of no
38 : * less than 10^-16 for numbers near 1.
39 : **/
40 :
41 : #define PROB_DISTR_PRIVATE
42 :
43 : #include "orconfig.h"
44 :
45 : #include "lib/math/prob_distr.h"
46 :
47 : #include "lib/crypt_ops/crypto_rand.h"
48 : #include "lib/cc/ctassert.h"
49 : #include "lib/log/util_bug.h"
50 :
51 : #include <float.h>
52 : #include <math.h>
53 : #include <stddef.h>
54 :
55 : #ifndef COCCI
56 : /** Declare a function that downcasts from a generic dist struct to the actual
57 : * subtype probablity distribution it represents. */
58 : #define DECLARE_PROB_DISTR_DOWNCAST_FN(name) \
59 : static inline \
60 : const struct name##_t * \
61 : dist_to_const_##name(const struct dist_t *obj) { \
62 : tor_assert(obj->ops == &name##_ops); \
63 : return SUBTYPE_P(obj, struct name ## _t, base); \
64 : }
65 601320 : DECLARE_PROB_DISTR_DOWNCAST_FN(uniform)
66 400100 : DECLARE_PROB_DISTR_DOWNCAST_FN(geometric)
67 400904 : DECLARE_PROB_DISTR_DOWNCAST_FN(logistic)
68 400904 : DECLARE_PROB_DISTR_DOWNCAST_FN(log_logistic)
69 701507 : DECLARE_PROB_DISTR_DOWNCAST_FN(genpareto)
70 501105 : DECLARE_PROB_DISTR_DOWNCAST_FN(weibull)
71 : #endif /* !defined(COCCI) */
72 :
73 : /**
74 : * Count number of one bits in 32-bit word.
75 : */
76 : static unsigned
77 3400714 : bitcount32(uint32_t x)
78 : {
79 :
80 : /* Count two-bit groups. */
81 3400714 : x -= (x >> 1) & UINT32_C(0x55555555);
82 :
83 : /* Count four-bit groups. */
84 3400714 : x = ((x >> 2) & UINT32_C(0x33333333)) + (x & UINT32_C(0x33333333));
85 :
86 : /* Count eight-bit groups. */
87 3400714 : x = (x + (x >> 4)) & UINT32_C(0x0f0f0f0f);
88 :
89 : /* Sum all eight-bit groups, and extract the sum. */
90 3400714 : return (x * UINT32_C(0x01010101)) >> 24;
91 : }
92 :
93 : /**
94 : * Count leading zeros in 32-bit word.
95 : */
96 : static unsigned
97 3400714 : clz32(uint32_t x)
98 : {
99 :
100 : /* Round up to a power of two. */
101 3400714 : x |= x >> 1;
102 3400714 : x |= x >> 2;
103 3400714 : x |= x >> 4;
104 3400714 : x |= x >> 8;
105 3400714 : x |= x >> 16;
106 :
107 : /* Subtract count of one bits from 32. */
108 3400714 : return (32 - bitcount32(x));
109 : }
110 :
111 : /*
112 : * Some lemmas that will be used throughout this file to prove various error
113 : * bounds:
114 : *
115 : * Lemma 1. If |d| <= 1/2, then 1/(1 + d) <= 2.
116 : *
117 : * Proof. If 0 <= d <= 1/2, then 1 + d >= 1, so that 1/(1 + d) <= 1.
118 : * If -1/2 <= d <= 0, then 1 + d >= 1/2, so that 1/(1 + d) <= 2. QED.
119 : *
120 : * Lemma 2. If b = a*(1 + d)/(1 + d') for |d'| < 1/2 and nonzero a, b,
121 : * then b = a*(1 + e) for |e| <= 2|d' - d|.
122 : *
123 : * Proof. |a - b|/|a|
124 : * = |a - a*(1 + d)/(1 + d')|/|a|
125 : * = |1 - (1 + d)/(1 + d')|
126 : * = |(1 + d' - 1 - d)/(1 + d')|
127 : * = |(d' - d)/(1 + d')|
128 : * <= 2|d' - d|, by Lemma 1,
129 : *
130 : * QED.
131 : *
132 : * Lemma 3. For |d|, |d'| < 1/4,
133 : *
134 : * |log((1 + d)/(1 + d'))| <= 4|d - d'|.
135 : *
136 : * Proof. Write
137 : *
138 : * log((1 + d)/(1 + d'))
139 : * = log(1 + (1 + d)/(1 + d') - 1)
140 : * = log(1 + (1 + d - 1 - d')/(1 + d')
141 : * = log(1 + (d - d')/(1 + d')).
142 : *
143 : * By Lemma 1, |(d - d')/(1 + d')| < 2|d' - d| < 1, so the Taylor
144 : * series of log(1 + x) converges absolutely for (d - d')/(1 + d'),
145 : * and thus we have
146 : *
147 : * |log(1 + (d - d')/(1 + d'))|
148 : * = |\sum_{n=1}^\infty ((d - d')/(1 + d'))^n/n|
149 : * <= \sum_{n=1}^\infty |(d - d')/(1 + d')|^n/n
150 : * <= \sum_{n=1}^\infty |2(d' - d)|^n/n
151 : * <= \sum_{n=1}^\infty |2(d' - d)|^n
152 : * = 1/(1 - |2(d' - d)|)
153 : * <= 4|d' - d|,
154 : *
155 : * QED.
156 : *
157 : * Lemma 4. If 1/e <= 1 + x <= e, then
158 : *
159 : * log(1 + (1 + d) x) = (1 + d') log(1 + x)
160 : *
161 : * for |d'| < 8|d|.
162 : *
163 : * Proof. Write
164 : *
165 : * log(1 + (1 + d) x)
166 : * = log(1 + x + x*d)
167 : * = log((1 + x) (1 + x + x*d)/(1 + x))
168 : * = log(1 + x) + log((1 + x + x*d)/(1 + x))
169 : * = log(1 + x) (1 + log((1 + x + x*d)/(1 + x))/log(1 + x)).
170 : *
171 : * The relative error is bounded by
172 : *
173 : * |log((1 + x + x*d)/(1 + x))/log(1 + x)|
174 : * <= 4|x + x*d - x|/|log(1 + x)|, by Lemma 3,
175 : * = 4|x*d|/|log(1 + x)|
176 : * < 8|d|,
177 : *
178 : * since in this range 0 < 1 - 1/e < x/log(1 + x) <= e - 1 < 2. QED.
179 : */
180 :
181 : /**
182 : * Compute the logistic function: f(x) = 1/(1 + e^{-x}) = e^x/(1 + e^x).
183 : * Maps a log-odds-space probability in [-infinity, +infinity] into a
184 : * direct-space probability in [0,1]. Inverse of logit.
185 : *
186 : * Ill-conditioned for large x; the identity logistic(-x) = 1 -
187 : * logistic(x) and the function logistichalf(x) = logistic(x) - 1/2 may
188 : * help to rearrange a computation.
189 : *
190 : * This implementation gives relative error bounded by 7 eps.
191 : */
192 : STATIC double
193 1833 : logistic(double x)
194 : {
195 1833 : if (x <= log(DBL_EPSILON/2)) {
196 : /*
197 : * If x <= log(DBL_EPSILON/2) = log(eps), then e^x <= eps. In this case
198 : * we will approximate the logistic() function with e^x because the
199 : * relative error is less than eps. Here is a calculation of the
200 : * relative error between the logistic() function and e^x and a proof
201 : * that it's less than eps:
202 : *
203 : * |e^x - e^x/(1 + e^x)|/|e^x/(1 + e^x)|
204 : * <= |1 - 1/(1 + e^x)|*|1 + e^x|
205 : * = |e^x/(1 + e^x)|*|1 + e^x|
206 : * = |e^x|
207 : * <= eps.
208 : */
209 141 : return exp(x); /* return e^x */
210 1692 : } else if (x <= -log(DBL_EPSILON/2)) {
211 : /*
212 : * e^{-x} > 0, so 1 + e^{-x} > 1, and 0 < 1/(1 +
213 : * e^{-x}) < 1; further, since e^{-x} < 1 + e^{-x}, we
214 : * also have 0 < 1/(1 + e^{-x}) < 1. Thus, if exp has
215 : * relative error d0, + has relative error d1, and /
216 : * has relative error d2, then we get
217 : *
218 : * (1 + d2)/[(1 + (1 + d0) e^{-x})(1 + d1)]
219 : * = (1 + d0)/[1 + e^{-x} + d0 e^{-x}
220 : * + d1 + d1 e^{-x} + d0 d1 e^{-x}]
221 : * = (1 + d0)/[(1 + e^{-x})
222 : * * (1 + d0 e^{-x}/(1 + e^{-x})
223 : * + d1/(1 + e^{-x})
224 : * + d0 d1 e^{-x}/(1 + e^{-x}))].
225 : * = (1 + d0)/[(1 + e^{-x})(1 + d')]
226 : * = [1/(1 + e^{-x})] (1 + d0)/(1 + d')
227 : *
228 : * where
229 : *
230 : * d' = d0 e^{-x}/(1 + e^{-x})
231 : * + d1/(1 + e^{-x})
232 : * + d0 d1 e^{-x}/(1 + e^{-x}).
233 : *
234 : * By Lemma 2 this relative error is bounded by
235 : *
236 : * 2|d0 - d'|
237 : * = 2|d0 - d0 e^{-x}/(1 + e^{-x})
238 : * - d1/(1 + e^{-x})
239 : * - d0 d1 e^{-x}/(1 + e^{-x})|
240 : * <= 2|d0| + 2|d0 e^{-x}/(1 + e^{-x})|
241 : * + 2|d1/(1 + e^{-x})|
242 : * + 2|d0 d1 e^{-x}/(1 + e^{-x})|
243 : * <= 2|d0| + 2|d0| + 2|d1| + 2|d0 d1|
244 : * <= 4|d0| + 2|d1| + 2|d0 d1|
245 : * <= 6 eps + 2 eps^2.
246 : */
247 1548 : return 1/(1 + exp(-x));
248 : } else {
249 : /*
250 : * e^{-x} <= eps, so the relative error of 1 from 1/(1
251 : * + e^{-x}) is
252 : *
253 : * |1/(1 + e^{-x}) - 1|/|1/(1 + e^{-x})|
254 : * = |e^{-x}/(1 + e^{-x})|/|1/(1 + e^{-x})|
255 : * = |e^{-x}|
256 : * <= eps.
257 : *
258 : * This computation avoids an intermediate overflow
259 : * exception, although the effect on the result is
260 : * harmless.
261 : *
262 : * XXX Should maybe raise inexact here.
263 : */
264 : return 1;
265 : }
266 : }
267 :
268 : /**
269 : * Compute the logit function: log p/(1 - p). Defined on [0,1]. Maps
270 : * a direct-space probability in [0,1] to a log-odds-space probability
271 : * in [-infinity, +infinity]. Inverse of logistic.
272 : *
273 : * Ill-conditioned near 1/2 and 1; the identity logit(1 - p) =
274 : * -logit(p) and the function logithalf(p0) = logit(1/2 + p0) may help
275 : * to rearrange a computation for p in [1/(1 + e), 1 - 1/(1 + e)].
276 : *
277 : * This implementation gives relative error bounded by 10 eps.
278 : */
279 : STATIC double
280 216482 : logit(double p)
281 : {
282 :
283 : /* logistic(-1) <= p <= logistic(+1) */
284 216482 : if (1/(1 + exp(1)) <= p && p <= 1/(1 + exp(-1))) {
285 : /*
286 : * For inputs near 1/2, we want to compute log1p(near
287 : * 0) rather than log(near 1), so write this as:
288 : *
289 : * log(p/(1 - p)) = -log((1 - p)/p)
290 : * = -log(1 + (1 - p)/p - 1)
291 : * = -log(1 + (1 - p - p)/p)
292 : * = -log(1 + (1 - 2p)/p).
293 : *
294 : * Since p = 2p/2 <= 1 <= 2*2p = 4p, the floating-point
295 : * evaluation of 1 - 2p is exact; the only error arises
296 : * from division and log1p. First, note that if
297 : * logistic(-1) <= p <= logistic(+1), (1 - 2p)/p lies
298 : * in the bounds of Lemma 4.
299 : *
300 : * If division has relative error d0 and log1p has
301 : * relative error d1, the outcome is
302 : *
303 : * -(1 + d1) log(1 + (1 - 2p) (1 + d0)/p)
304 : * = -(1 + d1) (1 + d') log(1 + (1 - 2p)/p)
305 : * = -(1 + d1 + d' + d1 d') log(1 + (1 - 2p)/p).
306 : *
307 : * where |d'| < 8|d0| by Lemma 4. The relative error
308 : * is then bounded by
309 : *
310 : * |d1 + d' + d1 d'|
311 : * <= |d1| + 8|d0| + 8|d1 d0|
312 : * <= 9 eps + 8 eps^2.
313 : */
314 446 : return -log1p((1 - 2*p)/p);
315 : } else {
316 : /*
317 : * For inputs near 0, although 1 - p may be rounded to
318 : * 1, it doesn't matter much because the magnitude of
319 : * the result is so much larger. For inputs near 1, we
320 : * can compute 1 - p exactly, although the precision on
321 : * the input is limited so we won't ever get more than
322 : * about 700 for the output.
323 : *
324 : * If - has relative error d0, / has relative error d1,
325 : * and log has relative error d2, then
326 : *
327 : * (1 + d2) log((1 + d0) p/[(1 - p)(1 + d1)])
328 : * = (1 + d2) [log(p/(1 - p)) + log((1 + d0)/(1 + d1))]
329 : * = log(p/(1 - p)) + d2 log(p/(1 - p))
330 : * + (1 + d2) log((1 + d0)/(1 + d1))
331 : * = log(p/(1 - p))*[1 + d2 +
332 : * + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))]
333 : *
334 : * Since 0 <= p < logistic(-1) or logistic(+1) < p <=
335 : * 1, we have |log(p/(1 - p))| > 1. Hence this error
336 : * is bounded by
337 : *
338 : * |d2 + (1 + d2) log((1 + d0)/(1 + d1))/log(p/(1 - p))|
339 : * <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))
340 : * / log(p/(1 - p))|
341 : * <= |d2| + |(1 + d2) log((1 + d0)/(1 + d1))|
342 : * <= |d2| + 4|(1 + d2) (d0 - d1)|, by Lemma 3,
343 : * <= |d2| + 4|d0 - d1 + d2 d0 - d1 d0|
344 : * <= |d2| + 4|d0| + 4|d1| + 4|d2 d0| + 4|d1 d0|
345 : * <= 9 eps + 8 eps^2.
346 : */
347 216036 : return log(p/(1 - p));
348 : }
349 : }
350 :
351 : /**
352 : * Compute the logit function, translated in input by 1/2: logithalf(p)
353 : * = logit(1/2 + p). Defined on [-1/2, 1/2]. Inverse of logistichalf.
354 : *
355 : * Ill-conditioned near +/-1/2. If |p0| > 1/2 - 1/(1 + e), it may be
356 : * better to compute 1/2 + p0 or -1/2 - p0 and to use logit instead.
357 : * This implementation gives relative error bounded by 34 eps.
358 : */
359 : STATIC double
360 186082 : logithalf(double p0)
361 : {
362 :
363 186082 : if (fabs(p0) <= 0.5 - 1/(1 + exp(1))) {
364 : /*
365 : * logit(1/2 + p0)
366 : * = log((1/2 + p0)/(1 - (1/2 + p0)))
367 : * = log((1/2 + p0)/(1/2 - p0))
368 : * = log(1 + (1/2 + p0)/(1/2 - p0) - 1)
369 : * = log(1 + (1/2 + p0 - (1/2 - p0))/(1/2 - p0))
370 : * = log(1 + (1/2 + p0 - 1/2 + p0)/(1/2 - p0))
371 : * = log(1 + 2 p0/(1/2 - p0))
372 : *
373 : * If the error of subtraction is d0, the error of
374 : * division is d1, and the error of log1p is d2, then
375 : * what we compute is
376 : *
377 : * (1 + d2) log(1 + (1 + d1) 2 p0/[(1 + d0) (1/2 - p0)])
378 : * = (1 + d2) log(1 + (1 + d') 2 p0/(1/2 - p0))
379 : * = (1 + d2) (1 + d'') log(1 + 2 p0/(1/2 - p0))
380 : * = (1 + d2 + d'' + d2 d'') log(1 + 2 p0/(1/2 - p0)),
381 : *
382 : * where |d'| < 2|d0 - d1| <= 4 eps by Lemma 2, and
383 : * |d''| < 8|d'| < 32 eps by Lemma 4 since
384 : *
385 : * 1/e <= 1 + 2*p0/(1/2 - p0) <= e
386 : *
387 : * when |p0| <= 1/2 - 1/(1 + e). Hence the relative
388 : * error is bounded by
389 : *
390 : * |d2 + d'' + d2 d''|
391 : * <= |d2| + |d''| + |d2 d''|
392 : * <= |d1| + 32 |d0| + 32 |d1 d0|
393 : * <= 33 eps + 32 eps^2.
394 : */
395 186062 : return log1p(2*p0/(0.5 - p0));
396 : } else {
397 : /*
398 : * We have a choice of computing logit(1/2 + p0) or
399 : * -logit(1 - (1/2 + p0)) = -logit(1/2 - p0). It
400 : * doesn't matter which way we do this: either way,
401 : * since 1/2 p0 <= 1/2 <= 2 p0, the sum and difference
402 : * are computed exactly. So let's do the one that
403 : * skips the final negation.
404 : *
405 : * The result is
406 : *
407 : * (1 + d1) log((1 + d0) (1/2 + p0)/[(1 + d2) (1/2 - p0)])
408 : * = (1 + d1) (1 + log((1 + d0)/(1 + d2))
409 : * / log((1/2 + p0)/(1/2 - p0)))
410 : * * log((1/2 + p0)/(1/2 - p0))
411 : * = (1 + d') log((1/2 + p0)/(1/2 - p0))
412 : * = (1 + d') logit(1/2 + p0)
413 : *
414 : * where
415 : *
416 : * d' = d1 + log((1 + d0)/(1 + d2))/logit(1/2 + p0)
417 : * + d1 log((1 + d0)/(1 + d2))/logit(1/2 + p0).
418 : *
419 : * For |p| > 1/2 - 1/(1 + e), logit(1/2 + p0) > 1.
420 : * Provided |d0|, |d2| < 1/4, by Lemma 3 we have
421 : *
422 : * |log((1 + d0)/(1 + d2))| <= 4|d0 - d2|.
423 : *
424 : * Hence the relative error is bounded by
425 : *
426 : * |d'| <= |d1| + 4|d0 - d2| + 4|d1| |d0 - d2|
427 : * <= |d1| + 4|d0| + 4|d2| + 4|d1 d0| + 4|d1 d2|
428 : * <= 9 eps + 8 eps^2.
429 : */
430 20 : return log((0.5 + p0)/(0.5 - p0));
431 : }
432 : }
433 :
434 : /*
435 : * The following random_uniform_01 is tailored for IEEE 754 binary64
436 : * floating-point or smaller. It can be adapted to larger
437 : * floating-point formats like i387 80-bit or IEEE 754 binary128, but
438 : * it may require sampling more bits.
439 : */
440 : CTASSERT(FLT_RADIX == 2);
441 : CTASSERT(-DBL_MIN_EXP <= 1021);
442 : CTASSERT(DBL_MANT_DIG <= 53);
443 :
444 : /**
445 : * Draw a floating-point number in [0, 1] with uniform distribution.
446 : *
447 : * Note that the probability of returning 0 is less than 2^-1074, so
448 : * callers need not check for it. However, callers that cannot handle
449 : * rounding to 1 must deal with that, because it occurs with
450 : * probability 2^-54, which is small but nonnegligible.
451 : */
452 : STATIC double
453 3400714 : random_uniform_01(void)
454 : {
455 3400714 : uint32_t z, x, hi, lo;
456 3400714 : double s;
457 :
458 : /*
459 : * Draw an exponent, geometrically distributed, but give up if
460 : * we get a run of more than 1088 zeros, which really means the
461 : * system is broken.
462 : */
463 3400714 : z = 0;
464 3400714 : while ((x = crypto_fast_rng_get_u32(get_thread_fast_rng())) == 0) {
465 0 : if (z >= 1088)
466 : /* Your bit sampler is broken. Go home. */
467 : return 0;
468 0 : z += 32;
469 : }
470 3400714 : z += clz32(x);
471 :
472 : /*
473 : * Pick 32-bit halves of an odd normalized significand.
474 : * Picking it odd breaks ties in the subsequent rounding, which
475 : * occur only with measure zero in the uniform distribution on
476 : * [0, 1].
477 : */
478 3400714 : hi = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x80000000);
479 3400714 : lo = crypto_fast_rng_get_u32(get_thread_fast_rng()) | UINT32_C(0x00000001);
480 :
481 : /* Round to nearest scaled significand in [2^63, 2^64]. */
482 3400714 : s = hi*(double)4294967296 + lo;
483 :
484 : /* Rescale into [1/2, 1] and apply exponent in one swell foop. */
485 3400714 : return s * ldexp(1, -(64 + z));
486 : }
487 :
488 : /*******************************************************************/
489 :
490 : /* Functions for specific probability distributions start here: */
491 :
492 : /*
493 : * Logistic(mu, sigma) distribution, supported on (-infinity,+infinity)
494 : *
495 : * This is the uniform distribution on [0,1] mapped into log-odds
496 : * space, scaled by sigma and translated by mu.
497 : *
498 : * pdf(x) = e^{-(x - mu)/sigma} sigma (1 + e^{-(x - mu)/sigma})^2
499 : * cdf(x) = 1/(1 + e^{-(x - mu)/sigma}) = logistic((x - mu)/sigma)
500 : * sf(x) = 1 - cdf(x) = 1 - logistic((x - mu)/sigma = logistic(-(x - mu)/sigma)
501 : * icdf(p) = mu + sigma log p/(1 - p) = mu + sigma logit(p)
502 : * isf(p) = mu + sigma log (1 - p)/p = mu - sigma logit(p)
503 : */
504 :
505 : /**
506 : * Compute the CDF of the Logistic(mu, sigma) distribution: the
507 : * logistic function. Well-conditioned for negative inputs and small
508 : * positive inputs; ill-conditioned for large positive inputs.
509 : */
510 : STATIC double
511 794 : cdf_logistic(double x, double mu, double sigma)
512 : {
513 794 : return logistic((x - mu)/sigma);
514 : }
515 :
516 : /**
517 : * Compute the SF of the Logistic(mu, sigma) distribution: the logistic
518 : * function reflected over the y axis. Well-conditioned for positive
519 : * inputs and small negative inputs; ill-conditioned for large negative
520 : * inputs.
521 : */
522 : STATIC double
523 810 : sf_logistic(double x, double mu, double sigma)
524 : {
525 810 : return logistic(-(x - mu)/sigma);
526 : }
527 :
528 : /**
529 : * Compute the inverse of the CDF of the Logistic(mu, sigma)
530 : * distribution: the logit function. Well-conditioned near 0;
531 : * ill-conditioned near 1/2 and 1.
532 : */
533 : STATIC double
534 314 : icdf_logistic(double p, double mu, double sigma)
535 : {
536 314 : return mu + sigma*logit(p);
537 : }
538 :
539 : /**
540 : * Compute the inverse of the SF of the Logistic(mu, sigma)
541 : * distribution: the -logit function. Well-conditioned near 0;
542 : * ill-conditioned near 1/2 and 1.
543 : */
544 : STATIC double
545 310 : isf_logistic(double p, double mu, double sigma)
546 : {
547 310 : return mu - sigma*logit(p);
548 : }
549 :
550 : /*
551 : * LogLogistic(alpha, beta) distribution, supported on (0, +infinity).
552 : *
553 : * This is the uniform distribution on [0,1] mapped into odds space,
554 : * scaled by positive alpha and shaped by positive beta.
555 : *
556 : * Equivalent to computing exp of a Logistic(log alpha, 1/beta) sample.
557 : * (Name arises because the pdf has LogLogistic(x; alpha, beta) =
558 : * Logistic(log x; log alpha, 1/beta) and mathematicians got their
559 : * covariance contravariant.)
560 : *
561 : * pdf(x) = (beta/alpha) (x/alpha)^{beta - 1}/(1 + (x/alpha)^beta)^2
562 : * = (1/e^mu sigma) (x/e^mu)^{1/sigma - 1} /
563 : * (1 + (x/e^mu)^{1/sigma})^2
564 : * cdf(x) = 1/(1 + (x/alpha)^-beta) = 1/(1 + (x/e^mu)^{-1/sigma})
565 : * = 1/(1 + (e^{log x}/e^mu)^{-1/sigma})
566 : * = 1/(1 + (e^{log x - mu})^{-1/sigma})
567 : * = 1/(1 + e^{-(log x - mu)/sigma})
568 : * = logistic((log x - mu)/sigma)
569 : * = logistic((log x - log alpha)/(1/beta))
570 : * sf(x) = 1 - 1/(1 + (x/alpha)^-beta)
571 : * = (x/alpha)^-beta/(1 + (x/alpha)^-beta)
572 : * = 1/((x/alpha)^beta + 1)
573 : * = 1/(1 + (x/alpha)^beta)
574 : * icdf(p) = alpha (p/(1 - p))^{1/beta}
575 : * = alpha e^{logit(p)/beta}
576 : * = e^{mu + sigma logit(p)}
577 : * isf(p) = alpha ((1 - p)/p)^{1/beta}
578 : * = alpha e^{-logit(p)/beta}
579 : * = e^{mu - sigma logit(p)}
580 : */
581 :
582 : /**
583 : * Compute the CDF of the LogLogistic(alpha, beta) distribution.
584 : * Well-conditioned for all x and alpha, and the condition number
585 : *
586 : * -beta/[1 + (x/alpha)^{-beta}]
587 : *
588 : * grows linearly with beta.
589 : *
590 : * Loosely, the relative error of this implementation is bounded by
591 : *
592 : * 4 eps + 2 eps^2 + O(beta eps),
593 : *
594 : * so don't bother trying this for beta anywhere near as large as
595 : * 1/eps, around which point it levels off at 1.
596 : */
597 : STATIC double
598 486 : cdf_log_logistic(double x, double alpha, double beta)
599 : {
600 : /*
601 : * Let d0 be the error of x/alpha; d1, of pow; d2, of +; and
602 : * d3, of the final quotient. The exponentiation gives
603 : *
604 : * ((1 + d0) x/alpha)^{-beta}
605 : * = (x/alpha)^{-beta} (1 + d0)^{-beta}
606 : * = (x/alpha)^{-beta} (1 + (1 + d0)^{-beta} - 1)
607 : * = (x/alpha)^{-beta} (1 + d')
608 : *
609 : * where d' = (1 + d0)^{-beta} - 1. If y = (x/alpha)^{-beta},
610 : * the denominator is
611 : *
612 : * (1 + d2) (1 + (1 + d1) (1 + d') y)
613 : * = (1 + d2) (1 + y + (d1 + d' + d1 d') y)
614 : * = 1 + y + (1 + d2) (d1 + d' + d1 d') y
615 : * = (1 + y) (1 + (1 + d2) (d1 + d' + d1 d') y/(1 + y))
616 : * = (1 + y) (1 + d''),
617 : *
618 : * where d'' = (1 + d2) (d1 + d' + d1 d') y/(1 + y). The
619 : * final result is
620 : *
621 : * (1 + d3) / [(1 + d2) (1 + d'') (1 + y)]
622 : * = (1 + d''') / (1 + y)
623 : *
624 : * for |d'''| <= 2|d3 - d''| by Lemma 2 as long as |d''| < 1/2
625 : * (which may not be the case for very large beta). This
626 : * relative error is therefore bounded by
627 : *
628 : * |d'''|
629 : * <= 2|d3 - d''|
630 : * <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d') y/(1 + y)|
631 : * <= 2|d3| + 2|(1 + d2) (d1 + d' + d1 d')|
632 : * = 2|d3| + 2|d1 + d' + d1 d' + d2 d1 + d2 d' + d2 d1 d'|
633 : * <= 2|d3| + 2|d1| + 2|d'| + 2|d1 d'| + 2|d2 d1| + 2|d2 d'|
634 : * + 2|d2 d1 d'|
635 : * <= 4 eps + 2 eps^2 + (2 + 2 eps + 2 eps^2) |d'|.
636 : *
637 : * Roughly, |d'| = |(1 + d0)^{-beta} - 1| grows like beta eps,
638 : * until it levels off at 1.
639 : */
640 486 : return 1/(1 + pow(x/alpha, -beta));
641 : }
642 :
643 : /**
644 : * Compute the SF of the LogLogistic(alpha, beta) distribution.
645 : * Well-conditioned for all x and alpha, and the condition number
646 : *
647 : * beta/[1 + (x/alpha)^beta]
648 : *
649 : * grows linearly with beta.
650 : *
651 : * Loosely, the relative error of this implementation is bounded by
652 : *
653 : * 4 eps + 2 eps^2 + O(beta eps)
654 : *
655 : * so don't bother trying this for beta anywhere near as large as
656 : * 1/eps, beyond which point it grows unbounded.
657 : */
658 : STATIC double
659 1122 : sf_log_logistic(double x, double alpha, double beta)
660 : {
661 : /*
662 : * The error analysis here is essentially the same as in
663 : * cdf_log_logistic, except that rather than levelling off at
664 : * 1, |(1 + d0)^beta - 1| grows unbounded.
665 : */
666 1122 : return 1/(1 + pow(x/alpha, beta));
667 : }
668 :
669 : /**
670 : * Compute the inverse of the CDF of the LogLogistic(alpha, beta)
671 : * distribution. Ill-conditioned for p near 1 and beta near 0 with
672 : * condition number 1/[beta (1 - p)].
673 : */
674 : STATIC double
675 350 : icdf_log_logistic(double p, double alpha, double beta)
676 : {
677 350 : return alpha*pow(p/(1 - p), 1/beta);
678 : }
679 :
680 : /**
681 : * Compute the inverse of the SF of the LogLogistic(alpha, beta)
682 : * distribution. Ill-conditioned for p near 1 and for large beta, with
683 : * condition number -1/[beta (1 - p)].
684 : */
685 : STATIC double
686 346 : isf_log_logistic(double p, double alpha, double beta)
687 : {
688 346 : return alpha*pow((1 - p)/p, 1/beta);
689 : }
690 :
691 : /*
692 : * Weibull(lambda, k) distribution, supported on (0, +infinity).
693 : *
694 : * pdf(x) = (k/lambda) (x/lambda)^{k - 1} e^{-(x/lambda)^k}
695 : * cdf(x) = 1 - e^{-(x/lambda)^k}
696 : * icdf(p) = lambda * (-log (1 - p))^{1/k}
697 : * sf(x) = e^{-(x/lambda)^k}
698 : * isf(p) = lambda * (-log p)^{1/k}
699 : */
700 :
701 : /**
702 : * Compute the CDF of the Weibull(lambda, k) distribution.
703 : * Well-conditioned for small x and k, and for large lambda --
704 : * condition number
705 : *
706 : * -k (x/lambda)^k exp(-(x/lambda)^k)/[exp(-(x/lambda)^k) - 1]
707 : *
708 : * grows linearly with k, x^k, and lambda^{-k}.
709 : */
710 : STATIC double
711 469 : cdf_weibull(double x, double lambda, double k)
712 : {
713 469 : return -expm1(-pow(x/lambda, k));
714 : }
715 :
716 : /**
717 : * Compute the SF of the Weibull(lambda, k) distribution.
718 : * Well-conditioned for small x and k, and for large lambda --
719 : * condition number
720 : *
721 : * -k (x/lambda)^k
722 : *
723 : * grows linearly with k, x^k, and lambda^{-k}.
724 : */
725 : STATIC double
726 1049 : sf_weibull(double x, double lambda, double k)
727 : {
728 1049 : return exp(-pow(x/lambda, k));
729 : }
730 :
731 : /**
732 : * Compute the inverse of the CDF of the Weibull(lambda, k)
733 : * distribution. Ill-conditioned for p near 1, and for k near 0;
734 : * condition number is
735 : *
736 : * (p/(1 - p))/(k log(1 - p)).
737 : */
738 : STATIC double
739 52 : icdf_weibull(double p, double lambda, double k)
740 : {
741 52 : return lambda*pow(-log1p(-p), 1/k);
742 : }
743 :
744 : /**
745 : * Compute the inverse of the SF of the Weibull(lambda, k)
746 : * distribution. Ill-conditioned for p near 0, and for k near 0;
747 : * condition number is
748 : *
749 : * 1/(k log(p)).
750 : */
751 : STATIC double
752 65 : isf_weibull(double p, double lambda, double k)
753 : {
754 65 : return lambda*pow(-log(p), 1/k);
755 : }
756 :
757 : /*
758 : * GeneralizedPareto(mu, sigma, xi), supported on (mu, +infinity) for
759 : * nonnegative xi, or (mu, mu - sigma/xi) for negative xi.
760 : *
761 : * Samples:
762 : * = mu - sigma log U, if xi = 0;
763 : * = mu + sigma (U^{-xi} - 1)/xi = mu + sigma*expm1(-xi log U)/xi, if xi =/= 0,
764 : * where U is uniform on (0,1].
765 : * = mu + sigma (e^{xi X} - 1)/xi,
766 : * where X has standard exponential distribution.
767 : *
768 : * pdf(x) = sigma^{-1} (1 + xi (x - mu)/sigma)^{-(1 + 1/xi)}
769 : * cdf(x) = 1 - (1 + xi (x - mu)/sigma)^{-1/xi}
770 : * = 1 - e^{-log(1 + xi (x - mu)/sigma)/xi}
771 : * --> 1 - e^{-(x - mu)/sigma} as xi --> 0
772 : * sf(x) = (1 + xi (x - mu)/sigma)^{-1/xi}
773 : * --> e^{-(x - mu)/sigma} as xi --> 0
774 : * icdf(p) = mu + sigma*(p^{-xi} - 1)/xi
775 : * = mu + sigma*expm1(-xi log p)/xi
776 : * --> mu + sigma*log p as xi --> 0
777 : * isf(p) = mu + sigma*((1 - p)^{xi} - 1)/xi
778 : * = mu + sigma*expm1(-xi log1p(-p))/xi
779 : * --> mu + sigma*log1p(-p) as xi --> 0
780 : */
781 :
782 : /**
783 : * Compute the CDF of the GeneralizedPareto(mu, sigma, xi)
784 : * distribution. Well-conditioned everywhere. For standard
785 : * distribution (mu=0, sigma=1), condition number
786 : *
787 : * (x/(1 + x xi)) / ((1 + x xi)^{1/xi} - 1)
788 : *
789 : * is bounded by 1, attained only at x = 0.
790 : */
791 : STATIC double
792 379 : cdf_genpareto(double x, double mu, double sigma, double xi)
793 : {
794 379 : double x_0 = (x - mu)/sigma;
795 :
796 : /*
797 : * log(1 + xi x_0)/xi
798 : * = (-1/xi) \sum_{n=1}^infinity (-xi x_0)^n/n
799 : * = (-1/xi) (-xi x_0 + \sum_{n=2}^infinity (-xi x_0)^n/n)
800 : * = x_0 - (1/xi) \sum_{n=2}^infinity (-xi x_0)^n/n
801 : * = x_0 - x_0 \sum_{n=2}^infinity (-xi x_0)^{n-1}/n
802 : * = x_0 (1 - d),
803 : *
804 : * where d = \sum_{n=2}^infinity (-xi x_0)^{n-1}/n. If |xi| <
805 : * eps/4|x_0|, then
806 : *
807 : * |d| <= \sum_{n=2}^infinity (eps/4)^{n-1}/n
808 : * <= \sum_{n=2}^infinity (eps/4)^{n-1}
809 : * = \sum_{n=1}^infinity (eps/4)^n
810 : * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
811 : * = (eps/4)/(1 - eps/4)
812 : * < eps/2
813 : *
814 : * for any 0 < eps < 2. Thus, the relative error of x_0 from
815 : * log(1 + xi x_0)/xi is bounded by eps.
816 : */
817 379 : if (fabs(xi) < 1e-17/x_0)
818 153 : return -expm1(-x_0);
819 : else
820 226 : return -expm1(-log1p(xi*x_0)/xi);
821 : }
822 :
823 : /**
824 : * Compute the SF of the GeneralizedPareto(mu, sigma, xi) distribution.
825 : * For standard distribution (mu=0, sigma=1), ill-conditioned for xi
826 : * near 0; condition number
827 : *
828 : * -x (1 + x xi)^{(-1 - xi)/xi}/(1 + x xi)^{-1/xi}
829 : * = -x (1 + x xi)^{-1/xi - 1}/(1 + x xi)^{-1/xi}
830 : * = -(x/(1 + x xi)) (1 + x xi)^{-1/xi}/(1 + x xi)^{-1/xi}
831 : * = -x/(1 + x xi)
832 : *
833 : * is bounded by 1/xi.
834 : */
835 : STATIC double
836 1343 : sf_genpareto(double x, double mu, double sigma, double xi)
837 : {
838 1343 : double x_0 = (x - mu)/sigma;
839 :
840 1343 : if (fabs(xi) < 1e-17/x_0)
841 573 : return exp(-x_0);
842 : else
843 770 : return exp(-log1p(xi*x_0)/xi);
844 : }
845 :
846 : /**
847 : * Compute the inverse of the CDF of the GeneralizedPareto(mu, sigma,
848 : * xi) distribution. Ill-conditioned for p near 1; condition number is
849 : *
850 : * xi (p/(1 - p))/(1 - (1 - p)^xi)
851 : */
852 : STATIC double
853 1542 : icdf_genpareto(double p, double mu, double sigma, double xi)
854 : {
855 : /*
856 : * To compute f(xi) = (U^{-xi} - 1)/xi = (e^{-xi log U} - 1)/xi
857 : * for xi near zero (note f(xi) --> -log U as xi --> 0), write
858 : * the absolutely convergent Taylor expansion
859 : *
860 : * f(xi) = (1/xi)*(-xi log U + \sum_{n=2}^infinity (-xi log U)^n/n!
861 : * = -log U + (1/xi)*\sum_{n=2}^infinity (-xi log U)^n/n!
862 : * = -log U + \sum_{n=2}^infinity xi^{n-1} (-log U)^n/n!
863 : * = -log U - log U \sum_{n=2}^infinity (-xi log U)^{n-1}/n!
864 : * = -log U (1 + \sum_{n=2}^infinity (-xi log U)^{n-1}/n!).
865 : *
866 : * Let d = \sum_{n=2}^infinity (-xi log U)^{n-1}/n!. What do we
867 : * lose if we discard it and use -log U as an approximation to
868 : * f(xi)? If |xi| < eps/-4log U, then
869 : *
870 : * |d| <= \sum_{n=2}^infinity |xi log U|^{n-1}/n!
871 : * <= \sum_{n=2}^infinity (eps/4)^{n-1}/n!
872 : * <= \sum_{n=1}^infinity (eps/4)^n
873 : * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
874 : * = (eps/4)/(1 - eps/4)
875 : * < eps/2,
876 : *
877 : * for any 0 < eps < 2. Hence, as long as |xi| < eps/-2log U,
878 : * f(xi) = -log U (1 + d) for |d| <= eps/2. |d| is the
879 : * relative error of f(xi) from -log U; from this bound, the
880 : * relative error of -log U from f(xi) is at most (eps/2)/(1 -
881 : * eps/2) = eps/2 + (eps/2)^2 + (eps/2)^3 + ... < eps for 0 <
882 : * eps < 1. Since -log U < 1000 for all U in (0, 1] in
883 : * binary64 floating-point, we can safely cut xi off at 1e-20 <
884 : * eps/4000 and attain <1ulp error from series truncation.
885 : */
886 1542 : if (fabs(xi) <= 1e-20)
887 636 : return mu - sigma*log1p(-p);
888 : else
889 906 : return mu + sigma*expm1(-xi*log1p(-p))/xi;
890 : }
891 :
892 : /**
893 : * Compute the inverse of the SF of the GeneralizedPareto(mu, sigma,
894 : * xi) distribution. Ill-conditioned for p near 1; condition number is
895 : *
896 : * -xi/(1 - p^{-xi})
897 : */
898 : STATIC double
899 1469 : isf_genpareto(double p, double mu, double sigma, double xi)
900 : {
901 1469 : if (fabs(xi) <= 1e-20)
902 621 : return mu - sigma*log(p);
903 : else
904 848 : return mu + sigma*expm1(-xi*log(p))/xi;
905 : }
906 :
907 : /*******************************************************************/
908 :
909 : /**
910 : * Deterministic samplers, parametrized by uniform integer and (0,1]
911 : * samples. No guarantees are made about _which_ mapping from the
912 : * integer and (0,1] samples these use; all that is guaranteed is the
913 : * distribution of the outputs conditioned on a uniform distribution on
914 : * the inputs. The automatic tests in test_prob_distr.c double-check
915 : * the particular mappings we use.
916 : *
917 : * Beware: Unlike random_uniform_01(), these are not guaranteed to be
918 : * supported on all possible outputs. See Ilya Mironov, `On the
919 : * Significance of the Least Significant Bits for Differential
920 : * Privacy', for an example of what can go wrong if you try to use
921 : * these to conceal information from an adversary but you expose the
922 : * specific full-precision floating-point values.
923 : *
924 : * Note: None of these samplers use rejection sampling; they are all
925 : * essentially inverse-CDF transforms with tweaks. If you were to add,
926 : * say, a Gamma sampler with the Marsaglia-Tsang method, you would have
927 : * to parametrize it by a potentially infinite stream of uniform (and
928 : * perhaps normal) samples rather than a fixed number, which doesn't
929 : * make for quite as nice automatic testing as for these.
930 : */
931 :
932 : /**
933 : * Deterministically sample from the interval [a, b], indexed by a
934 : * uniform random floating-point number p0 in (0, 1].
935 : *
936 : * Note that even if p0 is nonzero, the result may be equal to a, if
937 : * ulp(a)/2 is nonnegligible, e.g. if a = 1. For maximum resolution,
938 : * arrange |a| <= |b|.
939 : */
940 : STATIC double
941 600178 : sample_uniform_interval(double p0, double a, double b)
942 : {
943 : /*
944 : * XXX Prove that the distribution is, in fact, uniform on
945 : * [a,b], particularly around p0 = 1, or at least has very
946 : * small deviation from uniform, quantified appropriately
947 : * (e.g., like in Monahan 1984, or by KL divergence). It
948 : * almost certainly does but it would be nice to quantify the
949 : * error.
950 : */
951 600178 : if ((a <= 0 && 0 <= b) || (b <= 0 && 0 <= a)) {
952 : /*
953 : * When ab < 0, (1 - t) a + t b is monotonic, since for
954 : * a <= b it is a sum of nondecreasing functions of t,
955 : * and for b <= a, of nonincreasing functions of t.
956 : * Further, clearly at 0 and 1 it attains a and b,
957 : * respectively. Hence it is bounded within [a, b].
958 : */
959 300140 : return (1 - p0)*a + p0*b;
960 : } else {
961 : /*
962 : * a + (b - a) t is monotonic -- it is obviously a
963 : * nondecreasing function of t for a <= b. Further, it
964 : * attains a at 0, and while it may overshoot b at 1,
965 : * we have a
966 : *
967 : * Theorem. If 0 <= t < 1, then the floating-point
968 : * evaluation of a + (b - a) t is bounded in [a, b].
969 : *
970 : * Lemma 1. If 0 <= t < 1 is a floating-point number,
971 : * then for any normal floating-point number x except
972 : * the smallest in magnitude, |round(x*t)| < |x|.
973 : *
974 : * Proof. WLOG, assume x >= 0. Since the rounding
975 : * function and t |---> x*t are nondecreasing, their
976 : * composition t |---> round(x*t) is also
977 : * nondecreasing, so it suffices to consider the
978 : * largest floating-point number below 1, in particular
979 : * t = 1 - ulp(1)/2.
980 : *
981 : * Case I: If x is a power of two, then the next
982 : * floating-point number below x is x - ulp(x)/2 = x -
983 : * x*ulp(1)/2 = x*(1 - ulp(1)/2) = x*t, so, since x*t
984 : * is a floating-point number, multiplication is exact,
985 : * and thus round(x*t) = x*t < x.
986 : *
987 : * Case II: If x is not a power of two, then the
988 : * greatest lower bound of real numbers rounded to x is
989 : * x - ulp(x)/2 = x - ulp(T(x))/2 = x - T(x)*ulp(1)/2,
990 : * where T(X) is the largest power of two below x.
991 : * Anything below this bound is rounded to a
992 : * floating-point number smaller than x, and x*t = x*(1
993 : * - ulp(1)/2) = x - x*ulp(1)/2 < x - T(x)*ulp(1)/2
994 : * since T(x) < x, so round(x*t) < x*t < x. QED.
995 : *
996 : * Lemma 2. If x and y are subnormal, then round(x +
997 : * y) = x + y.
998 : *
999 : * Proof. It is a matter of adding the significands,
1000 : * since if we treat subnormals as having an implicit
1001 : * zero bit before the `binary' point, their exponents
1002 : * are all the same. There is at most one carry/borrow
1003 : * bit, which can always be accommodated either in a
1004 : * subnormal, or, at largest, in the implicit one bit
1005 : * of a normal.
1006 : *
1007 : * Lemma 3. Let x and y be floating-point numbers. If
1008 : * round(x - y) is subnormal or zero, then it is equal
1009 : * to x - y.
1010 : *
1011 : * Proof. Case I (equal): round(x - y) = 0 iff x = y;
1012 : * hence if round(x - y) = 0, then round(x - y) = 0 = x
1013 : * - y.
1014 : *
1015 : * Case II (subnormal/subnormal): If x and y are both
1016 : * subnormal, this follows directly from Lemma 2.
1017 : *
1018 : * Case IIIa (normal/subnormal): If x is normal and y
1019 : * is subnormal, then x and y must share sign, or else
1020 : * x - y would be larger than x and thus rounded to
1021 : * normal. If s is the smallest normal positive
1022 : * floating-point number, |x| < 2s since by
1023 : * construction 2s - |y| is normal for all subnormal y.
1024 : * This means that x and y must have the same exponent,
1025 : * so the difference is the difference of significands,
1026 : * which is exact.
1027 : *
1028 : * Case IIIb (subnormal/normal): Same as case IIIa for
1029 : * -(y - x).
1030 : *
1031 : * Case IV (normal/normal): If x and y are both normal,
1032 : * then they must share sign, or else x - y would be
1033 : * larger than x and thus rounded to normal. Note that
1034 : * |y| < 2|x|, for if |y| >= 2|x|, then |x| - |y| <=
1035 : * -|x| but -|x| is normal like x. Also, |x|/2 < |y|:
1036 : * if |x|/2 is subnormal, it must hold because y is
1037 : * normal; if |x|/2 is normal, then |x|/2 >= s, so
1038 : * since |x| - |y| < s,
1039 : *
1040 : * |x|/2 = |x| - |x|/2 <= |x| - s <= |y|;
1041 : *
1042 : * that is, |x|/2 < |y| < 2|x|, so by the Sterbenz
1043 : * lemma, round(x - y) = x - y. QED.
1044 : *
1045 : * Proof of theorem. WLOG, assume 0 <= a <= b. Since
1046 : * round(a + round(round(b - a)*t) is nondecreasing in
1047 : * t and attains a at 0, the lower end of the bound is
1048 : * trivial; we must show the upper end of the bound
1049 : * strictly. It suffices to show this for the largest
1050 : * floating-point number below 1, namely 1 - ulp(1)/2.
1051 : *
1052 : * Case I: round(b - a) is normal. Then it is at most
1053 : * the smallest floating-point number above b - a. By
1054 : * Lemma 1, round(round(b - a)*t) < round(b - a).
1055 : * Since the inequality is strict, and since
1056 : * round(round(b - a)*t) is a floating-point number
1057 : * below round(b - a), and since there are no
1058 : * floating-point numbers between b - a and round(b -
1059 : * a), we must have round(round(b - a)*t) < b - a.
1060 : * Then since y |---> round(a + y) is nondecreasing, we
1061 : * must have
1062 : *
1063 : * round(a + round(round(b - a)*t))
1064 : * <= round(a + (b - a))
1065 : * = round(b) = b.
1066 : *
1067 : * Case II: round(b - a) is subnormal. In this case,
1068 : * Lemma 1 falls apart -- we are not guaranteed the
1069 : * strict inequality. However, by Lemma 3, the
1070 : * difference is exact: round(b - a) = b - a. Thus,
1071 : *
1072 : * round(a + round(round(b - a)*t))
1073 : * <= round(a + round((b - a)*t))
1074 : * <= round(a + (b - a))
1075 : * = round(b)
1076 : * = b,
1077 : *
1078 : * QED.
1079 : */
1080 :
1081 : /* p0 is restricted to [0,1], but we use >= to silence -Wfloat-equal. */
1082 300038 : if (p0 >= 1)
1083 : return b;
1084 300030 : return a + (b - a)*p0;
1085 : }
1086 : }
1087 :
1088 : /**
1089 : * Deterministically sample from the standard logistic distribution,
1090 : * indexed by a uniform random 32-bit integer s and uniform random
1091 : * floating-point numbers t and p0 in (0, 1].
1092 : */
1093 : STATIC double
1094 400908 : sample_logistic(uint32_t s, double t, double p0)
1095 : {
1096 400908 : double sign = (s & 1) ? -1 : +1;
1097 400908 : double r;
1098 :
1099 : /*
1100 : * We carve up the interval (0, 1) into subregions to compute
1101 : * the inverse CDF precisely:
1102 : *
1103 : * A = (0, 1/(1 + e)] ---> (-infinity, -1]
1104 : * B = [1/(1 + e), 1/2] ---> [-1, 0]
1105 : * C = [1/2, 1 - 1/(1 + e)] ---> [0, 1]
1106 : * D = [1 - 1/(1 + e), 1) ---> [1, +infinity)
1107 : *
1108 : * Cases D and C are mirror images of cases A and B,
1109 : * respectively, so we choose between them by the sign chosen
1110 : * by a fair coin toss. We choose between cases A and B by a
1111 : * coin toss weighted by
1112 : *
1113 : * 2/(1 + e) = 1 - [1/2 - 1/(1 + e)]/(1/2):
1114 : *
1115 : * if it comes up heads, scale p0 into a uniform (0, 1/(1 + e)]
1116 : * sample p; if it comes up tails, scale p0 into a uniform (0,
1117 : * 1/2 - 1/(1 + e)] sample and compute the inverse CDF of p =
1118 : * 1/2 - p0.
1119 : */
1120 400908 : if (t <= 2/(1 + exp(1))) {
1121 : /* p uniform in (0, 1/(1 + e)], represented by p. */
1122 215312 : p0 /= 1 + exp(1);
1123 215312 : r = logit(p0);
1124 : } else {
1125 : /*
1126 : * p uniform in [1/(1 + e), 1/2), actually represented
1127 : * by p0 = 1/2 - p uniform in (0, 1/2 - 1/(1 + e)], so
1128 : * that p = 1/2 - p.
1129 : */
1130 185596 : p0 *= 0.5 - 1/(1 + exp(1));
1131 185596 : r = logithalf(p0);
1132 : }
1133 :
1134 : /*
1135 : * We have chosen from the negative half of the standard
1136 : * logistic distribution, which is symmetric with the positive
1137 : * half. Now use the sign to choose uniformly between them.
1138 : */
1139 400908 : return sign*r;
1140 : }
1141 :
1142 : /**
1143 : * Deterministically sample from the logistic distribution scaled by
1144 : * sigma and translated by mu.
1145 : */
1146 : static double
1147 400100 : sample_logistic_locscale(uint32_t s, double t, double p0, double mu,
1148 : double sigma)
1149 : {
1150 :
1151 400100 : return mu + sigma*sample_logistic(s, t, p0);
1152 : }
1153 :
1154 : /**
1155 : * Deterministically sample from the standard log-logistic
1156 : * distribution, indexed by a uniform random 32-bit integer s and a
1157 : * uniform random floating-point number p0 in (0, 1].
1158 : */
1159 : STATIC double
1160 400504 : sample_log_logistic(uint32_t s, double p0)
1161 : {
1162 :
1163 : /*
1164 : * Carve up the interval (0, 1) into (0, 1/2] and [1/2, 1); the
1165 : * condition numbers of the icdf and the isf coincide at 1/2.
1166 : */
1167 400504 : p0 *= 0.5;
1168 400504 : if ((s & 1) == 0) {
1169 : /* p = p0 in (0, 1/2] */
1170 200463 : return p0/(1 - p0);
1171 : } else {
1172 : /* p = 1 - p0 in [1/2, 1) */
1173 200041 : return (1 - p0)/p0;
1174 : }
1175 : }
1176 :
1177 : /**
1178 : * Deterministically sample from the log-logistic distribution with
1179 : * scale alpha and shape beta.
1180 : */
1181 : static double
1182 400100 : sample_log_logistic_scaleshape(uint32_t s, double p0, double alpha,
1183 : double beta)
1184 : {
1185 400100 : double x = sample_log_logistic(s, p0);
1186 :
1187 400100 : return alpha*pow(x, 1/beta);
1188 : }
1189 :
1190 : /**
1191 : * Deterministically sample from the standard exponential distribution,
1192 : * indexed by a uniform random 32-bit integer s and a uniform random
1193 : * floating-point number p0 in (0, 1].
1194 : */
1195 : static double
1196 1606360 : sample_exponential(uint32_t s, double p0)
1197 : {
1198 : /*
1199 : * We would like to evaluate log(p) for p near 0, and log1p(-p)
1200 : * for p near 1. Simply carve the interval into (0, 1/2] and
1201 : * [1/2, 1) by a fair coin toss.
1202 : */
1203 1606360 : p0 *= 0.5;
1204 1606360 : if ((s & 1) == 0)
1205 : /* p = p0 in (0, 1/2] */
1206 803352 : return -log(p0);
1207 : else
1208 : /* p = 1 - p0 in [1/2, 1) */
1209 803008 : return -log1p(-p0);
1210 : }
1211 :
1212 : /**
1213 : * Deterministically sample from a Weibull distribution with scale
1214 : * lambda and shape k -- just an exponential with a shape parameter in
1215 : * addition to a scale parameter. (Yes, lambda really is the scale,
1216 : * _not_ the rate.)
1217 : */
1218 : STATIC double
1219 500504 : sample_weibull(uint32_t s, double p0, double lambda, double k)
1220 : {
1221 :
1222 500504 : return lambda*pow(sample_exponential(s, p0), 1/k);
1223 : }
1224 :
1225 : /**
1226 : * Deterministically sample from the generalized Pareto distribution
1227 : * with shape xi, indexed by a uniform random 32-bit integer s and a
1228 : * uniform random floating-point number p0 in (0, 1].
1229 : */
1230 : STATIC double
1231 705756 : sample_genpareto(uint32_t s, double p0, double xi)
1232 : {
1233 705756 : double x = sample_exponential(s, p0);
1234 :
1235 : /*
1236 : * Write f(xi) = (e^{xi x} - 1)/xi for xi near zero as the
1237 : * absolutely convergent Taylor series
1238 : *
1239 : * f(x) = (1/xi) (xi x + \sum_{n=2}^infinity (xi x)^n/n!)
1240 : * = x + (1/xi) \sum_{n=2}^infinity (xi x)^n/n!
1241 : * = x + \sum_{n=2}^infinity xi^{n-1} x^n/n!
1242 : * = x + x \sum_{n=2}^infinity (xi x)^{n-1}/n!
1243 : * = x (1 + \sum_{n=2}^infinity (xi x)^{n-1}/n!).
1244 : *
1245 : * d = \sum_{n=2}^infinity (xi x)^{n-1}/n! is the relative error
1246 : * of f(x) from x. If |xi| < eps/4x, then
1247 : *
1248 : * |d| <= \sum_{n=2}^infinity |xi x|^{n-1}/n!
1249 : * <= \sum_{n=2}^infinity (eps/4)^{n-1}/n!
1250 : * <= \sum_{n=1}^infinity (eps/4)
1251 : * = (eps/4) \sum_{n=0}^infinity (eps/4)^n
1252 : * = (eps/4)/(1 - eps/4)
1253 : * < eps/2,
1254 : *
1255 : * for any 0 < eps < 2. Hence, as long as |xi| < eps/2x, f(xi)
1256 : * = x (1 + d) for |d| <= eps/2, so x = f(xi) (1 + d') for |d'|
1257 : * <= eps. What bound should we use for x?
1258 : *
1259 : * - If x is exponentially distributed, x > 200 with
1260 : * probability below e^{-200} << 2^{-256}, i.e. never.
1261 : *
1262 : * - If x is computed by -log(U) for U in (0, 1], x is
1263 : * guaranteed to be below 1000 in IEEE 754 binary64
1264 : * floating-point.
1265 : *
1266 : * We can safely cut xi off at 1e-20 < eps/4000 and attain an
1267 : * error bounded by 0.5 ulp for this expression.
1268 : */
1269 705756 : return (fabs(xi) < 1e-20 ? x : expm1(xi*x)/xi);
1270 : }
1271 :
1272 : /**
1273 : * Deterministically sample from a generalized Pareto distribution with
1274 : * shape xi, scaled by sigma and translated by mu.
1275 : */
1276 : static double
1277 700100 : sample_genpareto_locscale(uint32_t s, double p0, double mu, double sigma,
1278 : double xi)
1279 : {
1280 :
1281 700100 : return mu + sigma*sample_genpareto(s, p0, xi);
1282 : }
1283 :
1284 : /**
1285 : * Deterministically sample from the geometric distribution with
1286 : * per-trial success probability p.
1287 : **/
1288 : // clang-format off
1289 : /*
1290 : * XXX Quantify the error (KL divergence?) of this
1291 : * ceiling-of-exponential sampler from a true geometric distribution,
1292 : * which we could get by rejection sampling. Relevant papers:
1293 : *
1294 : * John F. Monahan, `Accuracy in Random Number Generation',
1295 : * Mathematics of Computation 45(172), October 1984, pp. 559--568.
1296 : https://pdfs.semanticscholar.org/aca6/74b96da1df77b2224e8cfc5dd6d61a471632.pdf
1297 : * Karl Bringmann and Tobias Friedrich, `Exact and Efficient
1298 : * Generation of Geometric Random Variates and Random Graphs', in
1299 : * Proceedings of the 40th International Colloaquium on Automata,
1300 : * Languages, and Programming -- ICALP 2013, Springer LNCS 7965,
1301 : * pp.267--278.
1302 : * https://doi.org/10.1007/978-3-642-39206-1_23
1303 : * https://people.mpi-inf.mpg.de/~kbringma/paper/2013ICALP-1.pdf
1304 : */
1305 : // clang-format on
1306 : static double
1307 400100 : sample_geometric(uint32_t s, double p0, double p)
1308 : {
1309 400100 : double x = sample_exponential(s, p0);
1310 :
1311 : /* This is actually a check against 1, but we do >= so that the compiler
1312 : does not raise a -Wfloat-equal */
1313 400100 : if (p >= 1)
1314 : return 1;
1315 :
1316 300100 : return ceil(-x/log1p(-p));
1317 : }
1318 :
1319 : /*******************************************************************/
1320 :
1321 : /** Public API for probability distributions:
1322 : *
1323 : * These are wrapper functions on top of the various probability distribution
1324 : * operations using the generic <b>dist</b> structure.
1325 :
1326 : * These are the functions that should be used by consumers of this API.
1327 : */
1328 :
1329 : /** Returns the name of the distribution in <b>dist</b>. */
1330 : const char *
1331 0 : dist_name(const struct dist_t *dist)
1332 : {
1333 0 : return dist->ops->name;
1334 : }
1335 :
1336 : /* Sample a value from <b>dist</b> and return it. */
1337 : double
1338 3000614 : dist_sample(const struct dist_t *dist)
1339 : {
1340 3000614 : return dist->ops->sample(dist);
1341 : }
1342 :
1343 : /** Compute the CDF of <b>dist</b> at <b>x</b>. */
1344 : double
1345 1448 : dist_cdf(const struct dist_t *dist, double x)
1346 : {
1347 1448 : return dist->ops->cdf(dist, x);
1348 : }
1349 :
1350 : /** Compute the SF (Survival function) of <b>dist</b> at <b>x</b>. */
1351 : double
1352 3700 : dist_sf(const struct dist_t *dist, double x)
1353 : {
1354 3700 : return dist->ops->sf(dist, x);
1355 : }
1356 :
1357 : /** Compute the iCDF (Inverse CDF) of <b>dist</b> at <b>x</b>. */
1358 : double
1359 52 : dist_icdf(const struct dist_t *dist, double p)
1360 : {
1361 52 : return dist->ops->icdf(dist, p);
1362 : }
1363 :
1364 : /** Compute the iSF (Inverse Survival function) of <b>dist</b> at <b>x</b>. */
1365 : double
1366 26 : dist_isf(const struct dist_t *dist, double p)
1367 : {
1368 26 : return dist->ops->isf(dist, p);
1369 : }
1370 :
1371 : /** Functions for uniform distribution */
1372 :
1373 : static double
1374 600114 : uniform_sample(const struct dist_t *dist)
1375 : {
1376 600114 : const struct uniform_t *U = dist_to_const_uniform(dist);
1377 600114 : double p0 = random_uniform_01();
1378 :
1379 600114 : return sample_uniform_interval(p0, U->a, U->b);
1380 : }
1381 :
1382 : static double
1383 584 : uniform_cdf(const struct dist_t *dist, double x)
1384 : {
1385 584 : const struct uniform_t *U = dist_to_const_uniform(dist);
1386 584 : if (x < U->a)
1387 : return 0;
1388 584 : else if (x < U->b)
1389 584 : return (x - U->a)/(U->b - U->a);
1390 : else
1391 : return 1;
1392 : }
1393 :
1394 : static double
1395 604 : uniform_sf(const struct dist_t *dist, double x)
1396 : {
1397 604 : const struct uniform_t *U = dist_to_const_uniform(dist);
1398 :
1399 604 : if (x > U->b)
1400 : return 0;
1401 604 : else if (x > U->a)
1402 604 : return (U->b - x)/(U->b - U->a);
1403 : else
1404 : return 1;
1405 : }
1406 :
1407 : static double
1408 12 : uniform_icdf(const struct dist_t *dist, double p)
1409 : {
1410 12 : const struct uniform_t *U = dist_to_const_uniform(dist);
1411 12 : double w = U->b - U->a;
1412 :
1413 12 : return (p < 0.5 ? (U->a + w*p) : (U->b - w*(1 - p)));
1414 : }
1415 :
1416 : static double
1417 6 : uniform_isf(const struct dist_t *dist, double p)
1418 : {
1419 6 : const struct uniform_t *U = dist_to_const_uniform(dist);
1420 6 : double w = U->b - U->a;
1421 :
1422 6 : return (p < 0.5 ? (U->b - w*p) : (U->a + w*(1 - p)));
1423 : }
1424 :
1425 : const struct dist_ops_t uniform_ops = {
1426 : .name = "uniform",
1427 : .sample = uniform_sample,
1428 : .cdf = uniform_cdf,
1429 : .sf = uniform_sf,
1430 : .icdf = uniform_icdf,
1431 : .isf = uniform_isf,
1432 : };
1433 :
1434 : /*******************************************************************/
1435 :
1436 : /** Private functions for each probability distribution. */
1437 :
1438 : /** Functions for logistic distribution: */
1439 :
1440 : static double
1441 400100 : logistic_sample(const struct dist_t *dist)
1442 : {
1443 400100 : const struct logistic_t *L = dist_to_const_logistic(dist);
1444 400100 : uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng());
1445 400100 : double t = random_uniform_01();
1446 400100 : double p0 = random_uniform_01();
1447 :
1448 400100 : return sample_logistic_locscale(s, t, p0, L->mu, L->sigma);
1449 : }
1450 :
1451 : static double
1452 388 : logistic_cdf(const struct dist_t *dist, double x)
1453 : {
1454 388 : const struct logistic_t *L = dist_to_const_logistic(dist);
1455 388 : return cdf_logistic(x, L->mu, L->sigma);
1456 : }
1457 :
1458 : static double
1459 404 : logistic_sf(const struct dist_t *dist, double x)
1460 : {
1461 404 : const struct logistic_t *L = dist_to_const_logistic(dist);
1462 404 : return sf_logistic(x, L->mu, L->sigma);
1463 : }
1464 :
1465 : static double
1466 8 : logistic_icdf(const struct dist_t *dist, double p)
1467 : {
1468 8 : const struct logistic_t *L = dist_to_const_logistic(dist);
1469 8 : return icdf_logistic(p, L->mu, L->sigma);
1470 : }
1471 :
1472 : static double
1473 4 : logistic_isf(const struct dist_t *dist, double p)
1474 : {
1475 4 : const struct logistic_t *L = dist_to_const_logistic(dist);
1476 4 : return isf_logistic(p, L->mu, L->sigma);
1477 : }
1478 :
1479 : const struct dist_ops_t logistic_ops = {
1480 : .name = "logistic",
1481 : .sample = logistic_sample,
1482 : .cdf = logistic_cdf,
1483 : .sf = logistic_sf,
1484 : .icdf = logistic_icdf,
1485 : .isf = logistic_isf,
1486 : };
1487 :
1488 : /** Functions for log-logistic distribution: */
1489 :
1490 : static double
1491 400100 : log_logistic_sample(const struct dist_t *dist)
1492 : {
1493 400100 : const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1494 400100 : uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng());
1495 400100 : double p0 = random_uniform_01();
1496 :
1497 400100 : return sample_log_logistic_scaleshape(s, p0, LL->alpha, LL->beta);
1498 : }
1499 :
1500 : static double
1501 78 : log_logistic_cdf(const struct dist_t *dist, double x)
1502 : {
1503 78 : const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1504 78 : return cdf_log_logistic(x, LL->alpha, LL->beta);
1505 : }
1506 :
1507 : static double
1508 714 : log_logistic_sf(const struct dist_t *dist, double x)
1509 : {
1510 714 : const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1511 714 : return sf_log_logistic(x, LL->alpha, LL->beta);
1512 : }
1513 :
1514 : static double
1515 8 : log_logistic_icdf(const struct dist_t *dist, double p)
1516 : {
1517 8 : const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1518 8 : return icdf_log_logistic(p, LL->alpha, LL->beta);
1519 : }
1520 :
1521 : static double
1522 4 : log_logistic_isf(const struct dist_t *dist, double p)
1523 : {
1524 4 : const struct log_logistic_t *LL = dist_to_const_log_logistic(dist);
1525 4 : return isf_log_logistic(p, LL->alpha, LL->beta);
1526 : }
1527 :
1528 : const struct dist_ops_t log_logistic_ops = {
1529 : .name = "log logistic",
1530 : .sample = log_logistic_sample,
1531 : .cdf = log_logistic_cdf,
1532 : .sf = log_logistic_sf,
1533 : .icdf = log_logistic_icdf,
1534 : .isf = log_logistic_isf,
1535 : };
1536 :
1537 : /** Functions for Weibull distribution */
1538 :
1539 : static double
1540 500100 : weibull_sample(const struct dist_t *dist)
1541 : {
1542 500100 : const struct weibull_t *W = dist_to_const_weibull(dist);
1543 500100 : uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng());
1544 500100 : double p0 = random_uniform_01();
1545 :
1546 500100 : return sample_weibull(s, p0, W->lambda, W->k);
1547 : }
1548 :
1549 : static double
1550 187 : weibull_cdf(const struct dist_t *dist, double x)
1551 : {
1552 187 : const struct weibull_t *W = dist_to_const_weibull(dist);
1553 187 : return cdf_weibull(x, W->lambda, W->k);
1554 : }
1555 :
1556 : static double
1557 803 : weibull_sf(const struct dist_t *dist, double x)
1558 : {
1559 803 : const struct weibull_t *W = dist_to_const_weibull(dist);
1560 803 : return sf_weibull(x, W->lambda, W->k);
1561 : }
1562 :
1563 : static double
1564 10 : weibull_icdf(const struct dist_t *dist, double p)
1565 : {
1566 10 : const struct weibull_t *W = dist_to_const_weibull(dist);
1567 10 : return icdf_weibull(p, W->lambda, W->k);
1568 : }
1569 :
1570 : static double
1571 5 : weibull_isf(const struct dist_t *dist, double p)
1572 : {
1573 5 : const struct weibull_t *W = dist_to_const_weibull(dist);
1574 5 : return isf_weibull(p, W->lambda, W->k);
1575 : }
1576 :
1577 : const struct dist_ops_t weibull_ops = {
1578 : .name = "Weibull",
1579 : .sample = weibull_sample,
1580 : .cdf = weibull_cdf,
1581 : .sf = weibull_sf,
1582 : .icdf = weibull_icdf,
1583 : .isf = weibull_isf,
1584 : };
1585 :
1586 : /** Functions for generalized Pareto distributions */
1587 :
1588 : static double
1589 700100 : genpareto_sample(const struct dist_t *dist)
1590 : {
1591 700100 : const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1592 700100 : uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng());
1593 700100 : double p0 = random_uniform_01();
1594 :
1595 700100 : return sample_genpareto_locscale(s, p0, GP->mu, GP->sigma, GP->xi);
1596 : }
1597 :
1598 : static double
1599 211 : genpareto_cdf(const struct dist_t *dist, double x)
1600 : {
1601 211 : const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1602 211 : return cdf_genpareto(x, GP->mu, GP->sigma, GP->xi);
1603 : }
1604 :
1605 : static double
1606 1175 : genpareto_sf(const struct dist_t *dist, double x)
1607 : {
1608 1175 : const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1609 1175 : return sf_genpareto(x, GP->mu, GP->sigma, GP->xi);
1610 : }
1611 :
1612 : static double
1613 14 : genpareto_icdf(const struct dist_t *dist, double p)
1614 : {
1615 14 : const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1616 14 : return icdf_genpareto(p, GP->mu, GP->sigma, GP->xi);
1617 : }
1618 :
1619 : static double
1620 7 : genpareto_isf(const struct dist_t *dist, double p)
1621 : {
1622 7 : const struct genpareto_t *GP = dist_to_const_genpareto(dist);
1623 7 : return isf_genpareto(p, GP->mu, GP->sigma, GP->xi);
1624 : }
1625 :
1626 : const struct dist_ops_t genpareto_ops = {
1627 : .name = "generalized Pareto",
1628 : .sample = genpareto_sample,
1629 : .cdf = genpareto_cdf,
1630 : .sf = genpareto_sf,
1631 : .icdf = genpareto_icdf,
1632 : .isf = genpareto_isf,
1633 : };
1634 :
1635 : /** Functions for geometric distribution on number of trials before success */
1636 :
1637 : static double
1638 400100 : geometric_sample(const struct dist_t *dist)
1639 : {
1640 400100 : const struct geometric_t *G = dist_to_const_geometric(dist);
1641 400100 : uint32_t s = crypto_fast_rng_get_u32(get_thread_fast_rng());
1642 400100 : double p0 = random_uniform_01();
1643 :
1644 400100 : return sample_geometric(s, p0, G->p);
1645 : }
1646 :
1647 : static double
1648 0 : geometric_cdf(const struct dist_t *dist, double x)
1649 : {
1650 0 : const struct geometric_t *G = dist_to_const_geometric(dist);
1651 :
1652 0 : if (x < 1)
1653 : return 0;
1654 : /* 1 - (1 - p)^floor(x) = 1 - e^{floor(x) log(1 - p)} */
1655 0 : return -expm1(floor(x)*log1p(-G->p));
1656 : }
1657 :
1658 : static double
1659 0 : geometric_sf(const struct dist_t *dist, double x)
1660 : {
1661 0 : const struct geometric_t *G = dist_to_const_geometric(dist);
1662 :
1663 0 : if (x < 1)
1664 : return 0;
1665 : /* (1 - p)^floor(x) = e^{ceil(x) log(1 - p)} */
1666 0 : return exp(floor(x)*log1p(-G->p));
1667 : }
1668 :
1669 : static double
1670 0 : geometric_icdf(const struct dist_t *dist, double p)
1671 : {
1672 0 : const struct geometric_t *G = dist_to_const_geometric(dist);
1673 :
1674 0 : return log1p(-p)/log1p(-G->p);
1675 : }
1676 :
1677 : static double
1678 0 : geometric_isf(const struct dist_t *dist, double p)
1679 : {
1680 0 : const struct geometric_t *G = dist_to_const_geometric(dist);
1681 :
1682 0 : return log(p)/log1p(-G->p);
1683 : }
1684 :
1685 : const struct dist_ops_t geometric_ops = {
1686 : .name = "geometric (1-based)",
1687 : .sample = geometric_sample,
1688 : .cdf = geometric_cdf,
1689 : .sf = geometric_sf,
1690 : .icdf = geometric_icdf,
1691 : .isf = geometric_isf,
1692 : };
|
