Reconditioning your quantile function

TL;DR

This paper addresses the numerical instability in the inversion method for generating random samples from distributions, introducing a robust reconditioning technique and providing an open-source implementation in pqRand.

Contribution

It proposes a novel reconditioning approach for the quantile function to improve numerical stability in Monte Carlo simulations, with a practical software package available.

Findings

Enhanced accuracy in tail sampling for various distributions

Open-source pqRand package implements the proposed method

Improved robustness of Monte Carlo simulations in nonlinear systems

Abstract

Monte Carlo simulation is an important tool for modeling highly nonlinear systems (like particle colliders and cellular membranes), and random, floating-point numbers are their fuel. These random samples are frequently generated via the inversion method, which harnesses the mapping of the quantile function Q(u) (e.g. to generate proposal variates for rejection sampling). Yet the increasingly large sample size of these simulations makes them vulnerable to a flaw in the inversion method; Q(u) is ill-conditioned in a distribution's tails, stripping precision from its sample. This flaw stems from limitations in machine arithmetic which are often overlooked during implementation (e.g. in popular C++ and Python libraries). This paper introduces a robust inversion method, which reconditions Q(u) by carefully drawing and using uniform variates. pqRand, a free C++ and Python package, implements this novel method for a number of popular distributions (exponential, normal, gamma, and more).

Figures

Equations

$$\{f\}=Q\big{(}\{U(0,1)\}\big{)}\;.$$

$$C(g)\equiv\left|\frac{g(x+\delta x)-g(x)}{g(x)}\Big{/}\frac{\delta x}{x}\right|=\left|x\,\frac{g^{\prime}(x)}{g(x)}\right|+\mathcal{O}(\delta x)\;.$$

$$f(x)=\lambda e^{-\lambda x}⟶$$

$$Q_1(u)=-\frac{1}{\lambda }\log (1-u)=-\frac{1}{\lambda }log1p(-u)⟶$$

$$Q_2(u)=-\frac{1}{\lambda }\log (u)⟶C(Q_2)=-\frac{1}{\log (u)}.$$

$$P^{\ast }(g)\equiv \left|\frac{g(x+\delta x)-g(x)}{g(x)}\right|=\delta x\left|\frac{g^{\prime }(x)}{g(x)}\right|+\mathcal{O}(\delta x^2).$$

$$P_{\text{E}}^{\ast }(Q)=ϵ\left|\frac{Q^{\prime }(u)}{Q(u)}\right|+\mathcal{O}({ϵ}^2).$$

$$H(X)\equiv -\sum _i\mathrm{Pr}(x_i){\log }_2\mathrm{Pr}(x_i).$$

$$H_{\text{E}}=-\sum _{i=1}^nϵ{\log }_2(ϵ)=-{\log }_2ϵ=P.$$

$$H_{\text{N}}=-\sum _{k=1}^K\left(\sum _{i=1}^{n/2}{2}^{-k}(2ϵ){\log }_2\left(2^{-k}(2ϵ)\right)\right)=\sum _{k=1}^K{2}^{-k}(P-1+k)\approx P+1(\text{for}K\gg 1).$$

$$H_{\text{E}}(k)=P-k(\text{for }u<2^{-k}).$$

$$1-\{U_{\text{N}}[0,1)\}\mapsto \{U_{\text{E}}(0,1]\}.$$

$$D_{\text{KL}}(\hat{P}||\hat{Q})=\sum _i\hat{P}(x_i){\log }_2\frac{\hat{P}(x_i)}{\hat{Q}(x_i)}.$$

$$\hat{P}(x_i)=c_i/N$$

$$\hat{Q}(x_i)={\int }_{x_{i,L}}^{x_{i,R}}f(x)\text{d}x=F(x_{i,R})-F(x_{i,L}).$$

Full text

Introduction

The inversion method samples from a probability distribution $f$ via its quantile function ${Q\equiv F^{-1}}$, the inverse of $f$’s cumulative distribution $F$ Devroye (1986); Press et al. (1992). $Q$ is used to transform a random sample from $U(0,1)$, the uniform distribution over the unit interval, into a random sample $\{f\}$;

[TABLE]

This scheme is powerful because quantile functions are formally exact. But any real-world implementation will be formally inexact because: (i) A source of true randomness is generally not practical (or even desirable), while a repeatable pseudo-random number generator (PRNG) is never perfect. (ii) The uniform variates $u$ and their mapping $Q(u)$ use finite-precision machine arithmetic. The first defect has received the lion’s share of attention, leaving the second largely ignored. As a result, common implementations of inversion sampling lose precision in the tails of $f$.

This leak must be subtle if no one has patched it. Nonetheless, the loss of precision commonly exceeds dozens of ULP (units in the last place) in a distribution’s tails. Contrast this to library math functions (sin, exp), which are painstakingly crafted to deliver no more than one ULP of systematic error. When the inversion method loses precision, it produces inferior, repetitive samples, to which Monte Carlo simulations may become sensitive as they grow more complex, drawing ever more random numbers. Proving that the effect is negligible is incredibly difficult, so the best alternative is to use the most numerically stable sampling scheme possible with floating point numbers — if it is not too slow. The robust inversion method proposed here is 80–100% as fast as the original.

To isolate the loss of precision, we examine the three independent steps of inversion sampling:

Generate random bits (i.i.d. coin flips) using a PRNG. 2. 2.

Convert those random bits into a uniform variate $u$ from $U(0,1)$. 3. 3.

Plug $u$ into $Q(u)$ to sample from the distribution $f$.

The first two steps do not depend on $f$, so they are totally generic (a major virtue of the method). Of them, step 1 has been exhaustively studied Knuth (1997); L’Ecuyer (1997); L’Ecuyer and Simard (2007), and is essentially a solved problem — when in doubt, use the Mersenne twister Matsumoto and Nishimura (1998); L’Ecuyer and Simard (2007). Step 3 has been validated using real analysis Devroye (1986); Steinbrecher and Shaw (2008), so that known quantile functions need only be translated into computer math functions.

This leaves step 2 which, at first glance, looks like a trivial coding task to port random bits into a real-valued $Q$. Yet computers cannot use real numbers, and neglecting this fact is dangerous — using this as its central maxim, this paper conducts a careful investigation of the inversion method from step 2 onward. Section I begins by using the condition number to probe step 3, finding that a distribution’s quantile function is numerically unstable in its tails. This provides a sound framework for Sec. II to find the subtle flaw in the canonical algorithm for drawing uniform variates (step 2). A robust inversion method is introduced to fix both problems, and is empirically validated in Sec. III by comparing the near-perfect sample obtained from the pqRand package to the deficient samples obtained from standard C++ and Python tools.

I $\boldsymbol{Q}$ are ill-conditioned, but they do not have to be

Real numbers are not countable, so computers cannot represent them. Machine arithmetic is limited to a countable set like rational numbers $\mathbb{Q}$. The most versatile rational approximation of $\mathbb{R}$ are floating point numbers, or “floats” — scientific notation in base-two $(m\times 2^{E})$. The precision of floats is limited to $P$, the number of binary digits in their mantissa $m$, which forces relative rounding errors of order ${\epsilon\equiv 2^{-P}}$ upon every floating point operation Goldberg (1991). The propagation of such errors makes floating point arithmetic formally inexact. In the worst case, subtle effects like cancellation can degrade the effective (or de facto) precision to just a handful of digits. Using floats with arbitrarily high $P$ mitigates such problems, but is usually emulated in software — an expensive cure. Prudence usually restricts calculations to the largest precision widely supported in hardware, binary64 ($P=53$), commonly called “double” precision.

Limited $P$ makes the intrinsic stability of a computation an important consideration; a result should not change dramatically when its input suffers from a pinch of rounding error. The numerical stability of a function $g(x)$ can be quantified via its condition number $C(g)$ — the relative change in $g(x)$ per the relative change in $x$ Einarsson (2005)

[TABLE]

When an $\mathcal{O}(\epsilon)$ rounding error causes $x$ to increment to the next representable value, $g(x)$ will increment by $C(g)$ representable values. So when $C(g)$ is large (i.e. ${\log_{2}C(g)\to P}$), $g(x)$ is ill-conditioned and imprecise; the tiniest shift in $x$ will cause $g(x)$ to hop over an enormous number of values — values through which the real-valued function passes, and which are representable with floats of precision $P$, but which cannot be attained via the floating point calculation $g(x)$. The condition number should be used to avoid such numerical catastrophes.

We now have a tool to uncover possible instability in the inversion method, specifically in its quantile function $Q$ (step 3). As a case study, we can examine the exponential distribution (the time between events in a Poisson process with rate $\lambda$, like radioactive decay);111${\tt log1p}(x)$ is an implementation of $\log(1+x)$ which sidesteps an unnecessary floating point cancellation ISO (2012).

[TABLE]

A well-conditioned sample from the exponential distribution will require $C(Q_{1})\leq\mathcal{O}(1)$ everywhere, but Fig. 1a clearly reveals that $C(Q_{1})$ (dashed) becomes large as $u\to 1$. Why is $Q_{1}$ ill-conditioned there? According to Eq. 2, a function can become ill-conditioned when it is steep $(\left|g^{\prime}/g\right|\gg 1)$, and $Q_{1}$ (solid) is clearly steep at both $u=0$ and $u=1$. These are $f$’s “tails” — a large range of sample space mapped by a thin, low probability slice of the unit interval. Yet in spite of its steepness, $Q_{1}$ remains well-conditioned throughout its small-value tail ($u\to 0$) because floats are denser near the origin — reusing the same set of mantissae, but with smaller exponents — and a denser set of $u$ allows a more continuous sampling of a rapidly changing $Q_{1}(u)$. This extra density manifests as the singularity-softening factor of $x$ in Eq. 2. Unfortunately, the same relief cannot occur as $u\to 1$, where representable $u$ are not dense enough to accommodate $Q_{1}$’s massive slope.

Because $Q_{1}$ is ill-conditioned near $u=1$, the large-$x$ portion of its sample $\{f\}$ will be imprecise; many large-$x$ floats which should be sampled are skipped-over by $Q_{1}$. This problem is not unique to the exponential distribution; it will occur whenever $f$ has two tails, because one of those tails will be located near $u=1$. Luckily, $U(0,1)$ is perfectly symmetric across the unit interval, so transforming $u\mapsto 1-u$ produces an equally valid quantile function;

[TABLE]

The virtue of using two valid $Q$’s is evident in Fig. 1b; for $u\leq\nicefrac{{1}}{{2}}$, each version is well-conditioned, with $Q_{1}$ sampling the small-value tail $(x\leq\text{median})$ and $Q_{2}$ the large-value tail $(x\geq\text{median})$. Since the pair collectively and stably spans $f$’s entire sample space, $f$ can be sampled via the composition method; for each variate, randomly choose one version of the quantile function (to avoid a high/low pattern), then feed that $Q$ a random $u$ from $U(0,\nicefrac{{1}}{{2}}]$.

This “quantile flip-flop” — a randomized, two-$Q$ composition split at the median — is a simple, general scheme to recondition a quantile function which becomes unstable as $u\to 1$. It is also immediately portable to antithetic variance reduction, a useful technique in Monte Carlo integration where, for every $x=Q(u)$ one also includes the opposite choice $x^{\prime}=Q(u^{\prime})$ Kroese et al. (2013). A common convention is $u^{\prime}\equiv 1-u$, which can create a negative covariance $\text{cov}(x,x^{\prime})$ that decreases the overall variance of the integral estimate. Generating antithetic variates with a quantile flip-flop is trivial; instead of randomly choosing $Q_{1}$ or $Q_{2}$ for each variate, always use both.

II An optimally uniform variate is maximally uneven

The condition number guided the development of the quantile flip-flop, a rather simple way to stabilize step 3 of the inversion method during machine implementation. Our investigation now proceeds to step 2 — sampling uniform variates. While steps 2 and 3 seem independent, we will find that there is an important interplay between them; a quantile function can be destabilized by sub-optimal uniform variates, but it can also wreck itself by mishandling optimal uniform variates.

The canonical method for generating uniform variates is Alg. 1 Knuth (1997); Press et al. (1992); L’Ecuyer (1997); ISO (2012); Python Software Foundation (2017a, b); NumPy Developers (2017a, b); an integer is randomly drawn from $[0,2^{B})$, then scaled to a float in the half-open unit interval $[0,1)$. Using $B\leq P$ produces a completely uniform sample space — each possible $u$ has the same probability, with a rigidly even spacing of $2^{-B}$ between each. Using $B=P$ gives the ultimate even sample $\{U_{\text{E}}[0,1)\}$, as depicted in Fig. 2E (which uses a ridiculously small $B=P=4$ to aide the eye). When $B>P$, line 5 will be forced to round many large $j$, as the mantissa of $a$ is not large enough to store every $j$ with full precision. As $B\to\infty$, this rounding saturates the floats available in $U[0,1)$, creating the uneven $\{U_{\text{N}}[0,1)\}$ depicted Fig. 2N. This uneven sample space is still uniform because large $u$ are more probable, absorbing more $j$ from rounding (due to their coarser spacing).

Depending on the choice of $B$, Alg. 1 can generate uniform variates which are either even or uneven, but which is better? There seem to be no definitive answers in the literature — which is likely why different implementations choose different $B$ — so we will have to find our own answer. We start by choosing the even uniform variate $\{U_{\text{E}}[0,1)\}$ as the null hypothesis, for two obvious reasons: (i) Fig. 2E certainly looks more uniform and (ii) taking $B\to\infty$ does not seem practical. However, we will soon find that perfect evenness has a subtle side effect — it forces all quantile functions to become ill-conditioned as $u\to 0$, even if they have an excellent condition number!

The condition number implicitly assumes that $\delta x$ is vanishingly small. This is true enough for a generic float, whose $\delta x=\mathcal{O}(\epsilon\,x)$ is much small than $x$. But the even uniform variates have an absolute spacing of ${\delta u=\epsilon}$. To account for a finite $\delta x$, we define a function’s effective precision

[TABLE]

Like $C(g)$, a large effective precision $P^{*}(g)$ indicates an ill-conditioned calculation. For a generic floating point calculation $\delta x=\mathcal{O}(\epsilon\,x)$, so $P^{*}$ reverts back to the condition number (${P^{*}(g)\approx\epsilon\,C(g)}$). But feeding even uniform variates into a quantile function gives $\delta u=\epsilon$, so

[TABLE]

Calculating $P^{*}_{\text{E}}(Q)$ for the quantile flip-flop of Fig. 1b indicates that both $Q$ become ill-conditioned as $u\to 0$ (where $Q$ becomes steep), in stark opposition to their excellent condition numbers. That using even uniform variates will break a quantile flip-flop is a problem not unique to the exponential distribution; it occurs whenever $f$ has a tail (so that $\left|Q^{\prime}/Q\right|\to\infty$ as $u\to 0$).

The reduced effective precision $P^{*}_{\text{E}}(Q)$ caused by even uniform variates creates sparsely populated tails; there are many extreme values which $\{f\}$ will never contain, and those which it does will be sampled too often. $\{U_{\text{E}}[0,1)\}$ is simply too finite; $2^{P}$ even uniform variates can supply no more than $2^{P}$ unique values. This implies that the uneven sample $\{U_{\text{N}}[0,1)\}$ will restore quantile stability, since its denser input space ($\delta u=\mathcal{O}(\epsilon u)$) will stabilize $P^{*}_{\text{N}}(Q)$ near the origin. These small $u$ expand the sample space of $\{f\}$ many times over, making its tails far less repetitive. And since uneven variates correspond to the limit where $B\to\infty$ in Alg. 1, they are equivalent to sampling $U[0,1-\epsilon)$ from $\mathbb{R}$ and rounding to the nearest float — the next best thing to a real-valued input for $Q$.

The virtue of using uneven uniform variates also follows from information theory. The Shannon entropy of a sample space $X$ counts how many bits of information are conveyed by each variate $x$;

[TABLE]

The sample space of the even uniform variates ($B=P$) has $n=\epsilon^{-1}$ equiprobable members, so

[TABLE]

This makes sense, since each even uniform variate originates from a $P$-bit pseudo-random integer.

The sample space of the uneven $\{U_{\text{N}}[0,1)\}$ contains every float in $[0,1)$, which is naturally partitioned into sub-domains $[2^{-k},2^{-k+1})$ with common exponent $-k$. Each domain comprises a fraction $2^{-k}$ of the unit interval, and the minimum exponent $-K$ depends on the floating point type (although $K\gg 1$ for binary32 and binary64). The uneven entropy is then the sum over sub-domains, each of which sums over the $n/2$ equiprobable mantissae222Ignoring the fact that exact powers of 2 are $\nicefrac{{3}}{{4}}$ as probable, which makes no difference once $P\gtrsim 10$.

[TABLE]

One more bit of information than even variates is not a windfall. But $H_{\text{E}}$ and $H_{\text{N}}$ are the entropies of the bulk sample $\{U[0,1)\}$. What is the entropy of the tail-sampling sub-space $U[0,2^{-k})$?

Rejecting all $u\geq 2^{-k}$ in the even sample $\{U_{\text{E}}[0,1)\}$, we find that smaller $u$ have less information

[TABLE]

This lack of information in even variates is inevitably mapped to the sample $\{f\}$, consistent with the deteriorating effective precision as $u\to 0$. But for uneven uniform variates, the sample space is fractal; each sub-space looks the same as the whole unit interval, so that ${H_{\text{N}}(k)=P+1}$ as before! Every $u$ has maximal information, and a high-entropy input should give a high-precision sample.

Both the effective precision $P^{*}(Q)$ and Shannon entropy $H$ predict that using even uniform variates will force a well-conditioned quantile function to become ill-conditioned, precluding a high-precision sample. Switching to uneven uniform variates will recondition it. But there is an important caveat; uneven variates are very delicate. Subtracting them from one mutates them back into even variates (with opposite boundary conditions);

[TABLE]

This is floating point cancellation. The subtraction erases any extra density in the uneven sample, because it maps the very dense region (near zero) to a region where floats are intrinsically sparse (near one). Conversely, the sparse region of the uneven sample (near one) has no extra information to convey when it is mapped near zero, and remains sparse. This is why $Q_{1}$ (Eq. 4) must use log1p.

III Precision: lost and found

In Sec. I we conditioned an intrinsically imprecise quantile function using a two-$Q$ composition. Then in Sec. II we determined that uneven uniform variates are required to keep $Q$ well-conditioned. These two practices comprise our robust inversion method, whose technical details we have deliberately left for the Appendix because we have yet to prove that it makes a material difference. If indiscreet sampling decimates the precision of $\{f\}$, it should be quite evident in an experiment!

The quality of a real-world sample $\{f\}$ can be assessed via its Kullback-Leibler divergence Ben-David et al. (2015)

[TABLE]

$D_{\text{KL}}$ quantifies the relative entropy between a posterior distribution $\widehat{P}$ and a prior distribution $\widehat{Q}$ (c.f. Eq. 8). The empirical $\widehat{P}$ is based on the count $c_{i}$ — the number of times $x_{i}$ appears in $\{f\}$

[TABLE]

(where $N$ is the sample size). The ideal density $\widehat{Q}$ is obtained by mapping $f$ onto floats, using the domain of real numbers $(x_{i,L},x_{i,R})$ that round to each $x_{i}$;

[TABLE]

$D_{\text{KL}}$ does not sum terms where $\widehat{P}(x_{i})=0$ (i.e. $x_{i}$ was not drawn), because ${\lim_{x\to 0}x\log x=0}$.

The $D_{\text{KL}}$ divergence is not a metric because it is not symmetric under exchange of $\widehat{P}$ and $\widehat{Q}$ Ben-David et al. (2015). And while $D_{\text{KL}}$ is frequently interpreted as the information gained when using distribution $\widehat{P}$ instead of $\widehat{Q}$, this is not true here. Consider a PRNG which samples from $\widehat{Q}=U(0,1)$, but samples so poorly that it always outputs $x=0.5$ (and thus emits zero information). Its $D_{\text{KL}}\approx P$ is clearly the precision lost by $\widehat{P}$ (the generator). In less extreme cases, since $\widehat{Q}$ is the most precise distribution possible given floats of precision $P$, any divergence denotes how many bits of precision were lost.

Our experiments calculate $D_{KL}$ for samples of the $\lambda=1$ exponential distribution generated via the inversion method. We use GNU’s std::mt19937 for our PRNG ($B=32$), fully seeding its state from the computer’s environmental noise (using GNU’s std::random_device). Calculating $D_{\text{KL}}$ requires recording the count for each unique float, and an accurate $D_{\text{KL}}$ requires a very large sample size ($N\gg P$, so that $\widehat{P}\to\widehat{Q}$ in the case of perfect agreement). To keep the experiments both exhaustive and tractable, and with no loss of generality, we use binary32 ($P=24$, or single precision). Since double precision is governed by the same IEEE 754 standard IEEE Task P754 (2008), and both types use library math functions with $\mathcal{O}(\epsilon)$ errors, the $D_{\text{KL}}$ results for binary64 will be identical.333A binary64 experiment is tractable, just not exhaustive. Memory constraints require intricate simulation of tiny sub-spaces of the unit interval, to act as a representative sample of the whole.

The first implementation we test is GNU’s std::exponential_distribution, a member of the C++11 <random> suite, which obtains its uniform variates from std::generate_canonical Free Software Foundation (2017a, b). Given our PRNG, these uniform variates are equivalent to calling Alg. 1 with ${B=32}$ and ${P=24}$. This creates a partially uneven sample $\{U_{\text{P}}[0,1)\}$, with ${B-P=8}$ bits more entropy than even variates. GNU’s implementation feeds these uniform variates into $Q_{1}$ (Eq. 4), but without removing its cancellation by using log1p. As predicted by Eq. 12, the cancellation strips any extra entropy from the partially uneven variates ($B>P$), converting then into even ones ($B=P$).

Figure 3 shows the bits of precision lost by three samples using the same PRNG seed, with the median at the center and increasingly improbable values near the edges — a format which becomes easier to understand by referring to the top axis, which shows the $x=Q(u)$ sampled by the various $u$. GNU’s std::exponential_distribution (▼) exhibits a clear and dramatic loss of precision as variates gets farther from the median (and more rare). This imprecision agrees exactly with the prediction of $H_{\text{E}}(k)$ (Eq. 11, solid line) — even uniform variates have limited information, and every time $u$ becomes half as small (so that $x$ is half as probable), one more bit of precision is lost. This loss of precision in the sample is clearly caused by using uniform variates, which will always happen if $Q_{1}$ neglects to use log1p internally. Since both Python’s random.expovariate Python Software Foundation (2017a) and Numpy’s numpy.random.exponential NumPy Developers (2017a) also commit this error, their samples are equally imprecise.

But GNU’s std::exponential_distribution could have done better; it drew partially uneven uniform variates ($B=32$, $P=24$), then spoiled them via cancellation. Enabling log1p in $Q_{1}$ and regenerating GNU’s sample (▲) permits Fig. 3 to isolate the two sources of imprecision identified in Secs. I & II. (i) Using log1p, $Q_{1}$ is allowed to be well-conditioned as $u\to 0$, so only the uniform variates themselves can degrade the small-value tail. Moving left from the median, the partially uneven variates maintain maximal precision until their 8-bit entropy buffer runs dry. (ii) Conversely, $Q_{1}$ is intrinsically ill-conditioned for $u>\nicefrac{{1}}{{2}}$ in the large-value tail, so the quality of the uniform variates is irrelevant; an ill-conditioned quantile function causes an immediate loss of precision.

pqRand generates its sample (❍) via our robust inversion method, feeding high-entropy, uneven uniform variates $U_{\text{N}}(0,\nicefrac{{1}}{{2}}]$ into a quantile flip-flop which is always well-conditioned ($Q_{1}$ samples ${x\leq\text{median}}$ and $Q_{2}$ samples ${x\geq\text{median}}$). Switching to a quantile flip-flop for this final data series means that, to the right of the median, the small values shown on the bottom horizontal axis are now $u$ instead of $1-u$. The sample’s tails exhibit ideal performance, in stark contrast to the standard inversion method, and precision is only lost near the median, where the composite $Q$ is a tad unstable ($C(Q)\gtrsim 1$, see Fig. 1b). That $D_{\text{KL}}\approx 0$ everywhere, and never exceeds 1 bit, is clear evidence that our robust inversion method fulfills its existential purpose, delivering the best sample possible with floats of precision $P$. Furthermore, this massive boost in quality arrives at ${\sim}80/100\%$ the speed of GNU’s std::exponential_distribution for binary32/64 (${\sim}30/40\,\mathrm{ns}$ per variate on an Intel i7 @ 2.9 GHz with GCC 6.3, optimization O2).

Similar samples for any rate $\lambda$, as well as many other distributions (uniform, normal, log-normal, Weibull, logistic, gamma) are available with pqRand, a free C++ and Python package hosted on GitHub Pedersen (2017). pqRand uses optimized C++ to generate uneven uniform variates (see the Appendix), with Cython wrappers for fast scripting. Yet the usefulness of pqRand is not restricted to the rarefied set of distributions with analytic quantile functions; pqRand uses rejection sampling for its own normal and gamma distributions. Rejection sampling gives access to any distribution $f(x)$, provided that one can more easily sample from the proposal distribution $g(x)\geq f(x)$. Since the final sample $\{f\}$ is merely a subset of the proposed sample $\{g\}$, a high-precision $\{f\}$ requires a high-precision $\{g\}$, which can be obtained via our robust inversion method.

IV Conclusion

Using the exponential distribution as a case study, we find two general sources of imprecision when sampling a probability distribution $f$ via the inversion method: (i) When $f$ has two tails (two places where ${Q^{\prime}/Q\gg 1}$), its quantile function $Q(u)$ becomes ill-conditioned as $u\to 1$. (ii) Drawing uniform random variates using the canonical algorithm (Alg. 1) gives too finite a sample space, making $Q(u)$ ill-conditioned as $u\to 0$ (even if $Q(u)$ has a good condition number there). Both problems can lose dozens of ULP of precision in a sample’s tails, and they are especially problematic for simulations using single precision — in the worst case, ${\sim 0.5\%}$ of variates will lose at least a third of their precision. This vulnerability is found in popular implementations of the inversion method (e.g. GNU’s implementation of C++11’s <random> suite Free Software Foundation (2017a), and the python.random Python Software Foundation (2017a) and numpy.random NumPy Developers (2017a) modules for Python, and more).

This paper introduces a robust inversion method which reconditions $Q$ by combining (i) uneven uniform variates (Alg. 2, see Appendix) with (ii) a quantile flip-flop (a two-$Q$ composition split at the median). Our method produces the best sample from $f$ possible with floats of precision $P$, and is significantly faster than schemes which “exactly” sample distributions to arbitrary precision Knuth and Yao (1976); Flajolet et al. (2011); Karney (2016). The precision of a random sample is especially important for large, non-linear Monte Carlo simulations, which can draw so many numbers that they may be sensitive to this vulnerability. Since it is generally difficult to exhaustively validate large simulations — in this case, to prove that a loss of precision in the tails has only negligible effects — the best strategy is to use the most numerically stable components at every step in the simulation chain, provided they are not prohibitively slow. To this end, we have released pqRand Pedersen (2017), a free C++ and Python implementation of our robust inversion method, which is 80–100% as fast as standard inversion sampling.

V Acknowledgements

Thanks to Zack Sullivan for his invaluable help and editorial suggestions, and to Andrew Webster for lowering the activation energy. This work was supported by the U.S. Department of Energy under award No. DE-SC0008347 and by Validate Health LLC.

Appendix A Drawing uneven uniform variates

In Sec. II we saw that the best uniform variates are uneven, obtained by taking $B\to\infty$ in Alg. 1. Since this will take an infinite amount of time, we must devlop an alternate scheme. A clue lies in the bitwise representation of the even uniform variate from Alg. 1, for which every $u<2^{-k}$ has a reduced entropy $H_{\text{E}}=P-k$ (Eq. 11). When $B=P$, Alg. 1 draws an integer $M$ from $[0,2^{P})$, then converts it to floating point. Inside the resulting float, the mantissa is stored as the integer $M^{*}$, which is just the original integer $M$ with its bits shifted left until $M^{*}\geq 2^{P-1}$. This bit-shift ensures that any $u<2^{-k}$ always has at least $k$ trailing zeroes in $M^{*}$; zeroes which contain no information. Filling this always-zero hole with new random bits will restore maximal entropy.

Given the domain required by a quantile flip-flop, Alg. 2 samples uneven $\{U_{\text{N}}(0,\nicefrac{{1}}{{2}}]\}$ from the half-open, half-unit interval. It works by taking $B\to\infty$, yet knowing that floating point arithmetic will truncate the infinite bit-stream to $P$ bits of precision. As soon as the RNG returns the first 1 (however many bits that takes), only the next $P+1$ bits are needed to convert to floating point; $P$ bits to fill the mantissa, and two extra bits for proper rounding. To fix $u$’s coarse location, the first loop (line 6) finds the first significant bit. The following conditional (line 7) requires $S\geq P+2$ significant bits. If $S$ is too small, $j$’s bits are shifted left until the most significant (leftmost) bit slides into the $P+2$ position (line 12). Then the vacated space on the right is filled with new random bits, and the leftward shift is factored into $n$, so that only $u$’s fine location changes (enhancing precision while preserving uniformity). Finally, the integer is rounded into $(0,\nicefrac{{1}}{{2}}]$.444We exclude zero from the output domain of Alg. 2 because, while theoretically possible, it will never happen (given a reliable RNG). Returning zero in binary32 (single precision) would require drawing more than 150 all-zero bits in the first loop. Given a billion cores drawing $B=32$ every nanosecond, that would take $\mathcal{O}(10^{55})\,\mathrm{years}$ (although the first variate with sub-maximal entropy would only take $\mathcal{O}(10^{41})\,\mathrm{years}$). For binary64, the numbers get ridiculous.

Algorithm 2 needs two extra bits to maintain uniformity when $j$ is converted to a float. With few exceptions, exact conversion of integers larger than $2^{P}$ is not possible because the mantissa lacks the necessary precision. Truncation $j$ won’t work because $j<2^{n-1}$, so Alg. 2 would never return $u=\nicefrac{{1}}{{2}}$, a value needed by a quantile flip-flop to sample the exact median. Since Alg. 2 must be able to round $j$ up, it uses round-to-nearest, ties-to-even (R2N-T2E). Being the most numerically stable IEEE 754 rounding mode, R2N-T2E is the default choice for most operating systems.

Yet R2N-T2E is slightly problematic because Alg. 2 is truncating a theoretically infinite bit stream to finite significance $S$. There are going to be rounding ties, and when T2E kicks in, it will pick even mantissae over odd ones, breaking uniformity. To defeat this bias, $j$ is made odd. This creates a systematic tie-breaker, because an odd $j$ is always closer to only one of the truncated options, without giving preference to the even option. This system only fails when $S=P+1$, and only the final bit needs removal. In this case, $j$ is equidistant from the two options, and T2E kicks in. Adding a random buffer bit (requiring $S\geq P+2$) precludes this failure.

An important property of Alg. 2 is that $u=\nicefrac{{1}}{{2}}$ is half as probable as its neighbor, $u=\frac{1}{2}(1-\epsilon)$. Imagine dividing the domain $[\nicefrac{{1}}{{4}},\nicefrac{{1}}{{2}}]$ into $2^{P-1}$ bins, with the bin edges depicting the representable $u$ in that domain. Uniformly filling the domain with $\mathbb{R}$, each $u$ absorbs a full bin of real numbers via rounding (a half bin to its left, a half bin to its right). The only exception is $u=\nicefrac{{1}}{{2}}$, which can only absorb a half bin from the left, making it half as probable. But recall that $\{U_{\text{N}}(0,\nicefrac{{1}}{{2}}]\}$ is intended for use in a quantile flip-flop — a regular quantile function folded in half at the median ($u=\nicefrac{{1}}{{2}}$). Since both $Q$ map to the median when they are fed $u=\nicefrac{{1}}{{2}}$, the median will be double-counted unless $u=\nicefrac{{1}}{{2}}$ is half as probable.

Not only can Alg. 2 produce better uniform variates than std::generate_canonical (see Fig. 3), it does so at equivalent computational speed (${\sim}5\,\mathrm{ns}$ per variate using MT19937 on an Intel i7 @ 2.9 GHz). This is possible because line 7 is rarely true (${\sim}0.1\%$ when $N=64$ and $P=53$), so the code to top-up entropy is rarely needed, and the main conditional branch is quite predictable. For most variates, the only extra overhead is verifying that $S\geq P+2$, then making $j$ odd, which take no time compared to the RNG and R2N-T2E operations.