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

ferray_random/distributions/
misc_continuous.rs
