1// Copyright 2018 Developers of the Rand project.
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
6// option. This file may not be copied, modified, or distributed
7// except according to those terms.
8
9//! The Bernoulli distribution `Bernoulli(p)`.
10
11use crate::distr::Distribution;
12use crate::Rng;
13use core::fmt;
14
15#[cfg(feature = "serde")]
16use serde::{Deserialize, Serialize};
17
18/// The [Bernoulli distribution](https://en.wikipedia.org/wiki/Bernoulli_distribution) `Bernoulli(p)`.
19///
20/// This distribution describes a single boolean random variable, which is true
21/// with probability `p` and false with probability `1 - p`.
22/// It is a special case of the Binomial distribution with `n = 1`.
23///
24/// # Plot
25///
26/// The following plot shows the Bernoulli distribution with `p = 0.1`,
27/// `p = 0.5`, and `p = 0.9`.
28///
29/// ![Bernoulli distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/bernoulli.svg)
30///
31/// # Example
32///
33/// ```rust
34/// use rand::distr::{Bernoulli, Distribution};
35///
36/// let d = Bernoulli::new(0.3).unwrap();
37/// let v = d.sample(&mut rand::rng());
38/// println!("{} is from a Bernoulli distribution", v);
39/// ```
40///
41/// # Precision
42///
43/// This `Bernoulli` distribution uses 64 bits from the RNG (a `u64`),
44/// so only probabilities that are multiples of 2<sup>-64</sup> can be
45/// represented.
46#[derive(Clone, Copy, Debug, PartialEq)]
47#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
48pub struct Bernoulli {
49 /// Probability of success, relative to the maximal integer.
50 p_int: u64,
51}
52
53// To sample from the Bernoulli distribution we use a method that compares a
54// random `u64` value `v < (p * 2^64)`.
55//
56// If `p == 1.0`, the integer `v` to compare against can not represented as a
57// `u64`. We manually set it to `u64::MAX` instead (2^64 - 1 instead of 2^64).
58// Note that value of `p < 1.0` can never result in `u64::MAX`, because an
59// `f64` only has 53 bits of precision, and the next largest value of `p` will
60// result in `2^64 - 2048`.
61//
62// Also there is a 100% theoretical concern: if someone consistently wants to
63// generate `true` using the Bernoulli distribution (i.e. by using a probability
64// of `1.0`), just using `u64::MAX` is not enough. On average it would return
65// false once every 2^64 iterations. Some people apparently care about this
66// case.
67//
68// That is why we special-case `u64::MAX` to always return `true`, without using
69// the RNG, and pay the performance price for all uses that *are* reasonable.
70// Luckily, if `new()` and `sample` are close, the compiler can optimize out the
71// extra check.
72const ALWAYS_TRUE: u64 = u64::MAX;
73
74// This is just `2.0.powi(64)`, but written this way because it is not available
75// in `no_std` mode.
76const SCALE: f64 = 2.0 * (1u64 << 63) as f64;
77
78/// Error type returned from [`Bernoulli::new`].
79#[derive(Clone, Copy, Debug, PartialEq, Eq)]
80pub enum BernoulliError {
81 /// `p < 0` or `p > 1`.
82 InvalidProbability,
83}
84
85impl fmt::Display for BernoulliError {
86 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
87 f.write_str(match self {
88 BernoulliError::InvalidProbability => "p is outside [0, 1] in Bernoulli distribution",
89 })
90 }
91}
92
93#[cfg(feature = "std")]
94impl std::error::Error for BernoulliError {}
95
96impl Bernoulli {
97 /// Construct a new `Bernoulli` with the given probability of success `p`.
98 ///
99 /// # Precision
100 ///
101 /// For `p = 1.0`, the resulting distribution will always generate true.
102 /// For `p = 0.0`, the resulting distribution will always generate false.
103 ///
104 /// This method is accurate for any input `p` in the range `[0, 1]` which is
105 /// a multiple of 2<sup>-64</sup>. (Note that not all multiples of
106 /// 2<sup>-64</sup> in `[0, 1]` can be represented as a `f64`.)
107 #[inline]
108 pub fn new(p: f64) -> Result<Bernoulli, BernoulliError> {
109 if !(0.0..1.0).contains(&p) {
110 if p == 1.0 {
111 return Ok(Bernoulli { p_int: ALWAYS_TRUE });
112 }
113 return Err(BernoulliError::InvalidProbability);
114 }
115 Ok(Bernoulli {
116 p_int: (p * SCALE) as u64,
117 })
118 }
119
120 /// Construct a new `Bernoulli` with the probability of success of
121 /// `numerator`-in-`denominator`. I.e. `new_ratio(2, 3)` will return
122 /// a `Bernoulli` with a 2-in-3 chance, or about 67%, of returning `true`.
123 ///
124 /// return `true`. If `numerator == 0` it will always return `false`.
125 /// For `numerator > denominator` and `denominator == 0`, this returns an
126 /// error. Otherwise, for `numerator == denominator`, samples are always
127 /// true; for `numerator == 0` samples are always false.
128 #[inline]
129 pub fn from_ratio(numerator: u32, denominator: u32) -> Result<Bernoulli, BernoulliError> {
130 if numerator > denominator || denominator == 0 {
131 return Err(BernoulliError::InvalidProbability);
132 }
133 if numerator == denominator {
134 return Ok(Bernoulli { p_int: ALWAYS_TRUE });
135 }
136 let p_int = ((f64::from(numerator) / f64::from(denominator)) * SCALE) as u64;
137 Ok(Bernoulli { p_int })
138 }
139
140 #[inline]
141 /// Returns the probability (`p`) of the distribution.
142 ///
143 /// This value may differ slightly from the input due to loss of precision.
144 pub fn p(&self) -> f64 {
145 if self.p_int == ALWAYS_TRUE {
146 1.0
147 } else {
148 (self.p_int as f64) / SCALE
149 }
150 }
151}
152
153impl Distribution<bool> for Bernoulli {
154 #[inline]
155 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> bool {
156 // Make sure to always return true for p = 1.0.
157 if self.p_int == ALWAYS_TRUE {
158 return true;
159 }
160 let v: u64 = rng.random();
161 v < self.p_int
162 }
163}
164
165#[cfg(test)]
166mod test {
167 use super::Bernoulli;
168 use crate::distr::Distribution;
169 use crate::Rng;
170
171 #[test]
172 #[cfg(feature = "serde")]
173 fn test_serializing_deserializing_bernoulli() {
174 let coin_flip = Bernoulli::new(0.5).unwrap();
175 let de_coin_flip: Bernoulli =
176 bincode::deserialize(&bincode::serialize(&coin_flip).unwrap()).unwrap();
177
178 assert_eq!(coin_flip.p_int, de_coin_flip.p_int);
179 }
180
181 #[test]
182 fn test_trivial() {
183 // We prefer to be explicit here.
184 #![allow(clippy::bool_assert_comparison)]
185
186 let mut r = crate::test::rng(1);
187 let always_false = Bernoulli::new(0.0).unwrap();
188 let always_true = Bernoulli::new(1.0).unwrap();
189 for _ in 0..5 {
190 assert_eq!(r.sample::<bool, _>(&always_false), false);
191 assert_eq!(r.sample::<bool, _>(&always_true), true);
192 assert_eq!(Distribution::<bool>::sample(&always_false, &mut r), false);
193 assert_eq!(Distribution::<bool>::sample(&always_true, &mut r), true);
194 }
195 }
196
197 #[test]
198 #[cfg_attr(miri, ignore)] // Miri is too slow
199 fn test_average() {
200 const P: f64 = 0.3;
201 const NUM: u32 = 3;
202 const DENOM: u32 = 10;
203 let d1 = Bernoulli::new(P).unwrap();
204 let d2 = Bernoulli::from_ratio(NUM, DENOM).unwrap();
205 const N: u32 = 100_000;
206
207 let mut sum1: u32 = 0;
208 let mut sum2: u32 = 0;
209 let mut rng = crate::test::rng(2);
210 for _ in 0..N {
211 if d1.sample(&mut rng) {
212 sum1 += 1;
213 }
214 if d2.sample(&mut rng) {
215 sum2 += 1;
216 }
217 }
218 let avg1 = (sum1 as f64) / (N as f64);
219 assert!((avg1 - P).abs() < 5e-3);
220
221 let avg2 = (sum2 as f64) / (N as f64);
222 assert!((avg2 - (NUM as f64) / (DENOM as f64)).abs() < 5e-3);
223 }
224
225 #[test]
226 fn value_stability() {
227 let mut rng = crate::test::rng(3);
228 let distr = Bernoulli::new(0.4532).unwrap();
229 let mut buf = [false; 10];
230 for x in &mut buf {
231 *x = rng.sample(distr);
232 }
233 assert_eq!(
234 buf,
235 [true, false, false, true, false, false, true, true, true, true]
236 );
237 }
238
239 #[test]
240 fn bernoulli_distributions_can_be_compared() {
241 assert_eq!(Bernoulli::new(1.0), Bernoulli::new(1.0));
242 }
243}
244