From f73cfcd10e3c0d8992438aa7b868f95eba322ee1 Mon Sep 17 00:00:00 2001 From: Huon Wilson Date: Sun, 27 Oct 2013 21:18:09 +1100 Subject: [PATCH] std::rand: substitute an unbiased Exp1 temporarily. The Ziggurat implementation of Exp1 is incorrect: it appears to have an expected value of around 0.994 rather than 1 and the density of a large number of sample means does not match the normal distribution expected by the central limit theorem (specifically, it's a normal distribution centred at 0.994). This replaces Ziggurat with a "plain" inverse-cdf implementation for now. This implementation has the expected properties, although it is approximately 4 times slower. cc #10084 (this does *not* fix it). --- src/libstd/rand/distributions.rs | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/src/libstd/rand/distributions.rs b/src/libstd/rand/distributions.rs index e7bcf8ce5d3a7..dfaf0e9a96346 100644 --- a/src/libstd/rand/distributions.rs +++ b/src/libstd/rand/distributions.rs @@ -341,20 +341,30 @@ impl IndependentSample for Normal { // has to be unwrapped before use as an `f64` (using either /// `*` or `cast::transmute` is safe). /// -/// Implemented via the ZIGNOR variant[1] of the Ziggurat method. The -/// exact description in the paper was adjusted to use tables for the -/// exponential distribution rather than normal. -/// -/// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to -/// Generate Normal Random -/// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield -/// College, Oxford +/// Implemented via the inverse cdf transform, that is, `-ln(U)` where +/// `U` is uniformly distributed 0 to 1. + +// FIXME #10084: The Ziggurat implementation is biased (mean of 0.994 +// rather than 1). This inverse cdf transform with logs is 4 times +// slower than the (incorrect) Ziggurat implementation. (Note: the +// Normal Ziggurat code appears to be unbiased.) +// +// Old docs: +// Implemented via the ZIGNOR variant[1] of the +// Ziggurat method. The exact description in the paper was adjusted to +// use tables for the exponential distribution rather than normal. +// +// [1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to +// Generate Normal Random +// Samples*](http://www.doornik.com/research/ziggurat.pdf). Nuffield +// College, Oxford pub struct Exp1(f64); // This could be done via `-rng.gen::().ln()` but that is slower. impl Rand for Exp1 { #[inline] fn rand(rng: &mut R) -> Exp1 { + /* FIXME #10084: see above. #[inline] fn pdf(x: f64) -> f64 { (-x).exp() @@ -368,6 +378,11 @@ impl Rand for Exp1 { &ziggurat_tables::ZIG_EXP_X, &ziggurat_tables::ZIG_EXP_F, &ziggurat_tables::ZIG_EXP_F_DIFF, pdf, zero_case)) + */ + + // we use 1 - U, because U is [0, 1), and log(0) is bad, so we + // convert it into (0, 1] + Exp1(-(1.0 - rng.gen::()).ln()) } }