statrs/distribution/
triangular.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::statistics::*;
3use crate::{Result, StatsError};
4use rand::Rng;
5use std::f64;
6
7/// Implements the
8/// [Triangular](https://en.wikipedia.org/wiki/Triangular_distribution)
9/// distribution
10///
11/// # Examples
12///
13/// ```
14/// use statrs::distribution::{Triangular, Continuous};
15/// use statrs::statistics::Distribution;
16///
17/// let n = Triangular::new(0.0, 5.0, 2.5).unwrap();
18/// assert_eq!(n.mean().unwrap(), 7.5 / 3.0);
19/// assert_eq!(n.pdf(2.5), 5.0 / 12.5);
20/// ```
21#[derive(Debug, Copy, Clone, PartialEq)]
22pub struct Triangular {
23    min: f64,
24    max: f64,
25    mode: f64,
26}
27
28impl Triangular {
29    /// Constructs a new triangular distribution with a minimum of `min`,
30    /// maximum of `max`, and a mode of `mode`.
31    ///
32    /// # Errors
33    ///
34    /// Returns an error if `min`, `max`, or `mode` are `NaN` or `±INF`.
35    /// Returns an error if `max < mode`, `mode < min`, or `max == min`.
36    ///
37    /// # Examples
38    ///
39    /// ```
40    /// use statrs::distribution::Triangular;
41    ///
42    /// let mut result = Triangular::new(0.0, 5.0, 2.5);
43    /// assert!(result.is_ok());
44    ///
45    /// result = Triangular::new(2.5, 1.5, 0.0);
46    /// assert!(result.is_err());
47    /// ```
48    pub fn new(min: f64, max: f64, mode: f64) -> Result<Triangular> {
49        if !min.is_finite() || !max.is_finite() || !mode.is_finite() {
50            return Err(StatsError::BadParams);
51        }
52        if max < mode || mode < min {
53            return Err(StatsError::BadParams);
54        }
55        if ulps_eq!(max, min, max_ulps = 0) {
56            return Err(StatsError::BadParams);
57        }
58        Ok(Triangular { min, max, mode })
59    }
60}
61
62impl ::rand::distributions::Distribution<f64> for Triangular {
63    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
64        sample_unchecked(rng, self.min, self.max, self.mode)
65    }
66}
67
68impl ContinuousCDF<f64, f64> for Triangular {
69    /// Calculates the cumulative distribution function for the triangular
70    /// distribution
71    /// at `x`
72    ///
73    /// # Formula
74    ///
75    /// ```ignore
76    /// if x == min {
77    ///     0
78    /// } if min < x <= mode {
79    ///     (x - min)^2 / ((max - min) * (mode - min))
80    /// } else if mode < x < max {
81    ///     1 - (max - min)^2 / ((max - min) * (max - mode))
82    /// } else {
83    ///     1
84    /// }
85    /// ```
86    fn cdf(&self, x: f64) -> f64 {
87        let a = self.min;
88        let b = self.max;
89        let c = self.mode;
90        if x <= a {
91            0.0
92        } else if x <= c {
93            (x - a) * (x - a) / ((b - a) * (c - a))
94        } else if x < b {
95            1.0 - (b - x) * (b - x) / ((b - a) * (b - c))
96        } else {
97            1.0
98        }
99    }
100
101    /// Calculates the survival function for the triangular
102    /// distribution at `x`
103    ///
104    /// # Formula
105    ///
106    /// ```ignore
107    /// if x == min {
108    ///     1
109    /// } if min < x <= mode {
110    ///     1 - (x - min)^2 / ((max - min) * (mode - min))
111    /// } else if mode < x < max {
112    ///     (max - min)^2 / ((max - min) * (max - mode))
113    /// } else {
114    ///     0
115    /// }
116    /// ```
117    fn sf(&self, x: f64) -> f64 {
118        let a = self.min;
119        let b = self.max;
120        let c = self.mode;
121        if x <= a {
122            1.0
123        } else if x <= c {
124            1.0 - ((x - a) * (x - a) / ((b - a) * (c - a)))
125        } else if x < b {
126            (b - x) * (b - x) / ((b - a) * (b - c))
127        } else {
128            0.0
129        }
130    }
131}
132
133impl Min<f64> for Triangular {
134    /// Returns the minimum value in the domain of the
135    /// triangular distribution representable by a double precision float
136    ///
137    /// # Remarks
138    ///
139    /// The return value is the same min used to construct the distribution
140    fn min(&self) -> f64 {
141        self.min
142    }
143}
144
145impl Max<f64> for Triangular {
146    /// Returns the maximum value in the domain of the
147    /// triangular distribution representable by a double precision float
148    ///
149    /// # Remarks
150    ///
151    /// The return value is the same max used to construct the distribution
152    fn max(&self) -> f64 {
153        self.max
154    }
155}
156
157impl Distribution<f64> for Triangular {
158    /// Returns the mean of the triangular distribution
159    ///
160    /// # Formula
161    ///
162    /// ```ignore
163    /// (min + max + mode) / 3
164    /// ```
165    fn mean(&self) -> Option<f64> {
166        Some((self.min + self.max + self.mode) / 3.0)
167    }
168    /// Returns the variance of the triangular distribution
169    ///
170    /// # Formula
171    ///
172    /// ```ignore
173    /// (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode) / 18
174    /// ```
175    fn variance(&self) -> Option<f64> {
176        let a = self.min;
177        let b = self.max;
178        let c = self.mode;
179        Some((a * a + b * b + c * c - a * b - a * c - b * c) / 18.0)
180    }
181    /// Returns the entropy of the triangular distribution
182    ///
183    /// # Formula
184    ///
185    /// ```ignore
186    /// 1 / 2 + ln((max - min) / 2)
187    /// ```
188    fn entropy(&self) -> Option<f64> {
189        Some(0.5 + ((self.max - self.min) / 2.0).ln())
190    }
191    /// Returns the skewness of the triangular distribution
192    ///
193    /// # Formula
194    ///
195    /// ```ignore
196    /// (sqrt(2) * (min + max - 2 * mode) * (2 * min - max - mode) * (min - 2 *
197    /// max + mode)) /
198    /// ( 5 * (min^2 + max^2 + mode^2 - min * max - min * mode - max * mode)^(3
199    /// / 2))
200    /// ```
201    fn skewness(&self) -> Option<f64> {
202        let a = self.min;
203        let b = self.max;
204        let c = self.mode;
205        let q = f64::consts::SQRT_2 * (a + b - 2.0 * c) * (2.0 * a - b - c) * (a - 2.0 * b + c);
206        let d = 5.0 * (a * a + b * b + c * c - a * b - a * c - b * c).powf(3.0 / 2.0);
207        Some(q / d)
208    }
209}
210
211impl Median<f64> for Triangular {
212    /// Returns the median of the triangular distribution
213    ///
214    /// # Formula
215    ///
216    /// ```ignore
217    /// if mode >= (min + max) / 2 {
218    ///     min + sqrt((max - min) * (mode - min) / 2)
219    /// } else {
220    ///     max - sqrt((max - min) * (max - mode) / 2)
221    /// }
222    /// ```
223    fn median(&self) -> f64 {
224        let a = self.min;
225        let b = self.max;
226        let c = self.mode;
227        if c >= (a + b) / 2.0 {
228            a + ((b - a) * (c - a) / 2.0).sqrt()
229        } else {
230            b - ((b - a) * (b - c) / 2.0).sqrt()
231        }
232    }
233}
234
235impl Mode<Option<f64>> for Triangular {
236    /// Returns the mode of the triangular distribution
237    ///
238    /// # Formula
239    ///
240    /// ```ignore
241    /// mode
242    /// ```
243    fn mode(&self) -> Option<f64> {
244        Some(self.mode)
245    }
246}
247
248impl Continuous<f64, f64> for Triangular {
249    /// Calculates the probability density function for the triangular
250    /// distribution
251    /// at `x`
252    ///
253    /// # Formula
254    ///
255    /// ```ignore
256    /// if x < min {
257    ///     0
258    /// } else if min <= x <= mode {
259    ///     2 * (x - min) / ((max - min) * (mode - min))
260    /// } else if mode < x <= max {
261    ///     2 * (max - x) / ((max - min) * (max - mode))
262    /// } else {
263    ///     0
264    /// }
265    /// ```
266    fn pdf(&self, x: f64) -> f64 {
267        let a = self.min;
268        let b = self.max;
269        let c = self.mode;
270        if a <= x && x <= c {
271            2.0 * (x - a) / ((b - a) * (c - a))
272        } else if c < x && x <= b {
273            2.0 * (b - x) / ((b - a) * (b - c))
274        } else {
275            0.0
276        }
277    }
278
279    /// Calculates the log probability density function for the triangular
280    /// distribution
281    /// at `x`
282    ///
283    /// # Formula
284    ///
285    /// ```ignore
286    /// ln( if x < min {
287    ///     0
288    /// } else if min <= x <= mode {
289    ///     2 * (x - min) / ((max - min) * (mode - min))
290    /// } else if mode < x <= max {
291    ///     2 * (max - x) / ((max - min) * (max - mode))
292    /// } else {
293    ///     0
294    /// } )
295    /// ```
296    fn ln_pdf(&self, x: f64) -> f64 {
297        self.pdf(x).ln()
298    }
299}
300
301fn sample_unchecked<R: Rng + ?Sized>(rng: &mut R, min: f64, max: f64, mode: f64) -> f64 {
302    let f: f64 = rng.gen();
303    if f < (mode - min) / (max - min) {
304        min + (f * (max - min) * (mode - min)).sqrt()
305    } else {
306        max - ((1.0 - f) * (max - min) * (max - mode)).sqrt()
307    }
308}
309
310#[rustfmt::skip]
311#[cfg(all(test, feature = "nightly"))]
312mod tests {
313    use std::fmt::Debug;
314    use crate::statistics::*;
315    use crate::distribution::{ContinuousCDF, Continuous, Triangular};
316    use crate::distribution::internal::*;
317    use crate::consts::ACC;
318
319    fn try_create(min: f64, max: f64, mode: f64) -> Triangular {
320        let n = Triangular::new(min, max, mode);
321        assert!(n.is_ok());
322        n.unwrap()
323    }
324
325    fn create_case(min: f64, max: f64, mode: f64) {
326        let n = try_create(min, max, mode);
327        assert_eq!(n.min(), min);
328        assert_eq!(n.max(), max);
329        assert_eq!(n.mode().unwrap(), mode);
330    }
331
332    fn bad_create_case(min: f64, max: f64, mode: f64) {
333        let n = Triangular::new(min, max, mode);
334        assert!(n.is_err());
335    }
336
337    fn get_value<T, F>(min: f64, max: f64, mode: f64, eval: F) -> T
338        where T: PartialEq + Debug,
339              F: Fn(Triangular) -> T
340    {
341        let n = try_create(min, max, mode);
342        eval(n)
343    }
344
345    fn test_case<F>(min: f64, max: f64, mode: f64, expected: f64, eval: F)
346        where F: Fn(Triangular) -> f64
347    {
348        let x = get_value(min, max, mode, eval);
349        assert_eq!(expected, x);
350    }
351
352    fn test_almost<F>(min: f64, max: f64, mode: f64, expected: f64, acc: f64, eval: F)
353        where F: Fn(Triangular) -> f64
354    {
355        let x = get_value(min, max, mode, eval);
356        assert_almost_eq!(expected, x, acc);
357    }
358
359    #[test]
360    fn test_create() {
361        create_case(-1.0, 1.0, 0.0);
362        create_case(1.0, 2.0, 1.0);
363        create_case(5.0, 25.0, 25.0);
364        create_case(1.0e-5, 1.0e5, 1.0e-3);
365        create_case(0.0, 1.0, 0.9);
366        create_case(-4.0, -0.5, -2.0);
367        create_case(-13.039, 8.42, 1.17);
368    }
369
370    #[test]
371    fn test_bad_create() {
372        bad_create_case(0.0, 0.0, 0.0);
373        bad_create_case(0.0, 1.0, -0.1);
374        bad_create_case(0.0, 1.0, 1.1);
375        bad_create_case(0.0, -1.0, 0.5);
376        bad_create_case(2.0, 1.0, 1.5);
377        bad_create_case(f64::NAN, 1.0, 0.5);
378        bad_create_case(0.2, f64::NAN, 0.5);
379        bad_create_case(0.5, 1.0, f64::NAN);
380        bad_create_case(f64::NAN, f64::NAN, f64::NAN);
381        bad_create_case(f64::NEG_INFINITY, 1.0, 0.5);
382        bad_create_case(0.0, f64::INFINITY, 0.5);
383    }
384
385    #[test]
386    fn test_variance() {
387        let variance = |x: Triangular| x.variance().unwrap();
388        test_case(0.0, 1.0, 0.5, 0.75 / 18.0, variance);
389        test_case(0.0, 1.0, 0.75, 0.8125 / 18.0, variance);
390        test_case(-5.0, 8.0, -3.5, 151.75 / 18.0, variance);
391        test_case(-5.0, 8.0, 5.0, 139.0 / 18.0, variance);
392        test_case(-5.0, -3.0, -4.0, 3.0 / 18.0, variance);
393        test_case(15.0, 134.0, 21.0, 13483.0 / 18.0, variance);
394    }
395
396    #[test]
397    fn test_entropy() {
398        let entropy = |x: Triangular| x.entropy().unwrap();
399        test_almost(0.0, 1.0, 0.5, -0.1931471805599453094172, 1e-16, entropy);
400        test_almost(0.0, 1.0, 0.75, -0.1931471805599453094172, 1e-16, entropy);
401        test_case(-5.0, 8.0, -3.5, 2.371802176901591426636, entropy);
402        test_case(-5.0, 8.0, 5.0, 2.371802176901591426636, entropy);
403        test_case(-5.0, -3.0, -4.0, 0.5, entropy);
404        test_case(15.0, 134.0, 21.0, 4.585976312551584075938, entropy);
405    }
406
407    #[test]
408    fn test_skewness() {
409        let skewness = |x: Triangular| x.skewness().unwrap();
410        test_case(0.0, 1.0, 0.5, 0.0, skewness);
411        test_case(0.0, 1.0, 0.75, -0.4224039833745502226059, skewness);
412        test_case(-5.0, 8.0, -3.5, 0.5375093589712976359809, skewness);
413        test_case(-5.0, 8.0, 5.0, -0.4445991743012595633537, skewness);
414        test_case(-5.0, -3.0, -4.0, 0.0, skewness);
415        test_case(15.0, 134.0, 21.0, 0.5605920922751860613217, skewness);
416    }
417
418    #[test]
419    fn test_mode() {
420        let mode = |x: Triangular| x.mode().unwrap();
421        test_case(0.0, 1.0, 0.5, 0.5, mode);
422        test_case(0.0, 1.0, 0.75, 0.75, mode);
423        test_case(-5.0, 8.0, -3.5, -3.5, mode);
424        test_case(-5.0, 8.0, 5.0, 5.0, mode);
425        test_case(-5.0, -3.0, -4.0, -4.0, mode);
426        test_case(15.0, 134.0, 21.0, 21.0, mode);
427    }
428
429    #[test]
430    fn test_median() {
431        let median = |x: Triangular| x.median();
432        test_case(0.0, 1.0, 0.5, 0.5, median);
433        test_case(0.0, 1.0, 0.75, 0.6123724356957945245493, median);
434        test_almost(-5.0, 8.0, -3.5, -0.6458082328952913226724, 1e-15, median);
435        test_almost(-5.0, 8.0, 5.0, 3.062257748298549652367, 1e-15, median);
436        test_case(-5.0, -3.0, -4.0, -4.0, median);
437        test_almost(15.0, 134.0, 21.0, 52.00304883716712238797, 1e-14, median);
438    }
439
440    #[test]
441    fn test_pdf() {
442        let pdf = |arg: f64| move |x: Triangular| x.pdf(arg);
443        test_case(0.0, 1.0, 0.5, 0.0, pdf(-1.0));
444        test_case(0.0, 1.0, 0.5, 0.0, pdf(1.1));
445        test_case(0.0, 1.0, 0.5, 1.0, pdf(0.25));
446        test_case(0.0, 1.0, 0.5, 2.0, pdf(0.5));
447        test_case(0.0, 1.0, 0.5, 1.0, pdf(0.75));
448        test_case(-5.0, 8.0, -3.5, 0.0, pdf(-5.1));
449        test_case(-5.0, 8.0, -3.5, 0.0, pdf(8.1));
450        test_case(-5.0, 8.0, -3.5, 0.1025641025641025641026, pdf(-4.0));
451        test_case(-5.0, 8.0, -3.5, 0.1538461538461538461538, pdf(-3.5));
452        test_case(-5.0, 8.0, -3.5, 0.05351170568561872909699, pdf(4.0));
453        test_case(-5.0, -3.0, -4.0, 0.0, pdf(-5.1));
454        test_case(-5.0, -3.0, -4.0, 0.0, pdf(-2.9));
455        test_case(-5.0, -3.0, -4.0, 0.5, pdf(-4.5));
456        test_case(-5.0, -3.0, -4.0, 1.0, pdf(-4.0));
457        test_case(-5.0, -3.0, -4.0, 0.5, pdf(-3.5));
458    }
459
460    #[test]
461    fn test_ln_pdf() {
462        let ln_pdf = |arg: f64| move |x: Triangular| x.ln_pdf(arg);
463        test_case(0.0, 1.0, 0.5, f64::NEG_INFINITY, ln_pdf(-1.0));
464        test_case(0.0, 1.0, 0.5, f64::NEG_INFINITY, ln_pdf(1.1));
465        test_case(0.0, 1.0, 0.5, 0.0, ln_pdf(0.25));
466        test_case(0.0, 1.0, 0.5, 2f64.ln(), ln_pdf(0.5));
467        test_case(0.0, 1.0, 0.5, 0.0, ln_pdf(0.75));
468        test_case(-5.0, 8.0, -3.5, f64::NEG_INFINITY, ln_pdf(-5.1));
469        test_case(-5.0, 8.0, -3.5, f64::NEG_INFINITY, ln_pdf(8.1));
470        test_case(-5.0, 8.0, -3.5, 0.1025641025641025641026f64.ln(), ln_pdf(-4.0));
471        test_case(-5.0, 8.0, -3.5, 0.1538461538461538461538f64.ln(), ln_pdf(-3.5));
472        test_case(-5.0, 8.0, -3.5, 0.05351170568561872909699f64.ln(), ln_pdf(4.0));
473        test_case(-5.0, -3.0, -4.0, f64::NEG_INFINITY, ln_pdf(-5.1));
474        test_case(-5.0, -3.0, -4.0, f64::NEG_INFINITY, ln_pdf(-2.9));
475        test_case(-5.0, -3.0, -4.0, 0.5f64.ln(), ln_pdf(-4.5));
476        test_case(-5.0, -3.0, -4.0, 0.0, ln_pdf(-4.0));
477        test_case(-5.0, -3.0, -4.0, 0.5f64.ln(), ln_pdf(-3.5));
478    }
479
480    #[test]
481    fn test_cdf() {
482        let cdf = |arg: f64| move |x: Triangular| x.cdf(arg);
483        test_case(0.0, 1.0, 0.5, 0.125, cdf(0.25));
484        test_case(0.0, 1.0, 0.5, 0.5, cdf(0.5));
485        test_case(0.0, 1.0, 0.5, 0.875, cdf(0.75));
486        test_case(-5.0, 8.0, -3.5, 0.05128205128205128205128, cdf(-4.0));
487        test_case(-5.0, 8.0, -3.5, 0.1153846153846153846154, cdf(-3.5));
488        test_case(-5.0, 8.0, -3.5, 0.892976588628762541806, cdf(4.0));
489        test_case(-5.0, -3.0, -4.0, 0.125, cdf(-4.5));
490        test_case(-5.0, -3.0, -4.0, 0.5, cdf(-4.0));
491        test_case(-5.0, -3.0, -4.0, 0.875, cdf(-3.5));
492    }
493
494    #[test]
495    fn test_cdf_lower_bound() {
496        let cdf = |arg: f64| move |x: Triangular| x.cdf(arg);
497        test_case(0.0, 3.0, 1.5, 0.0, cdf(-1.0));
498    }
499
500    #[test]
501    fn test_cdf_upper_bound() {
502        let cdf = |arg: f64| move |x: Triangular| x.cdf(arg);
503        test_case(0.0, 3.0, 1.5, 1.0, cdf(5.0));
504    }
505
506
507    #[test]
508    fn test_sf() {
509        let sf = |arg: f64| move |x: Triangular| x.sf(arg);
510        test_case(0.0, 1.0, 0.5, 0.875, sf(0.25));
511        test_case(0.0, 1.0, 0.5, 0.5, sf(0.5));
512        test_case(0.0, 1.0, 0.5, 0.125, sf(0.75));
513        test_case(-5.0, 8.0, -3.5, 0.9487179487179487, sf(-4.0));
514        test_case(-5.0, 8.0, -3.5, 0.8846153846153846, sf(-3.5));
515        test_case(-5.0, 8.0, -3.5, 0.10702341137123746, sf(4.0));
516        test_case(-5.0, -3.0, -4.0, 0.875, sf(-4.5));
517        test_case(-5.0, -3.0, -4.0, 0.5, sf(-4.0));
518        test_case(-5.0, -3.0, -4.0, 0.125, sf(-3.5));
519    }
520
521    #[test]
522    fn test_sf_lower_bound() {
523        let sf = |arg: f64| move |x: Triangular| x.sf(arg);
524        test_case(0.0, 3.0, 1.5, 1.0, sf(-1.0));
525    }
526
527    #[test]
528    fn test_sf_upper_bound() {
529        let sf = |arg: f64| move |x: Triangular| x.sf(arg);
530        test_case(0.0, 3.0, 1.5, 0.0, sf(5.0));
531    }
532
533    #[test]
534    fn test_continuous() {
535        test::check_continuous_distribution(&try_create(-5.0, 5.0, 0.0), -5.0, 5.0);
536        test::check_continuous_distribution(&try_create(-15.0, -2.0, -3.0), -15.0, -2.0);
537    }
538}
statrs/distribution/triangular.rs

statrs/distribution/
triangular.rs
