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            loc - scale * u.signum() * (1.0 - 2.0 * u.abs()).ln()
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 loc + scale * (u / (1.0 - u)).ln();
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() - 1.0
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 loc - scale * (-u.ln()).ln();
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.next_f64() * std::f64::consts::TAU - std::f64::consts::PI)
308            });
309            return vec_to_array_f64(data, &shape_vec);
310        }
311
312        // Best & Fisher algorithm
313        let tau = 1.0 + (1.0 + 4.0 * kappa * kappa).sqrt();
314        let rho = (tau - (2.0 * tau).sqrt()) / (2.0 * kappa);
315        let r = (1.0 + rho * rho) / (2.0 * rho);
316
317        let data = generate_vec(self, n, |bg| {
318            loop {
319                let u1 = bg.next_f64();
320                let z = (std::f64::consts::TAU * u1 - std::f64::consts::PI).cos();
321                let f_val = (1.0 + r * z) / (r + z);
322                let c = kappa * (r - f_val);
323                let u2 = bg.next_f64();
324                if u2 < c * (2.0 - c) || u2 <= c * (-c).exp() {
325                    let u3 = bg.next_f64();
326                    let theta = if u3 > 0.5 {
327                        mu + f_val.acos()
328                    } else {
329                        mu - f_val.acos()
330                    };
331                    return theta;
332                }
333            }
334        });
335        vec_to_array_f64(data, &shape_vec)
336    }
337
338    /// Generate an array of Wald (inverse Gaussian) distributed variates.
339    ///
340    /// # Arguments
341    /// * `mean` - Mean of the distribution, must be positive.
342    /// * `scale` - Scale parameter (lambda), must be positive.
343    /// * `size` - Number of values to generate.
344    ///
345    /// # Errors
346    /// Returns `FerrayError::InvalidValue` if `mean <= 0`, `scale <= 0`, or `size` is zero.
347    pub fn wald(
348        &mut self,
349        mean: f64,
350        scale: f64,
351        size: impl IntoShape,
352    ) -> Result<Array<f64, IxDyn>, FerrayError> {
353        if mean <= 0.0 {
354            return Err(FerrayError::invalid_value(format!(
355                "mean must be positive, got {mean}"
356            )));
357        }
358        if scale <= 0.0 {
359            return Err(FerrayError::invalid_value(format!(
360                "scale must be positive, got {scale}"
361            )));
362        }
363        // Michael, Schucany & Haas algorithm
364        let shape_vec = size.into_shape()?;
365        let n = shape_size(&shape_vec);
366        let data = generate_vec(self, n, |bg| {
367            let z = standard_normal_single(bg);
368            let v = z * z;
369            let mu = mean;
370            let lam = scale;
371            let x = mu + (mu * mu * v) / (2.0 * lam)
372                - (mu / (2.0 * lam)) * (4.0 * mu * lam * v + mu * mu * v * v).sqrt();
373            let u = bg.next_f64();
374            if u <= mu / (mu + x) { x } else { mu * mu / x }
375        });
376        vec_to_array_f64(data, &shape_vec)
377    }
378
379    /// Generate an array of standard Cauchy distributed variates.
380    ///
381    /// Uses the inverse CDF: tan(pi * (u - 0.5)).
382    ///
383    /// # Arguments
384    /// * `size` - Number of values to generate.
385    ///
386    /// # Errors
387    /// Returns `FerrayError::InvalidValue` if `size` is zero.
388    pub fn standard_cauchy(
389        &mut self,
390        size: impl IntoShape,
391    ) -> Result<Array<f64, IxDyn>, FerrayError> {
392        let shape_vec = size.into_shape()?;
393        let n = shape_size(&shape_vec);
394        let data = generate_vec(self, n, |bg| {
395            loop {
396                let u = bg.next_f64();
397                // Avoid u = 0.5 exactly (tan(0) = 0, which is fine, but avoid edge)
398                if (u - 0.5).abs() > f64::EPSILON {
399                    return (std::f64::consts::PI * (u - 0.5)).tan();
400                }
401            }
402        });
403        vec_to_array_f64(data, &shape_vec)
404    }
405}
406
407#[cfg(test)]
408mod tests {
409    use crate::default_rng_seeded;
410
411    #[test]
412    fn laplace_mean_variance() {
413        let mut rng = default_rng_seeded(42);
414        let n = 100_000;
415        let loc = 2.0;
416        let scale = 3.0;
417        let arr = rng.laplace(loc, scale, n).unwrap();
418        let slice = arr.as_slice().unwrap();
419        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
420        let var: f64 = slice.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / n as f64;
421        // Laplace(loc, scale): mean=loc, var=2*scale^2
422        let expected_var = 2.0 * scale * scale;
423        let se = (expected_var / n as f64).sqrt();
424        assert!(
425            (mean - loc).abs() < 3.0 * se,
426            "laplace mean {mean} too far from {loc}"
427        );
428        assert!(
429            (var - expected_var).abs() / expected_var < 0.05,
430            "laplace variance {var} too far from {expected_var}"
431        );
432    }
433
434    #[test]
435    fn logistic_mean() {
436        let mut rng = default_rng_seeded(42);
437        let n = 100_000;
438        let loc = 1.0;
439        let scale = 2.0;
440        let arr = rng.logistic(loc, scale, n).unwrap();
441        let slice = arr.as_slice().unwrap();
442        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
443        let expected_var = (std::f64::consts::PI * scale).powi(2) / 3.0;
444        let se = (expected_var / n as f64).sqrt();
445        assert!(
446            (mean - loc).abs() < 3.0 * se,
447            "logistic mean {mean} too far from {loc}"
448        );
449    }
450
451    #[test]
452    fn rayleigh_positive() {
453        let mut rng = default_rng_seeded(42);
454        let arr = rng.rayleigh(2.0, 10_000).unwrap();
455        for &v in arr.as_slice().unwrap() {
456            assert!(v > 0.0);
457        }
458    }
459
460    #[test]
461    fn rayleigh_mean() {
462        let mut rng = default_rng_seeded(42);
463        let n = 100_000;
464        let scale = 2.0;
465        let arr = rng.rayleigh(scale, n).unwrap();
466        let slice = arr.as_slice().unwrap();
467        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
468        // Rayleigh(sigma): mean = sigma * sqrt(pi/2)
469        let expected_mean = scale * (std::f64::consts::FRAC_PI_2).sqrt();
470        let expected_var = (4.0 - std::f64::consts::PI) / 2.0 * scale * scale;
471        let se = (expected_var / n as f64).sqrt();
472        assert!(
473            (mean - expected_mean).abs() < 3.0 * se,
474            "rayleigh mean {mean} too far from {expected_mean}"
475        );
476    }
477
478    #[test]
479    fn weibull_positive() {
480        let mut rng = default_rng_seeded(42);
481        let arr = rng.weibull(1.5, 10_000).unwrap();
482        for &v in arr.as_slice().unwrap() {
483            assert!(v > 0.0);
484        }
485    }
486
487    #[test]
488    fn pareto_ge_zero() {
489        let mut rng = default_rng_seeded(42);
490        let arr = rng.pareto(3.0, 10_000).unwrap();
491        for &v in arr.as_slice().unwrap() {
492            assert!(v >= 0.0, "pareto value {v} < 0");
493        }
494    }
495
496    #[test]
497    fn gumbel_mean() {
498        let mut rng = default_rng_seeded(42);
499        let n = 100_000;
500        let loc = 1.0;
501        let scale = 2.0;
502        let arr = rng.gumbel(loc, scale, n).unwrap();
503        let slice = arr.as_slice().unwrap();
504        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
505        // Gumbel: mean = loc + scale * gamma_euler (0.5772...)
506        let euler = 0.5772156649015329;
507        let expected_mean = loc + scale * euler;
508        let expected_var = std::f64::consts::PI * std::f64::consts::PI * scale * scale / 6.0;
509        let se = (expected_var / n as f64).sqrt();
510        assert!(
511            (mean - expected_mean).abs() < 3.0 * se,
512            "gumbel mean {mean} too far from {expected_mean}"
513        );
514    }
515
516    #[test]
517    fn power_in_range() {
518        let mut rng = default_rng_seeded(42);
519        let arr = rng.power(2.0, 10_000).unwrap();
520        for &v in arr.as_slice().unwrap() {
521            assert!((0.0..=1.0).contains(&v), "power value {v} out of [0,1]");
522        }
523    }
524
525    #[test]
526    fn triangular_in_range() {
527        let mut rng = default_rng_seeded(42);
528        let arr = rng.triangular(1.0, 3.0, 5.0, 10_000).unwrap();
529        for &v in arr.as_slice().unwrap() {
530            assert!(
531                (1.0..=5.0).contains(&v),
532                "triangular value {v} out of [1,5]"
533            );
534        }
535    }
536
537    #[test]
538    fn triangular_mean() {
539        let mut rng = default_rng_seeded(42);
540        let n = 100_000;
541        let (left, mode, right) = (1.0, 3.0, 5.0);
542        let arr = rng.triangular(left, mode, right, n).unwrap();
543        let slice = arr.as_slice().unwrap();
544        let mean: f64 = slice.iter().sum::<f64>() / n as f64;
545        let expected_mean = (left + mode + right) / 3.0;
546        assert!(
547            (mean - expected_mean).abs() < 0.05,
548            "triangular mean {mean} too far from {expected_mean}"
549        );
550    }
551
552    #[test]
553    fn vonmises_basic() {
554        let mut rng = default_rng_seeded(42);
555        let arr = rng.vonmises(0.0, 1.0, 10_000).unwrap();
556        assert_eq!(arr.shape(), &[10_000]);
557    }
558
559    #[test]
560    fn wald_positive() {
561        let mut rng = default_rng_seeded(42);
562        let arr = rng.wald(1.0, 1.0, 10_000).unwrap();
563        for &v in arr.as_slice().unwrap() {
564            assert!(v > 0.0, "wald value {v} must be positive");
565        }
566    }
567
568    #[test]
569    fn standard_cauchy_basic() {
570        let mut rng = default_rng_seeded(42);
571        let arr = rng.standard_cauchy(10_000).unwrap();
572        assert_eq!(arr.shape(), &[10_000]);
573        // Cauchy has no mean or variance, just check it runs
574    }
575
576    #[test]
577    fn bad_params() {
578        let mut rng = default_rng_seeded(42);
579        assert!(rng.laplace(0.0, 0.0, 10).is_err());
580        assert!(rng.logistic(0.0, -1.0, 10).is_err());
581        assert!(rng.rayleigh(0.0, 10).is_err());
582        assert!(rng.weibull(0.0, 10).is_err());
583        assert!(rng.pareto(0.0, 10).is_err());
584        assert!(rng.gumbel(0.0, 0.0, 10).is_err());
585        assert!(rng.power(0.0, 10).is_err());
586        assert!(rng.triangular(5.0, 3.0, 1.0, 10).is_err());
587        assert!(rng.vonmises(0.0, -1.0, 10).is_err());
588        assert!(rng.wald(0.0, 1.0, 10).is_err());
589        assert!(rng.wald(1.0, 0.0, 10).is_err());
590    }
591}
ferray_random/distributions/misc_continuous.rs

ferray_random/distributions/
misc_continuous.rs
