ferray_random/distributions/
misc_continuous.rs

1// ferray-random: Miscellaneous continuous distributions
2//
3// laplace, logistic, rayleigh, weibull, pareto, gumbel, power, triangular,
4// vonmises, wald, standard_cauchy
5
6use ferray_core::{Array, FerrayError, IxDyn};
7
8use crate::bitgen::BitGenerator;
9use crate::distributions::exponential::standard_exponential_single;
10use crate::distributions::normal::standard_normal_single;
11use crate::generator::{Generator, generate_vec, shape_size, vec_to_array_f64};
12use crate::shape::IntoShape;
13
14impl<B: BitGenerator> Generator<B> {
15    /// Generate an array of Laplace-distributed variates.
16    ///
17    /// PDF: f(x) = (1/(2*scale)) * exp(-|x - loc| / scale).
18    ///
19    /// # Arguments
20    /// * `loc` - Location parameter.
21    /// * `scale` - Scale parameter, must be positive.
22    /// * `size` - Number of values to generate.
23    ///
24    /// # Errors
25    /// Returns `FerrayError::InvalidValue` if `scale <= 0` or `size` is zero.
26    pub fn laplace(
27        &mut self,
28        loc: f64,
29        scale: f64,
30        size: impl IntoShape,
31    ) -> Result<Array<f64, IxDyn>, FerrayError> {
32        if scale <= 0.0 {
33            return Err(FerrayError::invalid_value(format!(
34                "scale must be positive, got {scale}"
35            )));
36        }
37        let shape_vec = size.into_shape()?;
38        let n = shape_size(&shape_vec);
39        let data = generate_vec(self, n, |bg| {
40            let mut u = bg.next_f64() - 0.5;
41            // Reject u.abs() >= 0.5 to avoid ln(0) = -inf
42            while u.abs() >= 0.5 {
43                u = bg.next_f64() - 0.5;
44            }
45            (scale * u.signum()).mul_add(-2.0f64.mul_add(-u.abs(), 1.0).ln(), loc)
46        });
47        vec_to_array_f64(data, &shape_vec)
48    }
49
50    /// Generate an array of logistic-distributed variates.
51    ///
52    /// Uses inverse CDF: loc + scale * ln(u / (1 - u)).
53    ///
54    /// # Arguments
55    /// * `loc` - Location parameter.
56    /// * `scale` - Scale parameter, must be positive.
57    /// * `size` - Number of values to generate.
58    ///
59    /// # Errors
60    /// Returns `FerrayError::InvalidValue` if `scale <= 0` or `size` is zero.
61    pub fn logistic(
62        &mut self,
63        loc: f64,
64        scale: f64,
65        size: impl IntoShape,
66    ) -> Result<Array<f64, IxDyn>, FerrayError> {
67        if scale <= 0.0 {
68            return Err(FerrayError::invalid_value(format!(
69                "scale must be positive, got {scale}"
70            )));
71        }
72        let shape_vec = size.into_shape()?;
73        let n = shape_size(&shape_vec);
74        let data = generate_vec(self, n, |bg| {
75            loop {
76                let u = bg.next_f64();
77                if u > f64::EPSILON && u < 1.0 - f64::EPSILON {
78                    return scale.mul_add((u / (1.0 - u)).ln(), loc);
79                }
80            }
81        });
82        vec_to_array_f64(data, &shape_vec)
83    }
84
85    /// Generate an array of Rayleigh-distributed variates.
86    ///
87    /// Uses inverse CDF: scale * sqrt(-2 * ln(1-u)).
88    ///
89    /// # Arguments
90    /// * `scale` - Scale parameter, must be positive.
91    /// * `size` - Number of values to generate.
92    ///
93    /// # Errors
94    /// Returns `FerrayError::InvalidValue` if `scale <= 0` or `size` is zero.
95    pub fn rayleigh(
96        &mut self,
97        scale: f64,
98        size: impl IntoShape,
99    ) -> Result<Array<f64, IxDyn>, FerrayError> {
100        if scale <= 0.0 {
101            return Err(FerrayError::invalid_value(format!(
102                "scale must be positive, got {scale}"
103            )));
104        }
105        let shape_vec = size.into_shape()?;
106        let n = shape_size(&shape_vec);
107        let data = generate_vec(self, n, |bg| {
108            scale * (2.0 * standard_exponential_single(bg)).sqrt()
109        });
110        vec_to_array_f64(data, &shape_vec)
111    }
112
113    /// Generate an array of Weibull-distributed variates.
114    ///
115    /// Uses inverse CDF: (-ln(1-u))^(1/a).
116    ///
117    /// # Arguments
118    /// * `a` - Shape parameter, must be positive.
119    /// * `size` - Number of values to generate.
120    ///
121    /// # Errors
122    /// Returns `FerrayError::InvalidValue` if `a <= 0` or `size` is zero.
123    pub fn weibull(
124        &mut self,
125        a: f64,
126        size: impl IntoShape,
127    ) -> Result<Array<f64, IxDyn>, FerrayError> {
128        if a <= 0.0 {
129            return Err(FerrayError::invalid_value(format!(
130                "a must be positive, got {a}"
131            )));
132        }
133        let shape_vec = size.into_shape()?;
134        let n = shape_size(&shape_vec);
135        let data = generate_vec(self, n, |bg| standard_exponential_single(bg).powf(1.0 / a));
136        vec_to_array_f64(data, &shape_vec)
137    }
138
139    /// Generate an array of Pareto (type II / Lomax) distributed variates.
140    ///
141    /// Uses inverse CDF: (1-u)^(-1/a) - 1 (then shifted by 1 to match `NumPy`).
142    /// `NumPy`'s Pareto: samples from Pareto(a) with `x_m=1`, so PDF = a / x^(a+1) for x >= 1.
143    ///
144    /// # Arguments
145    /// * `a` - Shape parameter, must be positive.
146    /// * `size` - Number of values to generate.
147    ///
148    /// # Errors
149    /// Returns `FerrayError::InvalidValue` if `a <= 0` or `size` is zero.
150    pub fn pareto(
151        &mut self,
152        a: f64,
153        size: impl IntoShape,
154    ) -> Result<Array<f64, IxDyn>, FerrayError> {
155        if a <= 0.0 {
156            return Err(FerrayError::invalid_value(format!(
157                "a must be positive, got {a}"
158            )));
159        }
160        let shape_vec = size.into_shape()?;
161        let n = shape_size(&shape_vec);
162        let data = generate_vec(self, n, |bg| {
163            let e = standard_exponential_single(bg);
164            (e / a).exp_m1()
165        });
166        vec_to_array_f64(data, &shape_vec)
167    }
168
169    /// Generate an array of Gumbel-distributed variates.
170    ///
171    /// Uses inverse CDF: loc - scale * ln(-ln(u)).
172    ///
173    /// # Arguments
174    /// * `loc` - Location parameter.
175    /// * `scale` - Scale parameter, must be positive.
176    /// * `size` - Number of values to generate.
177    ///
178    /// # Errors
179    /// Returns `FerrayError::InvalidValue` if `scale <= 0` or `size` is zero.
180    pub fn gumbel(
181        &mut self,
182        loc: f64,
183        scale: f64,
184        size: impl IntoShape,
185    ) -> Result<Array<f64, IxDyn>, FerrayError> {
186        if scale <= 0.0 {
187            return Err(FerrayError::invalid_value(format!(
188                "scale must be positive, got {scale}"
189            )));
190        }
191        let shape_vec = size.into_shape()?;
192        let n = shape_size(&shape_vec);
193        let data = generate_vec(self, n, |bg| {
194            loop {
195                let u = bg.next_f64();
196                if u > f64::EPSILON && u < 1.0 - f64::EPSILON {
197                    return scale.mul_add(-(-u.ln()).ln(), loc);
198                }
199            }
200        });
201        vec_to_array_f64(data, &shape_vec)
202    }
203
204    /// Generate an array of power-distributed variates.
205    ///
206    /// Power distribution with shape `a` on [0, 1]:
207    /// PDF: a * x^(a-1), CDF: x^a. Inverse CDF: u^(1/a).
208    ///
209    /// # Arguments
210    /// * `a` - Shape parameter, must be positive.
211    /// * `size` - Number of values to generate.
212    ///
213    /// # Errors
214    /// Returns `FerrayError::InvalidValue` if `a <= 0` or `size` is zero.
215    pub fn power(
216        &mut self,
217        a: f64,
218        size: impl IntoShape,
219    ) -> Result<Array<f64, IxDyn>, FerrayError> {
220        if a <= 0.0 {
221            return Err(FerrayError::invalid_value(format!(
222                "a must be positive, got {a}"
223            )));
224        }
225        let shape_vec = size.into_shape()?;
226        let n = shape_size(&shape_vec);
227        let data = generate_vec(self, n, |bg| {
228            let e = standard_exponential_single(bg);
229            // 1 - exp(-e) gives a U(0,1), then raise to 1/a
230            // More precisely: u^(1/a) where u ~ U(0,1)
231            // Using exponential: (1 - exp(-e))^(1/a)
232            (1.0 - (-e).exp()).powf(1.0 / a)
233        });
234        vec_to_array_f64(data, &shape_vec)
235    }
236
237    /// Generate an array of triangular-distributed variates.
238    ///
239    /// # Arguments
240    /// * `left` - Lower limit.
241    /// * `mode` - Mode (peak), must be in `[left, right]`.
242    /// * `right` - Upper limit, must be > `left`.
243    /// * `size` - Number of values to generate.
244    ///
245    /// # Errors
246    /// Returns `FerrayError::InvalidValue` if parameters are invalid or `size` is zero.
247    pub fn triangular(
248        &mut self,
249        left: f64,
250        mode: f64,
251        right: f64,
252        size: impl IntoShape,
253    ) -> Result<Array<f64, IxDyn>, FerrayError> {
254        if left >= right {
255            return Err(FerrayError::invalid_value(format!(
256                "left ({left}) must be less than right ({right})"
257            )));
258        }
259        if mode < left || mode > right {
260            return Err(FerrayError::invalid_value(format!(
261                "mode ({mode}) must be in [{left}, {right}]"
262            )));
263        }
264        let fc = (mode - left) / (right - left);
265        let shape_vec = size.into_shape()?;
266        let n = shape_size(&shape_vec);
267        let data = generate_vec(self, n, |bg| {
268            let u = bg.next_f64();
269            if u < fc {
270                left + ((right - left) * (mode - left) * u).sqrt()
271            } else {
272                right - ((right - left) * (right - mode) * (1.0 - u)).sqrt()
273            }
274        });
275        vec_to_array_f64(data, &shape_vec)
276    }
277
278    /// Generate an array of von Mises distributed variates.
279    ///
280    /// The von Mises distribution is a continuous distribution on the circle.
281    /// Uses Best & Fisher's algorithm.
282    ///
283    /// # Arguments
284    /// * `mu` - Mean direction (in radians).
285    /// * `kappa` - Concentration parameter, must be non-negative.
286    /// * `size` - Number of values to generate.
287    ///
288    /// # Errors
289    /// Returns `FerrayError::InvalidValue` if `kappa < 0` or `size` is zero.
290    pub fn vonmises(
291        &mut self,
292        mu: f64,
293        kappa: f64,
294        size: impl IntoShape,
295    ) -> Result<Array<f64, IxDyn>, FerrayError> {
296        if kappa < 0.0 {
297            return Err(FerrayError::invalid_value(format!(
298                "kappa must be non-negative, got {kappa}"
299            )));
300        }
301        let shape_vec = size.into_shape()?;
302        let n = shape_size(&shape_vec);
303
304        if kappa < 1e-6 {
305            // For very small kappa, the distribution is nearly uniform on [-pi, pi)
306            let data = generate_vec(self, n, |bg| {
307                mu + bg
308                    .next_f64()
309                    .mul_add(std::f64::consts::TAU, -std::f64::consts::PI)
310            });
311            return vec_to_array_f64(data, &shape_vec);
312        }
313
314        // Best & Fisher algorithm
315        let tau = 1.0 + (4.0 * kappa).mul_add(kappa, 1.0).sqrt();
316        let rho = (tau - (2.0 * tau).sqrt()) / (2.0 * kappa);
317        let r = (1.0 + rho * rho) / (2.0 * rho);
318
319        let data = generate_vec(self, n, |bg| {
320            loop {
321                let u1 = bg.next_f64();
322                let z = std::f64::consts::TAU
323                    .mul_add(u1, -std::f64::consts::PI)
324                    .cos();
325                let f_val = (1.0 + r * z) / (r + z);
326                let c = kappa * (r - f_val);
327                let u2 = bg.next_f64();
328                if u2 < c * (2.0 - c) || u2 <= c * (-c).exp() {
329                    let u3 = bg.next_f64();
330                    let theta = if u3 > 0.5 {
331                        mu + f_val.acos()
332                    } else {
333                        mu - f_val.acos()
334                    };
335                    return theta;
336                }
337            }
338        });
339        vec_to_array_f64(data, &shape_vec)
340    }
341
342    /// Generate an array of Wald (inverse Gaussian) distributed variates.
343    ///
344    /// # Arguments
345    /// * `mean` - Mean of the distribution, must be positive.
346    /// * `scale` - Scale parameter (lambda), must be positive.
347    /// * `size` - Number of values to generate.
348    ///
349    /// # Errors
350    /// Returns `FerrayError::InvalidValue` if `mean <= 0`, `scale <= 0`, or `size` is zero.
351    pub fn wald(
352        &mut self,
353        mean: f64,
354        scale: f64,
355        size: impl IntoShape,
356    ) -> Result<Array<f64, IxDyn>, FerrayError> {
357        if mean <= 0.0 {
358            return Err(FerrayError::invalid_value(format!(
359                "mean must be positive, got {mean}"
360            )));
361        }
362        if scale <= 0.0 {
363            return Err(FerrayError::invalid_value(format!(
364                "scale must be positive, got {scale}"
365            )));
366        }
367        // Michael, Schucany & Haas algorithm
368        let shape_vec = size.into_shape()?;
369        let n = shape_size(&shape_vec);
370        let data = generate_vec(self, n, |bg| {
371            let z = standard_normal_single(bg);
372            let v = z * z;
373            let mu = mean;
374            let lam = scale;
375            let x = (mu / (2.0 * lam)).mul_add(
376                -(4.0 * mu * lam).mul_add(v, mu * mu * v * v).sqrt(),
377                mu + (mu * mu * v) / (2.0 * lam),
378            );
379            let u = bg.next_f64();
380            if u <= mu / (mu + x) { x } else { mu * mu / x }
381        });
382        vec_to_array_f64(data, &shape_vec)
383    }
384
385    /// Generate an array of standard Cauchy distributed variates.
386    ///
387    /// Uses the inverse CDF: tan(pi * (u - 0.5)).
388    ///
389    /// # Arguments
390    /// * `size` - Number of values to generate.
391    ///
392    /// # Errors
393    /// Returns `FerrayError::InvalidValue` if `size` is zero.
394    pub fn standard_cauchy(
395        &mut self,
396        size: impl IntoShape,
397    ) -> Result<Array<f64, IxDyn>, FerrayError> {
398        let shape_vec = size.into_shape()?;
399        let n = shape_size(&shape_vec);
400        let data = generate_vec(self, n, |bg| {
401            loop {
402                let u = bg.next_f64();
403                // Avoid u = 0.5 exactly (tan(0) = 0, which is fine, but avoid edge)
404                if (u - 0.5).abs() > f64::EPSILON {
405                    return (std::f64::consts::PI * (u - 0.5)).tan();
406                }
407            }
408        });
409        vec_to_array_f64(data, &shape_vec)
410    }
411}
412
413#[cfg(test)]
414#[allow(clippy::unreadable_literal)] // Statistical reference constants are quoted to canonical 16-digit precision.
415mod tests {
416    use crate::default_rng_seeded;
417
418    #[test]
419    fn laplace_mean_variance() {
420        let mut rng = default_rng_seeded(42);
421        let n = 100_000;
422        let loc = 2.0;
423        let scale = 3.0;
424        let arr = rng.laplace(loc, scale, n).unwrap();
425        let slice = arr.as_slice().unwrap();
426        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
427        let var: f64 = slice.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n as f64;
428        // Laplace(loc, scale): mean=loc, var=2*scale^2
429        let expected_var = 2.0 * scale * scale;
430        let se = (expected_var / n as f64).sqrt();
431        assert!(
432            (mean - loc).abs() < 3.0 * se,
433            "laplace mean {mean} too far from {loc}"
434        );
435        assert!(
436            (var - expected_var).abs() / expected_var < 0.05,
437            "laplace variance {var} too far from {expected_var}"
438        );
439    }
440
441    #[test]
442    fn logistic_mean() {
443        let mut rng = default_rng_seeded(42);
444        let n = 100_000;
445        let loc = 1.0;
446        let scale = 2.0;
447        let arr = rng.logistic(loc, scale, n).unwrap();
448        let slice = arr.as_slice().unwrap();
449        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
450        let expected_var = (std::f64::consts::PI * scale).powi(2) / 3.0;
451        let se = (expected_var / n as f64).sqrt();
452        assert!(
453            (mean - loc).abs() < 3.0 * se,
454            "logistic mean {mean} too far from {loc}"
455        );
456    }
457
458    #[test]
459    fn rayleigh_positive() {
460        let mut rng = default_rng_seeded(42);
461        let arr = rng.rayleigh(2.0, 10_000).unwrap();
462        for &v in arr.as_slice().unwrap() {
463            assert!(v > 0.0);
464        }
465    }
466
467    #[test]
468    fn rayleigh_mean() {
469        let mut rng = default_rng_seeded(42);
470        let n = 100_000;
471        let scale = 2.0;
472        let arr = rng.rayleigh(scale, n).unwrap();
473        let slice = arr.as_slice().unwrap();
474        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
475        // Rayleigh(sigma): mean = sigma * sqrt(pi/2)
476        let expected_mean = scale * (std::f64::consts::FRAC_PI_2).sqrt();
477        let expected_var = (4.0 - std::f64::consts::PI) / 2.0 * scale * scale;
478        let se = (expected_var / n as f64).sqrt();
479        assert!(
480            (mean - expected_mean).abs() < 3.0 * se,
481            "rayleigh mean {mean} too far from {expected_mean}"
482        );
483    }
484
485    #[test]
486    fn weibull_positive() {
487        let mut rng = default_rng_seeded(42);
488        let arr = rng.weibull(1.5, 10_000).unwrap();
489        for &v in arr.as_slice().unwrap() {
490            assert!(v > 0.0);
491        }
492    }
493
494    #[test]
495    fn pareto_ge_zero() {
496        let mut rng = default_rng_seeded(42);
497        let arr = rng.pareto(3.0, 10_000).unwrap();
498        for &v in arr.as_slice().unwrap() {
499            assert!(v >= 0.0, "pareto value {v} < 0");
500        }
501    }
502
503    #[test]
504    fn gumbel_mean() {
505        let mut rng = default_rng_seeded(42);
506        let n = 100_000;
507        let loc = 1.0;
508        let scale = 2.0;
509        let arr = rng.gumbel(loc, scale, n).unwrap();
510        let slice = arr.as_slice().unwrap();
511        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
512        // Gumbel: mean = loc + scale * gamma_euler (0.5772...)
513        let euler = 0.5772156649015329;
514        let expected_mean = loc + scale * euler;
515        let expected_var = std::f64::consts::PI * std::f64::consts::PI * scale * scale / 6.0;
516        let se = (expected_var / n as f64).sqrt();
517        assert!(
518            (mean - expected_mean).abs() < 3.0 * se,
519            "gumbel mean {mean} too far from {expected_mean}"
520        );
521    }
522
523    #[test]
524    fn power_in_range() {
525        let mut rng = default_rng_seeded(42);
526        let arr = rng.power(2.0, 10_000).unwrap();
527        for &v in arr.as_slice().unwrap() {
528            assert!((0.0..=1.0).contains(&v), "power value {v} out of [0,1]");
529        }
530    }
531
532    #[test]
533    fn triangular_in_range() {
534        let mut rng = default_rng_seeded(42);
535        let arr = rng.triangular(1.0, 3.0, 5.0, 10_000).unwrap();
536        for &v in arr.as_slice().unwrap() {
537            assert!(
538                (1.0..=5.0).contains(&v),
539                "triangular value {v} out of [1,5]"
540            );
541        }
542    }
543
544    #[test]
545    fn triangular_mean() {
546        let mut rng = default_rng_seeded(42);
547        let n = 100_000;
548        let (left, mode, right) = (1.0, 3.0, 5.0);
549        let arr = rng.triangular(left, mode, right, n).unwrap();
550        let slice = arr.as_slice().unwrap();
551        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
552        let expected_mean = (left + mode + right) / 3.0;
553        assert!(
554            (mean - expected_mean).abs() < 0.05,
555            "triangular mean {mean} too far from {expected_mean}"
556        );
557    }
558
559    #[test]
560    fn vonmises_basic() {
561        let mut rng = default_rng_seeded(42);
562        let arr = rng.vonmises(0.0, 1.0, 10_000).unwrap();
563        assert_eq!(arr.shape(), &[10_000]);
564    }
565
566    #[test]
567    fn wald_positive() {
568        let mut rng = default_rng_seeded(42);
569        let arr = rng.wald(1.0, 1.0, 10_000).unwrap();
570        for &v in arr.as_slice().unwrap() {
571            assert!(v > 0.0, "wald value {v} must be positive");
572        }
573    }
574
575    #[test]
576    fn standard_cauchy_basic() {
577        let mut rng = default_rng_seeded(42);
578        let arr = rng.standard_cauchy(10_000).unwrap();
579        assert_eq!(arr.shape(), &[10_000]);
580        // Cauchy has no mean or variance, just check it runs
581    }
582
583    #[test]
584    fn bad_params() {
585        let mut rng = default_rng_seeded(42);
586        assert!(rng.laplace(0.0, 0.0, 10).is_err());
587        assert!(rng.logistic(0.0, -1.0, 10).is_err());
588        assert!(rng.rayleigh(0.0, 10).is_err());
589        assert!(rng.weibull(0.0, 10).is_err());
590        assert!(rng.pareto(0.0, 10).is_err());
591        assert!(rng.gumbel(0.0, 0.0, 10).is_err());
592        assert!(rng.power(0.0, 10).is_err());
593        assert!(rng.triangular(5.0, 3.0, 1.0, 10).is_err());
594        assert!(rng.vonmises(0.0, -1.0, 10).is_err());
595        assert!(rng.wald(0.0, 1.0, 10).is_err());
596        assert!(rng.wald(1.0, 0.0, 10).is_err());
597    }
598}
ferray_random/distributions/misc_continuous.rs

ferray_random/distributions/
misc_continuous.rs
