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

ferray_random/distributions/
misc_continuous.rs
