1 // Copyright 2018 Developers of the Rand project.
2 // Copyright 2013 The Rust Project Developers.
4 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5 // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6 // <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7 // option. This file may not be copied, modified, or distributed
8 // except according to those terms.
10 //! The normal and derived distributions.
13 use distributions
::{ziggurat_tables, Distribution, Open01}
;
14 use distributions
::utils
::ziggurat
;
16 /// Samples floating-point numbers according to the normal distribution
17 /// `N(0, 1)` (a.k.a. a standard normal, or Gaussian). This is equivalent to
18 /// `Normal::new(0.0, 1.0)` but faster.
20 /// See `Normal` for the general normal distribution.
22 /// Implemented via the ZIGNOR variant[^1] of the Ziggurat method.
24 /// [^1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to
25 /// Generate Normal Random Samples*](
26 /// https://www.doornik.com/research/ziggurat.pdf).
27 /// Nuffield College, Oxford
31 /// use rand::prelude::*;
32 /// use rand::distributions::StandardNormal;
34 /// let val: f64 = SmallRng::from_entropy().sample(StandardNormal);
35 /// println!("{}", val);
37 #[derive(Clone, Copy, Debug)]
38 pub struct StandardNormal
;
40 impl Distribution
<f64> for StandardNormal
{
41 fn sample
<R
: Rng
+ ?Sized
>(&self, rng
: &mut R
) -> f64 {
43 fn pdf(x
: f64) -> f64 {
47 fn zero_case
<R
: Rng
+ ?Sized
>(rng
: &mut R
, u
: f64) -> f64 {
48 // compute a random number in the tail by hand
50 // strange initial conditions, because the loop is not
51 // do-while, so the condition should be true on the first
52 // run, they get overwritten anyway (0 < 1, so these are
57 while -2.0 * y
< x
* x
{
58 let x_
: f64 = rng
.sample(Open01
);
59 let y_
: f64 = rng
.sample(Open01
);
61 x
= x_
.ln() / ziggurat_tables
::ZIG_NORM_R
;
65 if u
< 0.0 { x - ziggurat_tables::ZIG_NORM_R }
else { ziggurat_tables::ZIG_NORM_R - x }
68 ziggurat(rng
, true, // this is symmetric
69 &ziggurat_tables
::ZIG_NORM_X
,
70 &ziggurat_tables
::ZIG_NORM_F
,
75 /// The normal distribution `N(mean, std_dev**2)`.
77 /// This uses the ZIGNOR variant of the Ziggurat method, see [`StandardNormal`]
80 /// Note that [`StandardNormal`] is an optimised implementation for mean 0, and
81 /// standard deviation 1.
86 /// use rand::distributions::{Normal, Distribution};
88 /// // mean 2, standard deviation 3
89 /// let normal = Normal::new(2.0, 3.0);
90 /// let v = normal.sample(&mut rand::thread_rng());
91 /// println!("{} is from a N(2, 9) distribution", v)
94 /// [`StandardNormal`]: crate::distributions::StandardNormal
95 #[derive(Clone, Copy, Debug)]
102 /// Construct a new `Normal` distribution with the given mean and
103 /// standard deviation.
107 /// Panics if `std_dev < 0`.
109 pub fn new(mean
: f64, std_dev
: f64) -> Normal
{
110 assert
!(std_dev
>= 0.0, "Normal::new called with `std_dev` < 0");
117 impl Distribution
<f64> for Normal
{
118 fn sample
<R
: Rng
+ ?Sized
>(&self, rng
: &mut R
) -> f64 {
119 let n
= rng
.sample(StandardNormal
);
120 self.mean
+ self.std_dev
* n
125 /// The log-normal distribution `ln N(mean, std_dev**2)`.
127 /// If `X` is log-normal distributed, then `ln(X)` is `N(mean, std_dev**2)`
133 /// use rand::distributions::{LogNormal, Distribution};
135 /// // mean 2, standard deviation 3
136 /// let log_normal = LogNormal::new(2.0, 3.0);
137 /// let v = log_normal.sample(&mut rand::thread_rng());
138 /// println!("{} is from an ln N(2, 9) distribution", v)
140 #[derive(Clone, Copy, Debug)]
141 pub struct LogNormal
{
146 /// Construct a new `LogNormal` distribution with the given mean
147 /// and standard deviation.
151 /// Panics if `std_dev < 0`.
153 pub fn new(mean
: f64, std_dev
: f64) -> LogNormal
{
154 assert
!(std_dev
>= 0.0, "LogNormal::new called with `std_dev` < 0");
155 LogNormal { norm: Normal::new(mean, std_dev) }
158 impl Distribution
<f64> for LogNormal
{
159 fn sample
<R
: Rng
+ ?Sized
>(&self, rng
: &mut R
) -> f64 {
160 self.norm
.sample(rng
).exp()
166 use distributions
::Distribution
;
167 use super::{Normal, LogNormal}
;
171 let norm
= Normal
::new(10.0, 10.0);
172 let mut rng
= ::test
::rng(210);
174 norm
.sample(&mut rng
);
179 fn test_normal_invalid_sd() {
180 Normal
::new(10.0, -1.0);
185 fn test_log_normal() {
186 let lnorm
= LogNormal
::new(10.0, 10.0);
187 let mut rng
= ::test
::rng(211);
189 lnorm
.sample(&mut rng
);
194 fn test_log_normal_invalid_sd() {
195 LogNormal
::new(10.0, -1.0);