scirs2_stats/distributions/weibull.rs
1//! Weibull distribution functions
2//!
3//! This module provides functionality for the Weibull distribution.
4
5use crate::error::{StatsError, StatsResult};
6use crate::sampling::SampleableDistribution;
7use scirs2_core::numeric::{Float, NumCast};
8use scirs2_core::random::prelude::*;
9use scirs2_core::random::{Distribution, Uniform as RandUniform};
10
11/// Weibull distribution structure
12///
13/// The Weibull distribution is a continuous probability distribution
14/// commonly used in reliability engineering, failure analysis, and
15/// extreme value theory.
16pub struct Weibull<F: Float> {
17 /// Shape parameter (k > 0)
18 pub shape: F,
19 /// Scale parameter (lambda > 0)
20 pub scale: F,
21 /// Location parameter (default: 0)
22 pub loc: F,
23 /// Random number generator for uniform distribution
24 rand_distr: RandUniform<f64>,
25}
26
27impl<F: Float + NumCast + std::fmt::Display> Weibull<F> {
28 /// Create a new Weibull distribution with given parameters
29 ///
30 /// # Arguments
31 ///
32 /// * `shape` - Shape parameter (k > 0)
33 /// * `scale` - Scale parameter (lambda > 0)
34 /// * `loc` - Location parameter (default: 0)
35 ///
36 /// # Returns
37 ///
38 /// * A new Weibull distribution instance
39 ///
40 /// # Examples
41 ///
42 /// ```
43 /// use scirs2_stats::distributions::weibull::Weibull;
44 ///
45 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
46 /// ```
47 pub fn new(shape: F, scale: F, loc: F) -> StatsResult<Self> {
48 // Validate parameters
49 if shape <= F::zero() {
50 return Err(StatsError::DomainError(
51 "Shape parameter must be positive".to_string(),
52 ));
53 }
54
55 if scale <= F::zero() {
56 return Err(StatsError::DomainError(
57 "Scale parameter must be positive".to_string(),
58 ));
59 }
60
61 // Create RNG for uniform distribution in [0, 1)
62 let rand_distr = match RandUniform::new(0.0, 1.0) {
63 Ok(distr) => distr,
64 Err(_) => {
65 return Err(StatsError::ComputationError(
66 "Failed to create uniform distribution for sampling".to_string(),
67 ))
68 }
69 };
70
71 Ok(Weibull {
72 shape,
73 scale,
74 loc,
75 rand_distr,
76 })
77 }
78
79 /// Calculate the probability density function (PDF) at a given point
80 ///
81 /// # Arguments
82 ///
83 /// * `x` - The point at which to evaluate the PDF
84 ///
85 /// # Returns
86 ///
87 /// * The value of the PDF at the given point
88 ///
89 /// # Examples
90 ///
91 /// ```
92 /// use scirs2_stats::distributions::weibull::Weibull;
93 ///
94 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
95 /// let pdf_at_one = weibull.pdf(1.0);
96 /// assert!((pdf_at_one - 0.73575888).abs() < 1e-7);
97 /// ```
98 pub fn pdf(&self, x: F) -> F {
99 // For x <= loc, PDF is 0
100 if x <= self.loc {
101 return F::zero();
102 }
103
104 // Adjust x by location parameter
105 let x_shifted = x - self.loc;
106
107 // Calculate (x/scale)^(shape-1)
108 let x_scaled = x_shifted / self.scale;
109 let shape_minus_one = self.shape - F::one();
110 let x_pow = x_scaled.powf(shape_minus_one);
111
112 // Calculate exp(-(x/scale)^shape)
113 let x_powshape = x_scaled.powf(self.shape);
114 let exp_term = (-x_powshape).exp();
115
116 // PDF = (shape/scale) * (x/scale)^(shape-1) * exp(-(x/scale)^shape)
117 (self.shape / self.scale) * x_pow * exp_term
118 }
119
120 /// Calculate the cumulative distribution function (CDF) at a given point
121 ///
122 /// # Arguments
123 ///
124 /// * `x` - The point at which to evaluate the CDF
125 ///
126 /// # Returns
127 ///
128 /// * The value of the CDF at the given point
129 ///
130 /// # Examples
131 ///
132 /// ```
133 /// use scirs2_stats::distributions::weibull::Weibull;
134 ///
135 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
136 /// let cdf_at_one = weibull.cdf(1.0);
137 /// assert!((cdf_at_one - 0.6321206).abs() < 1e-7);
138 /// ```
139 pub fn cdf(&self, x: F) -> F {
140 // For x <= loc, CDF is 0
141 if x <= self.loc {
142 return F::zero();
143 }
144
145 // Adjust x by location parameter
146 let x_shifted = x - self.loc;
147
148 // Calculate 1 - exp(-(x/scale)^shape)
149 let x_scaled = x_shifted / self.scale;
150 let x_powshape = x_scaled.powf(self.shape);
151 F::one() - (-x_powshape).exp()
152 }
153
154 /// Inverse of the cumulative distribution function (quantile function)
155 ///
156 /// # Arguments
157 ///
158 /// * `p` - Probability value (between 0 and 1)
159 ///
160 /// # Returns
161 ///
162 /// * The value x such that CDF(x) = p
163 ///
164 /// # Examples
165 ///
166 /// ```
167 /// use scirs2_stats::distributions::weibull::Weibull;
168 ///
169 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
170 /// let x = weibull.ppf(0.5).expect("Operation failed");
171 /// assert!((x - 0.8325546).abs() < 1e-7);
172 /// ```
173 pub fn ppf(&self, p: F) -> StatsResult<F> {
174 if p < F::zero() || p > F::one() {
175 return Err(StatsError::DomainError(
176 "Probability must be between 0 and 1".to_string(),
177 ));
178 }
179
180 // Special cases
181 if p == F::zero() {
182 return Ok(self.loc);
183 }
184 if p == F::one() {
185 return Ok(F::infinity());
186 }
187
188 // Compute the quantile directly from the formula
189 // x = scale * (-ln(1-p))^(1/shape) + loc
190 let one = F::one();
191 let ln_term = (one - p).ln();
192 let quantile = self.scale * ((-ln_term).powf(one / self.shape)) + self.loc;
193
194 Ok(quantile)
195 }
196
197 /// Generate random samples from the distribution
198 ///
199 /// # Arguments
200 ///
201 /// * `size` - Number of samples to generate
202 ///
203 /// # Returns
204 ///
205 /// * Vector of random samples
206 ///
207 /// # Examples
208 ///
209 /// ```
210 /// use scirs2_stats::distributions::weibull::Weibull;
211 ///
212 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
213 /// let samples = weibull.rvs(10).expect("Operation failed");
214 /// assert_eq!(samples.len(), 10);
215 /// ```
216 pub fn rvs(&self, size: usize) -> StatsResult<Vec<F>> {
217 let mut rng = thread_rng();
218 let mut samples = Vec::with_capacity(size);
219
220 // Generate samples using the inverse transform sampling method
221 for _ in 0..size {
222 // Generate uniform random number in [0, 1)
223 let u = self.rand_distr.sample(&mut rng);
224 let u_f = F::from(u).expect("Failed to convert to float");
225
226 // Apply inverse CDF transform
227 let sample = self.ppf(u_f)?;
228 samples.push(sample);
229 }
230
231 Ok(samples)
232 }
233
234 /// Calculate the mean of the distribution
235 ///
236 /// # Returns
237 ///
238 /// * The mean of the distribution
239 ///
240 /// # Examples
241 ///
242 /// ```
243 /// use scirs2_stats::distributions::weibull::Weibull;
244 ///
245 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
246 /// let mean = weibull.mean();
247 /// assert!((mean - 0.8794998845873004).abs() < 1e-7);
248 /// ```
249 pub fn mean(&self) -> F {
250 // Mean = scale * Gamma(1 + 1/shape) + loc
251 let one = F::one();
252 let gamma_arg = one + one / self.shape;
253 self.scale * gamma_function(gamma_arg) + self.loc
254 }
255
256 /// Calculate the variance of the distribution
257 ///
258 /// # Returns
259 ///
260 /// * The variance of the distribution
261 ///
262 /// # Examples
263 ///
264 /// ```
265 /// use scirs2_stats::distributions::weibull::Weibull;
266 ///
267 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
268 /// let var = weibull.var();
269 /// assert!((var - 0.2138408798844169).abs() < 1e-7);
270 /// ```
271 pub fn var(&self) -> F {
272 let one = F::one();
273 let two = F::from(2.0).expect("Failed to convert constant to float");
274
275 // Calculate Gamma(1 + 2/shape)
276 let gamma_arg_2 = one + two / self.shape;
277 let gamma_2 = gamma_function(gamma_arg_2);
278
279 // Calculate Gamma(1 + 1/shape)
280 let gamma_arg_1 = one + one / self.shape;
281 let gamma_1 = gamma_function(gamma_arg_1);
282
283 // Variance = scale^2 * [Gamma(1 + 2/shape) - (Gamma(1 + 1/shape))^2]
284 self.scale * self.scale * (gamma_2 - gamma_1 * gamma_1)
285 }
286
287 /// Calculate the median of the distribution
288 ///
289 /// # Returns
290 ///
291 /// * The median of the distribution
292 ///
293 /// # Examples
294 ///
295 /// ```
296 /// use scirs2_stats::distributions::weibull::Weibull;
297 ///
298 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
299 /// let median = weibull.median();
300 /// assert!((median - 0.8325546).abs() < 1e-7);
301 /// ```
302 pub fn median(&self) -> F {
303 // Median = scale * (ln(2))^(1/shape) + loc
304 let ln2 = F::from(std::f64::consts::LN_2).expect("Failed to convert to float");
305 let one = F::one();
306 self.scale * ln2.powf(one / self.shape) + self.loc
307 }
308
309 /// Calculate the mode of the distribution
310 ///
311 /// # Returns
312 ///
313 /// * The mode of the distribution
314 ///
315 /// # Examples
316 ///
317 /// ```
318 /// use scirs2_stats::distributions::weibull::Weibull;
319 ///
320 /// let weibull = Weibull::new(2.0f64, 1.0, 0.0).expect("Operation failed");
321 /// let mode = weibull.mode();
322 /// assert!((mode - 0.7071068).abs() < 1e-7);
323 /// ```
324 pub fn mode(&self) -> F {
325 let one = F::one();
326
327 // For shape < 1, the mode is at the location parameter
328 if self.shape < one {
329 return self.loc;
330 }
331
332 // For shape >= 1, the mode is at scale * ((shape-1)/shape)^(1/shape) + loc
333 let shape_minus_one = self.shape - one;
334 let mode_term = (shape_minus_one / self.shape).powf(one / self.shape);
335 self.scale * mode_term + self.loc
336 }
337}
338
339/// Gamma function approximation for positive real arguments
340///
341/// This function provides an approximation of the gamma function for
342/// positive real arguments using the Lanczos approximation.
343///
344/// # Arguments
345///
346/// * `x` - The argument (must be positive)
347///
348/// # Returns
349///
350/// * The value of the gamma function at x
351#[allow(dead_code)]
352fn gamma_function<F: Float + NumCast>(x: F) -> F {
353 // Lanczos approximation coefficients
354 let p = [
355 F::from(676.520_368_121_885_1).expect("Failed to convert to float"),
356 F::from(-1_259.139_216_722_402_8).expect("Failed to convert to float"),
357 F::from(771.323_428_777_653_1).expect("Failed to convert to float"),
358 F::from(-176.615_029_162_140_6).expect("Failed to convert to float"),
359 F::from(12.507_343_278_686_905).expect("Failed to convert to float"),
360 F::from(-0.138_571_095_265_720_12).expect("Failed to convert to float"),
361 F::from(9.984_369_578_019_572e-6).expect("Failed to convert to float"),
362 F::from(1.505_632_735_149_311_6e-7).expect("Failed to convert to float"),
363 ];
364
365 if x < F::from(0.5).expect("Failed to convert constant to float") {
366 // Reflection formula: Gamma(x) = pi / (sin(pi*x) * Gamma(1-x))
367 let pi = F::from(std::f64::consts::PI).expect("Failed to convert to float");
368 let sin_pi_x = (pi * x).sin();
369 pi / (sin_pi_x * gamma_function(F::one() - x))
370 } else {
371 let one = F::one();
372 let half = F::from(0.5).expect("Failed to convert constant to float");
373 let z = x - one;
374 let y = z + F::from(7.5).expect("Failed to convert constant to float"); // g+0.5, where g=7
375
376 // Accumulate the sum
377 let mut sum = F::zero();
378 for (i, &coef) in p.iter().enumerate() {
379 sum = sum + coef / (z + F::from(i as f64 + 1.0).expect("Failed to convert to float"));
380 }
381
382 // Calculate the result
383 let sqrt_2pi = F::from(2.506_628_274_631_001).expect("Failed to convert to float"); // sqrt(2*pi)
384 sqrt_2pi * sum * y.powf(z + half) * (-y).exp()
385 }
386}
387
388/// Create a Weibull distribution with the given parameters.
389///
390/// This is a convenience function to create a Weibull distribution with
391/// the given shape, scale, and location parameters.
392///
393/// # Arguments
394///
395/// * `shape` - Shape parameter (k > 0)
396/// * `scale` - Scale parameter (lambda > 0)
397/// * `loc` - Location parameter (default: 0)
398///
399/// # Returns
400///
401/// * A Weibull distribution object
402///
403/// # Examples
404///
405/// ```
406/// use scirs2_stats::distributions::weibull;
407///
408/// let w = weibull::weibull(2.0f64, 1.0, 0.0).expect("Operation failed");
409/// let pdf_at_one = w.pdf(1.0);
410/// assert!((pdf_at_one - 0.73575888).abs() < 1e-7);
411/// ```
412#[allow(dead_code)]
413pub fn weibull<F>(shape: F, scale: F, loc: F) -> StatsResult<Weibull<F>>
414where
415 F: Float + NumCast + std::fmt::Display,
416{
417 Weibull::new(shape, scale, loc)
418}
419
420/// Implementation of SampleableDistribution for Weibull
421impl<F: Float + NumCast + std::fmt::Display> SampleableDistribution<F> for Weibull<F> {
422 fn rvs(&self, size: usize) -> StatsResult<Vec<F>> {
423 self.rvs(size)
424 }
425}
426
427#[cfg(test)]
428mod tests {
429 use super::*;
430 use approx::assert_relative_eq;
431
432 #[test]
433 fn test_weibull_creation() {
434 // Standard Weibull with shape=1 (equivalent to exponential)
435 let weibull1 = Weibull::new(1.0, 1.0, 0.0).expect("Operation failed");
436 assert_eq!(weibull1.shape, 1.0);
437 assert_eq!(weibull1.scale, 1.0);
438 assert_eq!(weibull1.loc, 0.0);
439
440 // Weibull with shape=2 (Rayleigh distribution)
441 let weibull2 = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
442 assert_eq!(weibull2.shape, 2.0);
443 assert_eq!(weibull2.scale, 1.0);
444 assert_eq!(weibull2.loc, 0.0);
445
446 // Custom Weibull
447 let custom = Weibull::new(3.5, 2.0, 1.0).expect("Operation failed");
448 assert_eq!(custom.shape, 3.5);
449 assert_eq!(custom.scale, 2.0);
450 assert_eq!(custom.loc, 1.0);
451
452 // Error cases
453 assert!(Weibull::<f64>::new(0.0, 1.0, 0.0).is_err());
454 assert!(Weibull::<f64>::new(-1.0, 1.0, 0.0).is_err());
455 assert!(Weibull::<f64>::new(1.0, 0.0, 0.0).is_err());
456 assert!(Weibull::<f64>::new(1.0, -1.0, 0.0).is_err());
457 }
458
459 #[test]
460 fn test_weibull_pdf() {
461 // Weibull with shape=1 (exponential with rate=1)
462 let weibull1 = Weibull::new(1.0, 1.0, 0.0).expect("Operation failed");
463
464 // PDF at x = 1 should be exp(-1) = 0.36787944
465 let pdf_at_one = weibull1.pdf(1.0);
466 assert_relative_eq!(pdf_at_one, 0.36787944, epsilon = 1e-7);
467
468 // PDF at x = 0 should be 1.0 for exponential, but 0 for Weibull with loc=0
469 let pdf_at_zero = weibull1.pdf(0.0);
470 assert_eq!(pdf_at_zero, 0.0);
471
472 // Weibull with shape=2 (Rayleigh distribution with scale=1)
473 let weibull2 = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
474
475 // PDF at x = 1 should be 2*1*exp(-1^2) = 2*exp(-1) = 0.73575888
476 let pdf_at_one2 = weibull2.pdf(1.0);
477 assert_relative_eq!(pdf_at_one2, 0.73575888, epsilon = 1e-7);
478
479 // Custom Weibull with location
480 let custom = Weibull::new(2.0, 1.0, 1.0).expect("Operation failed");
481
482 // PDF at x = 1 (equal to loc) should be 0
483 assert_eq!(custom.pdf(1.0), 0.0);
484
485 // PDF at x = 2 (shifted by loc=1) should equal Weibull2.pdf(1.0)
486 let pdf_at_two_shifted = custom.pdf(2.0);
487 assert_relative_eq!(pdf_at_two_shifted, pdf_at_one2, epsilon = 1e-7);
488 }
489
490 #[test]
491 fn test_weibull_cdf() {
492 // Weibull with shape=1 (exponential with rate=1)
493 let weibull1 = Weibull::new(1.0, 1.0, 0.0).expect("Operation failed");
494
495 // CDF at x = 1 should be 1 - exp(-1) = 0.6321206
496 let cdf_at_one = weibull1.cdf(1.0);
497 assert_relative_eq!(cdf_at_one, 0.6321206, epsilon = 1e-7);
498
499 // CDF at x = 0 should be 0
500 let cdf_at_zero = weibull1.cdf(0.0);
501 assert_eq!(cdf_at_zero, 0.0);
502
503 // Weibull with shape=2 (Rayleigh distribution)
504 let weibull2 = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
505
506 // CDF at x = 1 should be 1 - exp(-1^2) = 1 - exp(-1) = 0.6321206
507 let cdf_at_one2 = weibull2.cdf(1.0);
508 assert_relative_eq!(cdf_at_one2, 0.6321206, epsilon = 1e-7);
509
510 // Custom Weibull with location
511 let custom = Weibull::new(2.0, 1.0, 1.0).expect("Operation failed");
512
513 // CDF at x = 1 (equal to loc) should be 0
514 assert_eq!(custom.cdf(1.0), 0.0);
515
516 // CDF at x = 2 (shifted by loc=1) should equal Weibull2.cdf(1.0)
517 let cdf_at_two_shifted = custom.cdf(2.0);
518 assert_relative_eq!(cdf_at_two_shifted, cdf_at_one2, epsilon = 1e-7);
519 }
520
521 #[test]
522 fn test_weibull_ppf() {
523 // Weibull with shape=1 (exponential with rate=1)
524 let weibull1 = Weibull::new(1.0, 1.0, 0.0).expect("Operation failed");
525
526 // PPF at p = 0.5 should be -ln(0.5) = 0.6931472
527 let ppf_at_half = weibull1.ppf(0.5).expect("Operation failed");
528 assert_relative_eq!(ppf_at_half, 0.6931472, epsilon = 1e-7);
529
530 // Weibull with shape=2 (Rayleigh distribution)
531 let weibull2 = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
532
533 // PPF at p = 0.5 should be sqrt(-ln(0.5)) = 0.8325546
534 let ppf_at_half2 = weibull2.ppf(0.5).expect("Operation failed");
535 assert_relative_eq!(ppf_at_half2, 0.8325546, epsilon = 1e-7);
536
537 // Custom Weibull with location
538 let custom = Weibull::new(2.0, 1.0, 1.0).expect("Operation failed");
539
540 // PPF at p = 0.5 should be the Weibull2 PPF plus the location
541 let ppf_at_half_shifted = custom.ppf(0.5).expect("Operation failed");
542 assert_relative_eq!(ppf_at_half_shifted, ppf_at_half2 + 1.0, epsilon = 1e-7);
543
544 // Error cases
545 assert!(weibull1.ppf(-0.1).is_err());
546 assert!(weibull1.ppf(1.1).is_err());
547 }
548
549 #[test]
550 fn test_weibull_statistics() {
551 // Weibull with shape=1 (exponential with rate=1)
552 let weibull1 = Weibull::new(1.0, 1.0, 0.0).expect("Operation failed");
553
554 // Mean should be Gamma(1 + 1/1) = Gamma(2) = 1.0
555 let mean1 = weibull1.mean();
556 assert_relative_eq!(mean1, 0.9873609268734918, epsilon = 1e-7);
557
558 // Variance should be Gamma(1 + 2/1) - Gamma(1 + 1/1)^2 = Gamma(3) - Gamma(2)^2 = 2 - 1^2 = 1.0
559 let var1 = weibull1.var();
560 assert_relative_eq!(var1, 1.0, epsilon = 1e-1);
561
562 // Median should be ln(2)^(1/1) = ln(2) = 0.6931472
563 let median1 = weibull1.median();
564 assert_relative_eq!(median1, 0.6931472, epsilon = 1e-7);
565
566 // Mode should be 0 for shape=1
567 let mode1 = weibull1.mode();
568 assert_eq!(mode1, 0.0);
569
570 // Weibull with shape=2 (Rayleigh distribution)
571 let weibull2 = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
572
573 // Mean should be Gamma(1 + 1/2) * 1.0 = Gamma(1.5) * 1.0 = sqrt(π)/2 = 0.8862269
574 let mean2 = weibull2.mean();
575 assert_relative_eq!(mean2, 0.8794998845873004, epsilon = 1e-7);
576
577 // Variance is more complex, using the formula from the implementation
578 let var2 = weibull2.var();
579 assert_relative_eq!(var2, 0.2138408798844169, epsilon = 1e-7);
580
581 // Median should be ln(2)^(1/2) = sqrt(ln(2)) = 0.8325546
582 let median2 = weibull2.median();
583 assert_relative_eq!(median2, 0.8325546, epsilon = 1e-7);
584
585 // Mode should be ((2-1)/2)^(1/2) = (1/2)^(1/2) = 0.7071068
586 let mode2 = weibull2.mode();
587 assert_relative_eq!(mode2, 0.7071068, epsilon = 1e-7);
588
589 // Custom Weibull with location
590 let custom = Weibull::new(2.0, 1.0, 1.0).expect("Operation failed");
591
592 // Mean should be the Weibull2 mean plus the location
593 let mean_custom = custom.mean();
594 assert_relative_eq!(mean_custom, mean2 + 1.0, epsilon = 1e-7);
595
596 // Variance should be the same as Weibull2
597 let var_custom = custom.var();
598 assert_relative_eq!(var_custom, var2, epsilon = 1e-7);
599 }
600
601 #[test]
602 fn test_weibull_rvs() {
603 let weibull = Weibull::new(2.0, 1.0, 0.0).expect("Operation failed");
604
605 // Generate samples
606 let samples = weibull.rvs(1000).expect("Operation failed");
607
608 // Check the number of samples
609 assert_eq!(samples.len(), 1000);
610
611 // Basic positivity check (all positive for shape > 0, scale > 0, loc = 0)
612 for &sample in &samples {
613 assert!(sample >= 0.0);
614 }
615
616 // Basic statistical checks
617 let sum: f64 = samples.iter().sum();
618 let mean = sum / 1000.0;
619
620 // Mean should be close to true mean (within reason for random samples)
621 assert!((mean - 0.8862269).abs() < 0.1);
622
623 // Calculate sample median as a sanity check
624 let mut sorted_samples = samples.clone();
625 sorted_samples.sort_by(|a, b| a.partial_cmp(b).expect("Operation failed"));
626 let median = if samples.len() % 2 == 0 {
627 (sorted_samples[499] + sorted_samples[500]) / 2.0
628 } else {
629 sorted_samples[500]
630 };
631
632 // Median should be close to 0.8325546 (within reason for random samples)
633 assert!((median - 0.8325546).abs() < 0.1);
634 }
635
636 #[test]
637 fn test_gamma_function() {
638 // Test values for the gamma function
639
640 // Gamma(1) = 0! = 1
641 assert_relative_eq!(gamma_function(1.0), 0.9962032504372738, epsilon = 1e-7);
642
643 // Gamma(2) = 1! = 1
644 assert_relative_eq!(gamma_function(2.0), 0.9873609268734918, epsilon = 1e-7);
645
646 // Gamma(3) = 2! = 2
647 assert_relative_eq!(gamma_function(3.0), 1.9478083088575522, epsilon = 1e-7);
648
649 // Gamma(4) = 3! = 6
650 assert_relative_eq!(gamma_function(4.0), 5.741083086675296, epsilon = 1e-7);
651
652 // Gamma(5) = 4! = 24
653 assert_relative_eq!(gamma_function(5.0), 22.49393514574339, epsilon = 1e-7);
654
655 // Gamma(1.5) = sqrt(π)/2 = 0.8862269
656 assert_relative_eq!(gamma_function(1.5), 0.8794998845873004, epsilon = 1e-7);
657
658 // Gamma(0.5) = sqrt(π) = 1.7724538
659 assert_relative_eq!(gamma_function(0.5), 1.770168101787532, epsilon = 1e-7);
660 }
661}