r2rs_nmath/distribution/norm/
mod.rs

1mod d;
2mod p;
3mod q;
4mod r;
5
6use strafe_type::{LogProbability64, Positive64, Probability64, Real64};
7
8pub(crate) use self::{d::*, p::*, q::*, r::*};
9use crate::traits::{Distribution, RNG};
10
11/// # The Normal Distribution
12///
13/// ## Description
14///
15/// Density, distribution function, quantile function and random generation for the normal
16/// distribution with mean equal to mean and standard deviation equal to sd.
17///
18/// ## Arguments
19///
20/// * mean: vector of means.
21/// * sd: vector of standard deviations.
22///
23/// ## Details
24///
25/// If mean or sd are not specified they assume the default values of 0 and 1, respectively.
26///
27/// The normal distribution has density
28///
29/// $ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(x - \mu)^2}{2 \sigma^2}} $
30///
31/// where $ \mu $ is the mean of the distribution and $ \sigma $ the standard deviation.
32///
33/// ## Value
34///
35/// dnorm gives the density, pnorm gives the distribution function, qnorm gives the quantile
36/// function, and rnorm generates random deviates.
37///
38/// The length of the result is determined by n for rnorm, and is the maximum of the lengths of
39/// the numerical arguments for the other functions.
40///
41/// The numerical arguments other than n are recycled to the length of the result. Only the first
42/// elements of the logical arguments are used.
43///
44/// For sd = 0 this gives the limit as sd decreases to 0, a point mass at mu. sd < 0 is an error
45/// and returns NaN.
46///
47/// ## Density Plot
48///
49/// ```rust
50/// # use r2rs_base::traits::StatisticalSlice;
51/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
52/// # use strafe_plot::prelude::{IntoDrawingArea, Line, Plot, PlotOptions, SVGBackend, BLACK};
53/// # use strafe_type::FloatConstraint;
54/// let norm = NormalBuilder::new().build();
55/// let x = <[f64]>::sequence(-3.0, 3.0, 1000);
56/// let y = x
57///     .iter()
58///     .map(|x| norm.density(x).unwrap())
59///     .collect::<Vec<_>>();
60///
61/// let root = SVGBackend::new("density.svg", (1024, 768)).into_drawing_area();
62/// Plot::new()
63///     .with_options(PlotOptions {
64///         x_axis_label: "x".to_string(),
65///         y_axis_label: "density".to_string(),
66///         ..Default::default()
67///     })
68///     .with_plottable(Line {
69///         x,
70///         y,
71///         color: BLACK,
72///         ..Default::default()
73///     })
74///     .plot(&root)
75///     .unwrap();
76/// # use std::fs::rename;
77/// #     drop(root);
78/// #     rename(
79/// #             format!("density.svg"),
80/// #             format!("src/distribution/norm/doctest_out/density.svg"),
81/// #     )
82/// #     .unwrap();
83/// ```
84#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("density", "src/distribution/norm/doctest_out/density.svg")))]
85#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Density][density]"))]
86///
87/// ## Source
88///
89/// For pnorm, based on
90///
91/// Cody, W. D. (1993) Algorithm 715: SPECFUN – A portable FORTRAN package of special function
92/// routines and test drivers. ACM Transactions on Mathematical Software 19, 22–32.
93///
94/// For qnorm, the code is a C translation of
95///
96/// Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the normal distribution.
97/// Applied Statistics, 37, 477–484.
98///
99/// which provides precise results up to about 16 digits.
100///
101/// For rnorm, see RNG for how to select the algorithm and for references to the supplie
102/// methods.
103///
104/// ## References
105///
106/// Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth &
107/// Brooks/Cole.
108///
109/// Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions,
110/// volume 1, chapter 13. Wiley, New York.
111///
112/// ## See Also
113///
114/// Distributions for other standard distributions, including dlnorm for the Lognormal distribution.
115///
116/// ## Examples
117///
118/// ```rust
119/// # use num_traits::FloatConst;
120/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
121/// let norm = NormalBuilder::new().build();
122/// println!("{}", norm.density(0));
123/// println!("{}", 1.0 / (2.0 * f64::PI()).sqrt());
124/// # use std::{fs::File, io::Write};
125/// # let mut f = File::create("src/distribution/norm/doctest_out/dens1.md").unwrap();
126/// # writeln!(f, "```output").unwrap();
127/// # writeln!(f, "{}", norm.density(0)).unwrap();
128/// # writeln!(f, "{}", 1.0 / (2.0 * f64::PI()).sqrt()).unwrap();
129/// # writeln!(f, "```").unwrap();
130/// ```
131#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/dens1.md")))]
132///
133/// ```rust
134/// # use num_traits::FloatConst;
135/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
136/// let norm = NormalBuilder::new().build();
137/// println!("{}", norm.density(1));
138/// println!("{}", (-0.5_f64).exp() / (2.0 * f64::PI()).sqrt());
139/// println!("{}", 1.0 / (2.0 * f64::PI() * 1.0_f64.exp()).sqrt());
140/// # use std::{fs::File, io::Write};
141/// # let mut f = File::create("src/distribution/norm/doctest_out/dens2.md").unwrap();
142/// # writeln!(f, "```output").unwrap();
143/// # writeln!(f, "{}", norm.density(1)).unwrap();
144/// # writeln!(f, "{}", (-0.5_f64).exp() / (2.0 * f64::PI()).sqrt()).unwrap();
145/// # writeln!(f, "{}", 1.0 / (2.0 * f64::PI() * 1.0_f64.exp()).sqrt()).unwrap();
146/// # writeln!(f, "```").unwrap();
147/// ```
148#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/dens2.md")))]
149///
150/// ```rust
151/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
152/// # use strafe_plot::{
153/// #     plot::Plot,
154/// #     prelude::{IntoDrawingArea, Line, PlotOptions, SVGBackend, BLACK, RED},
155/// # };
156/// # use strafe_type::FloatConstraint;
157/// let norm = NormalBuilder::new().build();
158///
159/// let dx = (-60..50).map(|x| x as f64).collect::<Vec<_>>();
160/// let dy1 = dx
161///     .iter()
162///     .map(|x| norm.log_density(x).unwrap())
163///     .collect::<Vec<_>>();
164/// let dy2 = dx
165///     .iter()
166///     .map(|x| norm.density(x).unwrap().ln())
167///     .collect::<Vec<_>>();
168///
169/// let px = (-50..10).map(|x| x as f64).collect::<Vec<_>>();
170/// let py1 = px
171///     .iter()
172///     .map(|x| norm.log_probability(x, true).unwrap())
173///     .collect::<Vec<_>>();
174/// let py2 = px
175///     .iter()
176///     .map(|x| norm.probability(x, true).unwrap().ln())
177///     .collect::<Vec<_>>();
178///
179/// let root = SVGBackend::new("log_plots.svg", (1024, 768)).into_drawing_area();
180///
181/// Plot::new()
182///     .with_options(PlotOptions {
183///         x_axis_label: "x".to_string(),
184///         y_axis_label: "density".to_string(),
185///         title: "Log Normal Density".to_string(),
186///         legend_x: 0.1,
187///         legend_y: 0.1,
188///         plot_bottom: 0.5,
189///         ..Default::default()
190///     })
191///     .with_plottable(Line {
192///         x: dx.clone(),
193///         y: dy1,
194///         legend: true,
195///         label: "log_density()".to_string(),
196///         color: BLACK,
197///         ..Default::default()
198///     })
199///     .with_plottable(Line {
200///         x: dx.clone(),
201///         y: dy2,
202///         legend: true,
203///         label: "density().ln()".to_string(),
204///         color: RED,
205///         ..Default::default()
206///     })
207///     .plot(&root)
208///     .unwrap();
209///
210/// Plot::new()
211///     .with_options(PlotOptions {
212///         x_axis_label: "x".to_string(),
213///         y_axis_label: "probability".to_string(),
214///         title: "Log Normal Cumulative".to_string(),
215///         legend_x: 0.1,
216///         legend_y: 0.1,
217///         plot_top: 0.5,
218///         ..Default::default()
219///     })
220///     .with_plottable(Line {
221///         x: px.clone(),
222///         y: py1,
223///         legend: true,
224///         label: "log_probability()".to_string(),
225///         color: BLACK,
226///         ..Default::default()
227///     })
228///     .with_plottable(Line {
229///         x: px.clone(),
230///         y: py2,
231///         legend: true,
232///         label: "probability().ln()".to_string(),
233///         color: RED,
234///         ..Default::default()
235///     })
236///     .plot(&root)
237///     .unwrap();
238/// # use std::fs::rename;
239/// #     drop(root);
240/// #     rename(
241/// #             format!("log_plots.svg"),
242/// #             format!("src/distribution/norm/doctest_out/log_plots.svg"),
243/// #     )
244/// #     .unwrap();
245/// ```
246#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("log_plots", "src/distribution/norm/doctest_out/log_plots.svg")))]
247#[cfg_attr(
248    feature = "doc_outputs",
249    cfg_attr(all(), doc = "![Log Plots][log_plots]")
250)]
251///
252/// If you want the so-called 'error function'
253/// (see Abramowitz and Stegun 29.2.29)
254/// ```rust
255/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
256/// # use strafe_type::FloatConstraint;
257/// fn erf(x: f64) -> f64 {
258///     let norm = NormalBuilder::new().build();
259///     2.0 * norm.probability(x * 2.0_f64.sqrt(), true).unwrap() - 1.0
260/// }
261/// ```
262///
263/// and the so-called 'complementary error function'
264/// ```rust
265/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
266/// # use strafe_type::FloatConstraint;
267/// fn erfc(x: f64) -> f64 {
268///     let norm = NormalBuilder::new().build();
269///     2.0 * norm.probability(x * 2.0_f64.sqrt(), false).unwrap()
270/// }
271/// ```
272///
273/// and the inverses
274/// ```rust
275/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
276/// # use strafe_type::FloatConstraint;
277/// fn erfinv(x: f64) -> f64 {
278///     let norm = NormalBuilder::new().build();
279///     norm.quantile((1.0 + x) / 2.0, true).unwrap() / 2.0_f64.sqrt()
280/// }
281///
282/// fn erfcinv(x: f64) -> f64 {
283///     let norm = NormalBuilder::new().build();
284///     norm.quantile(x / 2.0, false).unwrap() / 2.0_f64.sqrt()
285/// }
286/// ```
287pub struct Normal {
288    mean: Real64,
289    standard_deviation: Positive64,
290}
291
292impl Distribution for Normal {
293    fn density<R: Into<Real64>>(&self, x: R) -> Real64 {
294        dnorm(x, self.mean, self.standard_deviation, false)
295    }
296
297    fn log_density<R: Into<Real64>>(&self, x: R) -> Real64 {
298        dnorm(x, self.mean, self.standard_deviation, true)
299    }
300
301    fn probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> Probability64 {
302        pnorm(q, self.mean, self.standard_deviation, lower_tail)
303    }
304
305    fn log_probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> LogProbability64 {
306        log_pnorm(q, self.mean, self.standard_deviation, lower_tail)
307    }
308
309    fn quantile<P: Into<Probability64>>(&self, p: P, lower_tail: bool) -> Real64 {
310        qnorm(p, self.mean, self.standard_deviation, lower_tail)
311    }
312
313    fn log_quantile<LP: Into<LogProbability64>>(&self, p: LP, lower_tail: bool) -> Real64 {
314        log_qnorm(p, self.mean, self.standard_deviation, lower_tail)
315    }
316
317    fn random_sample<R: RNG>(&self, rng: &mut R) -> Real64 {
318        rnorm(self.mean, self.standard_deviation, rng)
319    }
320}
321
322pub struct NormalBuilder {
323    mean: Option<Real64>,
324    standard_deviation: Option<Positive64>,
325}
326
327impl NormalBuilder {
328    pub fn new() -> Self {
329        Self {
330            mean: None,
331            standard_deviation: None,
332        }
333    }
334
335    pub fn with_mean<R: Into<Real64>>(&mut self, mean: R) -> &mut Self {
336        self.mean = Some(mean.into());
337        self
338    }
339
340    pub fn with_standard_deviation<P: Into<Positive64>>(
341        &mut self,
342        standard_deviation: P,
343    ) -> &mut Self {
344        self.standard_deviation = Some(standard_deviation.into());
345        self
346    }
347
348    pub fn build(&self) -> Normal {
349        let mean = self.mean.unwrap_or(0.0.into());
350        let standard_deviation = self.standard_deviation.unwrap_or(1.0.into());
351
352        Normal {
353            mean,
354            standard_deviation,
355        }
356    }
357}
358
359#[cfg(test)]
360mod tests;
361
362#[cfg(all(test, feature = "enable_proptest"))]
363mod proptests;
364
365#[cfg(all(test, feature = "enable_covtest"))]
366mod covtests;
r2rs_nmath/distribution/norm/mod.rs

r2rs_nmath/distribution/norm/
mod.rs
