scirs2_stats/distributions/laplace.rs
1//! Laplace distribution functions
2//!
3//! This module provides functionality for the Laplace (double exponential) distribution.
4
5use crate::error::{StatsError, StatsResult};
6use crate::sampling::SampleableDistribution;
7use crate::traits::{ContinuousDistribution, Distribution as ScirsDist};
8use scirs2_core::ndarray::Array1;
9use scirs2_core::numeric::{Float, NumCast};
10use scirs2_core::random::{Distribution, Uniform as RandUniform};
11
12/// Laplace distribution structure
13///
14/// The Laplace distribution, also known as the double exponential distribution,
15/// is a continuous probability distribution that resembles a symmetric version of the
16/// exponential distribution, placed back-to-back. It has heavier tails than the normal distribution.
17pub struct Laplace<F: Float> {
18 /// Location parameter (mean, median, and mode of the distribution)
19 pub loc: F,
20 /// Scale parameter (diversity) > 0
21 pub scale: F,
22 /// Random number generator for uniform distribution
23 rand_distr: RandUniform<f64>,
24}
25
26impl<F: Float + NumCast + std::fmt::Display> Laplace<F> {
27 /// Create a new Laplace distribution with given parameters
28 ///
29 /// # Arguments
30 ///
31 /// * `loc` - Location parameter (mean, median, and mode of the distribution)
32 /// * `scale` - Scale parameter (diversity) > 0
33 ///
34 /// # Returns
35 ///
36 /// * A new Laplace distribution instance
37 ///
38 /// # Examples
39 ///
40 /// ```
41 /// use scirs2_stats::distributions::laplace::Laplace;
42 ///
43 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
44 /// ```
45 pub fn new(loc: F, scale: F) -> StatsResult<Self> {
46 // Validate parameters
47 if scale <= F::zero() {
48 return Err(StatsError::DomainError(
49 "Scale parameter must be positive".to_string(),
50 ));
51 }
52
53 // Create RNG for uniform distribution in [0, 1)
54 let rand_distr = match RandUniform::new(0.0, 1.0) {
55 Ok(distr) => distr,
56 Err(_) => {
57 return Err(StatsError::ComputationError(
58 "Failed to create uniform distribution for sampling".to_string(),
59 ))
60 }
61 };
62
63 Ok(Laplace {
64 loc,
65 scale,
66 rand_distr,
67 })
68 }
69
70 /// Calculate the probability density function (PDF) at a given point
71 ///
72 /// # Arguments
73 ///
74 /// * `x` - The point at which to evaluate the PDF
75 ///
76 /// # Returns
77 ///
78 /// * The value of the PDF at the given point
79 ///
80 /// # Examples
81 ///
82 /// ```
83 /// use scirs2_stats::distributions::laplace::Laplace;
84 ///
85 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
86 /// let pdf_at_zero = laplace.pdf(0.0);
87 /// assert!((pdf_at_zero - 0.5).abs() < 1e-7);
88 /// ```
89 pub fn pdf(&self, x: F) -> F {
90 let half = F::from(0.5).expect("Failed to convert constant to float");
91 let abs_value = if x >= self.loc {
92 x - self.loc
93 } else {
94 self.loc - x
95 };
96
97 // PDF = (1/(2*scale)) * exp(-|x-loc|/scale)
98 half / self.scale * (-abs_value / self.scale).exp()
99 }
100
101 /// Calculate the cumulative distribution function (CDF) at a given point
102 ///
103 /// # Arguments
104 ///
105 /// * `x` - The point at which to evaluate the CDF
106 ///
107 /// # Returns
108 ///
109 /// * The value of the CDF at the given point
110 ///
111 /// # Examples
112 ///
113 /// ```
114 /// use scirs2_stats::distributions::laplace::Laplace;
115 ///
116 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
117 /// let cdf_at_zero = laplace.cdf(0.0);
118 /// assert!((cdf_at_zero - 0.5).abs() < 1e-7);
119 /// ```
120 pub fn cdf(&self, x: F) -> F {
121 let half = F::from(0.5).expect("Failed to convert constant to float");
122
123 if x < self.loc {
124 // CDF = (1/2) * exp((x-loc)/scale)
125 half * ((x - self.loc) / self.scale).exp()
126 } else {
127 // CDF = 1 - (1/2) * exp(-(x-loc)/scale)
128 F::one() - half * (-(x - self.loc) / self.scale).exp()
129 }
130 }
131
132 /// Inverse of the cumulative distribution function (quantile function)
133 ///
134 /// # Arguments
135 ///
136 /// * `p` - Probability value (between 0 and 1)
137 ///
138 /// # Returns
139 ///
140 /// * The value x such that CDF(x) = p
141 ///
142 /// # Examples
143 ///
144 /// ```
145 /// use scirs2_stats::distributions::laplace::Laplace;
146 ///
147 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
148 /// let x = laplace.ppf(0.75).expect("Operation failed");
149 /// assert!((x - 0.693147).abs() < 1e-6);
150 /// ```
151 pub fn ppf(&self, p: F) -> StatsResult<F> {
152 if p < F::zero() || p > F::one() {
153 return Err(StatsError::DomainError(
154 "Probability must be between 0 and 1".to_string(),
155 ));
156 }
157
158 let half = F::from(0.5).expect("Failed to convert constant to float");
159
160 let quantile = if p < half {
161 // Q(p) = loc + scale * ln(2p)
162 self.loc + self.scale * (p + p).ln()
163 } else {
164 // Q(p) = loc - scale * ln(2(1-p))
165 self.loc - self.scale * ((F::one() - p) + (F::one() - p)).ln()
166 };
167
168 Ok(quantile)
169 }
170
171 /// Generate random samples from the distribution
172 ///
173 /// # Arguments
174 ///
175 /// * `size` - Number of samples to generate
176 ///
177 /// # Returns
178 ///
179 /// * Vector of random samples
180 ///
181 /// # Examples
182 ///
183 /// ```
184 /// use scirs2_stats::distributions::laplace::Laplace;
185 ///
186 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
187 /// let samples = laplace.rvs_vec(10).expect("Operation failed");
188 /// assert_eq!(samples.len(), 10);
189 /// ```
190 pub fn rvs_vec(&self, size: usize) -> StatsResult<Vec<F>> {
191 let mut rng = scirs2_core::random::thread_rng();
192 let mut samples = Vec::with_capacity(size);
193
194 for _ in 0..size {
195 // Generate uniform random number in [0, 1)
196 let u = self.rand_distr.sample(&mut rng);
197 let u_f = F::from(u).expect("Failed to convert to float");
198
199 // Apply inverse CDF transform
200 let sample = match self.ppf(u_f) {
201 Ok(s) => s,
202 Err(_) => continue, // Skip invalid samples
203 };
204
205 samples.push(sample);
206 }
207
208 // Ensure we have exactly 'size' samples
209 while samples.len() < size {
210 let u = self.rand_distr.sample(&mut rng);
211 let u_f = F::from(u).expect("Failed to convert to float");
212
213 let sample = match self.ppf(u_f) {
214 Ok(s) => s,
215 Err(_) => continue,
216 };
217
218 samples.push(sample);
219 }
220
221 Ok(samples)
222 }
223
224 /// Generate random samples from the distribution
225 ///
226 /// # Arguments
227 ///
228 /// * `size` - Number of samples to generate
229 ///
230 /// # Returns
231 ///
232 /// * Array of random samples
233 ///
234 /// # Examples
235 ///
236 /// ```
237 /// use scirs2_stats::distributions::laplace::Laplace;
238 ///
239 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
240 /// let samples = laplace.rvs(10).expect("Operation failed");
241 /// assert_eq!(samples.len(), 10);
242 /// ```
243 pub fn rvs(&self, size: usize) -> StatsResult<Array1<F>> {
244 let samples_vec = self.rvs_vec(size)?;
245 Ok(Array1::from(samples_vec))
246 }
247
248 /// Calculate the mean of the distribution
249 ///
250 /// # Returns
251 ///
252 /// * The mean of the distribution
253 ///
254 /// # Examples
255 ///
256 /// ```
257 /// use scirs2_stats::distributions::laplace::Laplace;
258 ///
259 /// let laplace = Laplace::new(2.0f64, 1.0).expect("Operation failed");
260 /// let mean = laplace.mean();
261 /// assert_eq!(mean, 2.0);
262 /// ```
263 pub fn mean(&self) -> F {
264 // Mean = loc
265 self.loc
266 }
267
268 /// Calculate the variance of the distribution
269 ///
270 /// # Returns
271 ///
272 /// * The variance of the distribution
273 ///
274 /// # Examples
275 ///
276 /// ```
277 /// use scirs2_stats::distributions::laplace::Laplace;
278 ///
279 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
280 /// let variance = laplace.var();
281 /// assert!((variance - 2.0).abs() < 1e-7);
282 /// ```
283 pub fn var(&self) -> F {
284 // Variance = 2 * scale^2
285 let two = F::from(2.0).expect("Failed to convert constant to float");
286 two * self.scale * self.scale
287 }
288
289 /// Calculate the standard deviation of the distribution
290 ///
291 /// # Returns
292 ///
293 /// * The standard deviation of the distribution
294 ///
295 /// # Examples
296 ///
297 /// ```
298 /// use scirs2_stats::distributions::laplace::Laplace;
299 ///
300 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
301 /// let std_dev = laplace.std();
302 /// assert!((std_dev - 1.414213).abs() < 1e-6);
303 /// ```
304 pub fn std(&self) -> F {
305 // Std = sqrt(var) = sqrt(2) * scale
306 self.var().sqrt()
307 }
308
309 /// Calculate the median of the distribution
310 ///
311 /// # Returns
312 ///
313 /// * The median of the distribution
314 ///
315 /// # Examples
316 ///
317 /// ```
318 /// use scirs2_stats::distributions::laplace::Laplace;
319 ///
320 /// let laplace = Laplace::new(2.0f64, 1.0).expect("Operation failed");
321 /// let median = laplace.median();
322 /// assert_eq!(median, 2.0);
323 /// ```
324 pub fn median(&self) -> F {
325 // Median = loc
326 self.loc
327 }
328
329 /// Calculate the mode of the distribution
330 ///
331 /// # Returns
332 ///
333 /// * The mode of the distribution
334 ///
335 /// # Examples
336 ///
337 /// ```
338 /// use scirs2_stats::distributions::laplace::Laplace;
339 ///
340 /// let laplace = Laplace::new(2.0f64, 1.0).expect("Operation failed");
341 /// let mode = laplace.mode();
342 /// assert_eq!(mode, 2.0);
343 /// ```
344 pub fn mode(&self) -> F {
345 // Mode = loc
346 self.loc
347 }
348
349 /// Calculate the skewness of the distribution
350 ///
351 /// # Returns
352 ///
353 /// * The skewness (which is always 0 for the Laplace distribution)
354 ///
355 /// # Examples
356 ///
357 /// ```
358 /// use scirs2_stats::distributions::laplace::Laplace;
359 ///
360 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
361 /// let skewness = laplace.skewness();
362 /// assert_eq!(skewness, 0.0);
363 /// ```
364 pub fn skewness(&self) -> F {
365 // Skewness = 0 (symmetric distribution)
366 F::zero()
367 }
368
369 /// Calculate the kurtosis of the distribution
370 ///
371 /// # Returns
372 ///
373 /// * The kurtosis of the distribution
374 ///
375 /// # Examples
376 ///
377 /// ```
378 /// use scirs2_stats::distributions::laplace::Laplace;
379 ///
380 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
381 /// let kurtosis = laplace.kurtosis();
382 /// assert!((kurtosis - 3.0).abs() < 1e-7);
383 /// ```
384 pub fn kurtosis(&self) -> F {
385 // Excess kurtosis = 3 (higher than normal distribution's 0)
386 F::from(3.0).expect("Failed to convert constant to float")
387 }
388
389 /// Calculate the entropy of the distribution
390 ///
391 /// # Returns
392 ///
393 /// * The entropy value
394 ///
395 /// # Examples
396 ///
397 /// ```
398 /// use scirs2_stats::distributions::laplace::Laplace;
399 ///
400 /// let laplace = Laplace::new(0.0f64, 1.0).expect("Operation failed");
401 /// let entropy = laplace.entropy();
402 /// assert!((entropy - 1.693147).abs() < 1e-6);
403 /// ```
404 pub fn entropy(&self) -> F {
405 // Entropy = 1 + ln(2*scale)
406 let one = F::one();
407 let two = F::from(2.0).expect("Failed to convert constant to float");
408 one + (two * self.scale).ln()
409 }
410}
411
412/// Create a Laplace distribution with the given parameters.
413///
414/// This is a convenience function to create a Laplace distribution with
415/// the given location and scale parameters.
416///
417/// # Arguments
418///
419/// * `loc` - Location parameter (mean, median, and mode of the distribution)
420/// * `scale` - Scale parameter (diversity) > 0
421///
422/// # Returns
423///
424/// * A Laplace distribution object
425///
426/// # Examples
427///
428/// ```
429/// use scirs2_stats::distributions::laplace;
430///
431/// let l = laplace::laplace(0.0f64, 1.0).expect("Operation failed");
432/// let pdf_at_zero = l.pdf(0.0);
433/// assert!((pdf_at_zero - 0.5).abs() < 1e-7);
434/// ```
435#[allow(dead_code)]
436pub fn laplace<F>(loc: F, scale: F) -> StatsResult<Laplace<F>>
437where
438 F: Float + NumCast + std::fmt::Display,
439{
440 Laplace::new(loc, scale)
441}
442
443/// Implementation of SampleableDistribution for Laplace
444impl<F: Float + NumCast + std::fmt::Display> SampleableDistribution<F> for Laplace<F> {
445 fn rvs(&self, size: usize) -> StatsResult<Vec<F>> {
446 self.rvs_vec(size)
447 }
448}
449
450/// Implementation of Distribution trait for Laplace
451impl<F: Float + NumCast + std::fmt::Display> ScirsDist<F> for Laplace<F> {
452 fn mean(&self) -> F {
453 self.mean()
454 }
455
456 fn var(&self) -> F {
457 self.var()
458 }
459
460 fn std(&self) -> F {
461 self.std()
462 }
463
464 fn rvs(&self, size: usize) -> StatsResult<Array1<F>> {
465 self.rvs(size)
466 }
467
468 fn entropy(&self) -> F {
469 self.entropy()
470 }
471}
472
473/// Implementation of ContinuousDistribution trait for Laplace
474impl<F: Float + NumCast + std::fmt::Display> ContinuousDistribution<F> for Laplace<F> {
475 fn pdf(&self, x: F) -> F {
476 self.pdf(x)
477 }
478
479 fn cdf(&self, x: F) -> F {
480 self.cdf(x)
481 }
482
483 fn ppf(&self, p: F) -> StatsResult<F> {
484 self.ppf(p)
485 }
486}
487
488#[cfg(test)]
489mod tests {
490 use super::*;
491 use approx::assert_relative_eq;
492
493 #[test]
494 fn test_laplace_creation() {
495 // Standard Laplace (loc=0, scale=1)
496 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
497 assert_eq!(laplace.loc, 0.0);
498 assert_eq!(laplace.scale, 1.0);
499
500 // Custom Laplace
501 let custom = Laplace::new(-2.0, 0.5).expect("Operation failed");
502 assert_eq!(custom.loc, -2.0);
503 assert_eq!(custom.scale, 0.5);
504
505 // Error case: non-positive scale
506 assert!(Laplace::<f64>::new(0.0, 0.0).is_err());
507 assert!(Laplace::<f64>::new(0.0, -1.0).is_err());
508 }
509
510 #[test]
511 fn test_laplace_pdf() {
512 // Standard Laplace (loc=0, scale=1)
513 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
514
515 // PDF at x = 0 should be 1/(2*1) = 0.5
516 let pdf_at_zero = laplace.pdf(0.0);
517 assert_relative_eq!(pdf_at_zero, 0.5, epsilon = 1e-7);
518
519 // PDF at x = 1 should be 0.5 * exp(-1/1) = 0.5 * exp(-1) = 0.5 * 0.36787944 ≈ 0.1839397
520 let pdf_at_one = laplace.pdf(1.0);
521 assert_relative_eq!(pdf_at_one, 0.1839397, epsilon = 1e-7);
522
523 // PDF at x = -1 should be same as x = 1 due to symmetry
524 let pdf_at_neg_one = laplace.pdf(-1.0);
525 assert_relative_eq!(pdf_at_neg_one, pdf_at_one, epsilon = 1e-10);
526
527 // Custom Laplace with loc=-2, scale=0.5
528 let custom = Laplace::new(-2.0, 0.5).expect("Operation failed");
529
530 // PDF at x = -2 should be 1/(2*0.5) = 1
531 let pdf_at_loc = custom.pdf(-2.0);
532 assert_relative_eq!(pdf_at_loc, 1.0, epsilon = 1e-7);
533
534 // PDF at x = -1.5 should be 0.5/0.5 * exp(-|-1.5-(-2)|/0.5) = exp(-0.5/0.5) = exp(-1) ≈ 0.36787944
535 let pdf_at_custom = custom.pdf(-1.5);
536 assert_relative_eq!(pdf_at_custom, 0.36787944, epsilon = 1e-7);
537 }
538
539 #[test]
540 fn test_laplace_cdf() {
541 // Standard Laplace (loc=0, scale=1)
542 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
543
544 // CDF at x = 0 should be 0.5
545 let cdf_at_zero = laplace.cdf(0.0);
546 assert_relative_eq!(cdf_at_zero, 0.5, epsilon = 1e-10);
547
548 // CDF at x = 1 should be 1 - 0.5*exp(-1/1) = 1 - 0.5*exp(-1) ≈ 1 - 0.5*0.36787944 ≈ 0.8160603
549 let cdf_at_one = laplace.cdf(1.0);
550 assert_relative_eq!(cdf_at_one, 0.8160603, epsilon = 1e-7);
551
552 // CDF at x = -1 should be 0.5*exp((-1-0)/1) = 0.5*exp(-1) ≈ 0.5*0.36787944 ≈ 0.1839397
553 let cdf_at_neg_one = laplace.cdf(-1.0);
554 assert_relative_eq!(cdf_at_neg_one, 0.1839397, epsilon = 1e-7);
555
556 // Custom Laplace with loc=-2, scale=0.5
557 let custom = Laplace::new(-2.0, 0.5).expect("Operation failed");
558
559 // CDF at x = -2 should be 0.5
560 let cdf_at_loc = custom.cdf(-2.0);
561 assert_relative_eq!(cdf_at_loc, 0.5, epsilon = 1e-10);
562
563 // CDF at x = -1.5 should be 1 - 0.5*exp(-(-1.5-(-2))/0.5) = 1 - 0.5*exp(-0.5/0.5) = 1 - 0.5*exp(-1) ≈ 0.8160603
564 let cdf_at_custom = custom.cdf(-1.5);
565 assert_relative_eq!(cdf_at_custom, 0.8160603, epsilon = 1e-7);
566 }
567
568 #[test]
569 fn test_laplace_ppf() {
570 // Standard Laplace (loc=0, scale=1)
571 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
572
573 // PPF at p = 0.5 should be 0
574 let ppf_at_half = laplace.ppf(0.5).expect("Operation failed");
575 assert_relative_eq!(ppf_at_half, 0.0, epsilon = 1e-10);
576
577 // PPF at p = 0.75 should be log(2*0.75) ≈ log(1.5) ≈ 0.693147
578 let ppf_at_75 = laplace.ppf(0.75).expect("Operation failed");
579 assert_relative_eq!(ppf_at_75, 0.693147, epsilon = 1e-6);
580
581 // PPF at p = 0.25 should be -log(2*0.75) ≈ -log(1.5) ≈ -0.693147
582 let ppf_at_25 = laplace.ppf(0.25).expect("Operation failed");
583 assert_relative_eq!(ppf_at_25, -0.693147, epsilon = 1e-6);
584
585 // Custom Laplace with loc=-2, scale=0.5
586 let custom = Laplace::new(-2.0, 0.5).expect("Operation failed");
587
588 // PPF at p = 0.5 should be -2.0
589 let ppf_at_half_custom = custom.ppf(0.5).expect("Operation failed");
590 assert_relative_eq!(ppf_at_half_custom, -2.0, epsilon = 1e-10);
591
592 // PPF at p = 0.75 should be -2.0 + 0.5*log(1.5) ≈ -2.0 + 0.5*0.693147 ≈ -1.653426
593 let ppf_at_75_custom = custom.ppf(0.75).expect("Operation failed");
594 assert_relative_eq!(ppf_at_75_custom, -1.653426, epsilon = 1e-6);
595
596 // Error cases
597 assert!(laplace.ppf(-0.1).is_err());
598 assert!(laplace.ppf(1.1).is_err());
599 }
600
601 #[test]
602 fn test_laplace_properties() {
603 // Standard Laplace (loc=0, scale=1)
604 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
605
606 // Mean = loc = 0
607 let mean = laplace.mean();
608 assert_eq!(mean, 0.0);
609
610 // Variance = 2 * scale^2 = 2 * 1^2 = 2
611 let variance = laplace.var();
612 assert_eq!(variance, 2.0);
613
614 // Std = sqrt(variance) = sqrt(2) ≈ 1.414213
615 let std_dev = laplace.std();
616 assert_relative_eq!(std_dev, 1.414213, epsilon = 1e-6);
617
618 // Median = loc = 0
619 let median = laplace.median();
620 assert_eq!(median, 0.0);
621
622 // Mode = loc = 0
623 let mode = laplace.mode();
624 assert_eq!(mode, 0.0);
625
626 // Skewness = 0 (symmetric)
627 let skewness = laplace.skewness();
628 assert_eq!(skewness, 0.0);
629
630 // Kurtosis = 3
631 let kurtosis = laplace.kurtosis();
632 assert_eq!(kurtosis, 3.0);
633
634 // Entropy = 1 + ln(2*scale) = 1 + ln(2) ≈ 1.693147
635 let entropy = laplace.entropy();
636 assert_relative_eq!(entropy, 1.693147, epsilon = 1e-6);
637
638 // Custom Laplace with loc=-2, scale=0.5
639 let custom = Laplace::new(-2.0, 0.5).expect("Operation failed");
640
641 // Mean = loc = -2
642 let mean_custom = custom.mean();
643 assert_eq!(mean_custom, -2.0);
644
645 // Variance = 2 * scale^2 = 2 * 0.5^2 = 0.5
646 let variance_custom = custom.var();
647 assert_eq!(variance_custom, 0.5);
648
649 // Std = sqrt(variance) = sqrt(0.5) ≈ 0.707107
650 let std_dev_custom = custom.std();
651 assert_relative_eq!(std_dev_custom, 0.707107, epsilon = 1e-6);
652
653 // Entropy = 1 + ln(2*scale) = 1 + ln(2*0.5) = 1 + ln(1) = 1
654 let entropy_custom = custom.entropy();
655 assert_relative_eq!(entropy_custom, 1.0, epsilon = 1e-10);
656 }
657
658 #[test]
659 fn test_laplace_rvs() {
660 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
661
662 // Generate samples
663 let samples = laplace.rvs(100).expect("Operation failed");
664
665 // Check the number of samples
666 assert_eq!(samples.len(), 100);
667
668 // Calculate sample mean and check it's reasonably close to loc = 0
669 // (with large enough samples)
670 let sum: f64 = samples.iter().sum();
671 let mean = sum / samples.len() as f64;
672
673 // The mean could be off due to randomness, but should be within a reasonable range
674 assert!(mean.abs() < 0.5);
675 }
676
677 #[test]
678 fn test_laplace_inverse_cdf() {
679 // Test that cdf(ppf(p)) == p and ppf(cdf(x)) == x
680 let laplace = Laplace::new(0.0, 1.0).expect("Operation failed");
681
682 // Test various probability values
683 let probabilities = [0.1, 0.25, 0.5, 0.75, 0.9];
684 for &p in &probabilities {
685 let x = laplace.ppf(p).expect("Operation failed");
686 let p_back = laplace.cdf(x);
687 assert_relative_eq!(p_back, p, epsilon = 1e-7);
688 }
689
690 // Test various x values
691 let x_values = [-3.0, -1.0, 0.0, 1.0, 3.0];
692 for &x in &x_values {
693 let p = laplace.cdf(x);
694 let x_back = laplace.ppf(p).expect("Operation failed");
695 assert_relative_eq!(x_back, x, epsilon = 1e-7);
696 }
697 }
698}