|
1 | 1 | //! Sequential Importance Sampling (SIS) primitives — multinomial |
2 | | -//! sampling with replacement, in $O(n)$ time. |
| 2 | +//! sampling with replacement, in O(n) time. |
3 | 3 | //! |
4 | 4 | //! See `README.md` for a tutorial introduction and `INTERNALS.md` for |
5 | 5 | //! the algorithm specification, math proofs, and floating-point |
|
8 | 8 | //! # API at a glance |
9 | 9 | //! |
10 | 10 | //! - [`sample_indices`] — the main entry point. Returns an iterator |
11 | | -//! yielding $n$ indices into `weights` iid with replacement, each |
| 11 | +//! yielding `n` indices into `weights` iid with replacement, each |
12 | 12 | //! with probability proportional to its weight. Output is in |
13 | 13 | //! ascending order. Streaming: one `powf` call per yielded index. |
14 | 14 | //! - [`sample_indices_buffered`] — buffered variant, typically |
15 | 15 | //! ~1.32× faster on x86 (more on hardware with a slow `powf`). |
16 | 16 | //! Takes an `&mut [u32]` buffer rather than returning an iterator |
17 | 17 | //! (it uses the buffer as f32 scratch). |
18 | | -//! - [`SortedUniforms`] — iterator yielding $n$ $\mathrm{Uniform}(0, 1)$ |
19 | | -//! variates in ascending order in $O(n)$ time. Useful in its own |
20 | | -//! right outside sampling (e.g. inverse-CDF sampling where you |
21 | | -//! want sorted output). |
| 18 | +//! - [`SortedUniforms`] — iterator yielding `n` Uniform(0, 1) variates |
| 19 | +//! in ascending order in O(n) time. Useful in its own right outside |
| 20 | +//! sampling (e.g. inverse-CDF sampling where you want sorted |
| 21 | +//! output). |
22 | 22 | //! - [`first_uniform`] — low-level per-step primitive used by |
23 | | -//! [`SortedUniforms`]. Samples $\min(U_1, \ldots, U_k)$ for $k$ |
24 | | -//! iid $\mathrm{Uniform}(0, 1)$ (equivalently, $\mathrm{Beta}(1, k)$). |
25 | | -//! Most callers won't touch this directly. |
| 23 | +//! [`SortedUniforms`]. Samples the minimum of `k` iid Uniform(0, 1) |
| 24 | +//! draws (≡ `Beta(1, k)`). Most callers won't touch this directly. |
26 | 25 | //! |
27 | 26 | //! Indices are yielded/written as `u32`, not `usize`, so the API has |
28 | 27 | //! the same layout on every platform. Callers cast to `usize` at the |
@@ -69,13 +68,12 @@ use rand_distr::{Distribution, Exp1}; |
69 | 68 | // first_uniform |
70 | 69 | // --------------------------------------------------------------------------- |
71 | 70 |
|
72 | | -/// Draw a sample distributed as $\min(U_1, \ldots, U_k)$ for $k$ |
73 | | -/// iid $\mathrm{Uniform}(0, 1)$ variates (equivalently, |
74 | | -/// $\mathrm{Beta}(1, k)$ in standard notation). |
| 71 | +/// Draw a sample distributed as the minimum of `k` iid Uniform(0, 1) |
| 72 | +/// variates (equivalently, `Beta(1, k)` in standard notation). |
75 | 73 | /// |
76 | 74 | /// Computed in closed form via the inverse CDF |
77 | | -/// $F^{-1}(u) = 1 - (1 - u)^{1/k}$. One `powf` call per invocation; |
78 | | -/// constant time in $k$. |
| 75 | +/// `F⁻¹(u) = 1 − (1 − u)^(1/k)`. One `powf` call per invocation; |
| 76 | +/// constant time in `k`. |
79 | 77 | /// |
80 | 78 | /// This is the per-step primitive driving the order-statistic |
81 | 79 | /// recurrence in [`SortedUniforms`]. Most callers won't need this |
@@ -104,18 +102,17 @@ pub fn first_uniform<R: Rng + ?Sized>(rng: &mut R, k: u32) -> f32 { |
104 | 102 | // SortedUniforms |
105 | 103 | // --------------------------------------------------------------------------- |
106 | 104 |
|
107 | | -/// Streaming iterator yielding $n$ $\mathrm{Uniform}(0, 1)$ variates |
108 | | -/// in ascending order in $O(n)$ time. |
| 105 | +/// Streaming iterator yielding `n` Uniform(0, 1) variates in ascending |
| 106 | +/// order in O(n) time. |
109 | 107 | /// |
110 | 108 | /// The yielded values are distributed exactly as the order statistics |
111 | | -/// $U_{(1)} \le U_{(2)} \le \cdots \le U_{(n)}$ of $n$ iid |
112 | | -/// $\mathrm{Uniform}(0, 1)$ draws — i.e. the same as drawing $n$ iid |
| 109 | +/// of `n` iid Uniform(0, 1) draws — i.e. the same as drawing `n` iid |
113 | 110 | /// uniforms and sorting them, but produced one at a time without a |
114 | 111 | /// sort. Internally uses the Bentley–Saxe spacings recurrence |
115 | | -/// $U_{(i)} = U_{(i-1)} + (1 - U_{(i-1)}) \cdot Z$ with |
116 | | -/// $Z \sim \mathrm{Beta}(1, n - i + 1)$ supplied by [`first_uniform`]. |
| 112 | +/// `U_(i) = U_(i-1) + (1 − U_(i-1)) · Z` with `Z ~ Beta(1, n − i + 1)` |
| 113 | +/// supplied by [`first_uniform`]. |
117 | 114 | /// |
118 | | -/// Holds a mutable reference to the RNG. Yields exactly $n$ values, |
| 115 | +/// Holds a mutable reference to the RNG. Yields exactly `n` values, |
119 | 116 | /// then `None` thereafter. |
120 | 117 | /// |
121 | 118 | /// See `INTERNALS.md` §3.1 / §5.1 for the algorithm and its |
@@ -193,15 +190,15 @@ fn kahan_add(sum: &mut f32, c: &mut f32, x: f32) { |
193 | 190 | *sum = t; |
194 | 191 | } |
195 | 192 |
|
196 | | -/// Build an iterator yielding $n$ indices into `weights`, each |
| 193 | +/// Build an iterator yielding `n` indices into `weights`, each |
197 | 194 | /// drawn iid (with replacement) with probability proportional to |
198 | 195 | /// its weight ("multinomial sampling"). Indices are yielded in |
199 | 196 | /// ascending order; collect into a `Vec` and shuffle afterward if |
200 | 197 | /// you need them in random order. |
201 | 198 | /// |
202 | | -/// Streaming variant: runs in $O(m + n)$ total (where |
203 | | -/// $m$ = `weights.len()`), allocates nothing, and uses one |
204 | | -/// [`first_uniform`] call (one `powf`) per yielded index. |
| 199 | +/// Streaming variant: runs in O(`weights.len()` + `n`) total, |
| 200 | +/// allocates nothing, and uses one [`first_uniform`] call (one |
| 201 | +/// `powf`) per yielded index. |
205 | 202 | /// |
206 | 203 | /// # Preconditions |
207 | 204 | /// - `weights` is nonempty. |
@@ -307,16 +304,12 @@ impl<'a, R: Rng + ?Sized> FusedIterator for SampleIndices<'a, R> {} |
307 | 304 | /// hardware with a slow `powf`). |
308 | 305 | /// |
309 | 306 | /// Generates sorted uniforms via the Gamma-ratio identity |
310 | | -/// |
311 | | -/// $$ |
312 | | -/// U_{(i)} = \frac{E_1 + \cdots + E_i}{E_1 + \cdots + E_{n+1}}, |
313 | | -/// \qquad E_j \sim \mathrm{Exp}(1) \text{ iid} |
314 | | -/// $$ |
315 | | -/// |
316 | | -/// rather than via [`first_uniform`], avoiding the per-element |
317 | | -/// `powf`. Internally repurposes `out` as scratch (each `u32` slot |
318 | | -/// temporarily holds the f32 bit pattern of an $E_j$ draw via |
319 | | -/// [`f32::to_bits`], later overwritten with the output index). |
| 307 | +/// (`U_(i) = (E_1 + ... + E_i) / (E_1 + ... + E_(n+1))` for `E_j` |
| 308 | +/// iid Exp(1)) rather than via [`first_uniform`], avoiding the |
| 309 | +/// per-element `powf`. Internally repurposes `out` as scratch |
| 310 | +/// (each `u32` slot temporarily holds the f32 bit pattern of an |
| 311 | +/// `E_j` draw via [`f32::to_bits`], later overwritten with the |
| 312 | +/// output index). |
320 | 313 | /// |
321 | 314 | /// See `INTERNALS.md` §4.4 for the algorithm and §5.2 for the |
322 | 315 | /// `target.min(total)` clip that keeps the merge bounded. |
|
0 commit comments