1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
/// Weibull inverse cumulative distribution function (quantile function) for scalar values.
///
/// Evaluates the inverse CDF of the Weibull distribution with scale parameter `a`
/// and shape parameter `b` at the probability value `p`.
///
/// # Arguments
/// * `p` - The probability value at which to evaluate the inverse CDF (0 ≤ p ≤ 1).
/// * `a` - Scale parameter (a > 0).
/// * `b` - Shape parameter (b > 0).
///
/// # Returns
/// The quantile value `x` at the probability `p`. Returns `f64::NAN` if inputs or parameters are invalid.
///
/// # Mathematical Formula
/// For <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mn>0</mn><mo>≤</mo><mi>p</mi><mo>≤</mo><mn>1</mn></math>,
/// <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mi>a</mi><mo>></mo><mn>0</mn></math>,
/// <math xmlns="http://www.w3.org/1998/Math/MathML" display="inline"><mi>b</mi><mo>></mo><mn>0</mn></math>:
///
/// <math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
/// <mi>x</mi>
/// <mo>=</mo>
/// <msup>
/// <mrow>
/// <mi>a</mi>
/// <mo>⋅</mo>
/// <mo>[</mo>
/// <mo>−</mo>
/// <mi>ln</mi>
/// <mo stretchy="false">(</mo>
/// <mn>1</mn>
/// <mo>−</mo>
/// <mi>p</mi>
/// <mo stretchy="false">)</mo>
/// <mo>]</mo>
/// </mrow>
/// <mrow>
/// <mfrac>
/// <mn>1</mn>
/// <mi>b</mi>
/// </mfrac>
/// </mrow>
/// </msup>
/// </math>
pub fn wblinv(p: f64, a: f64, b: f64) -> f64 {
// 1. Return NaN for incoming NaNs or out of range parameters
if p.is_nan() || a.is_nan() || b.is_nan() || a <= 0.0 || b <= 0.0 {
return f64::NAN;
}
// 2. Out of bounds checking for probability domain
if p < 0.0 || p > 1.0 {
return f64::NAN;
}
// 3. Exact evaluation boundary cases
if p == 0.0 {
return 0.0;
}
if p == 1.0 {
return f64::INFINITY;
}
// 4. Transform probability to the intermediate quantile variable 'q'
// log1p(-p) evaluates ln(1 + (-p)) = ln(1 - p) precisely for values near 0.
let q = -(-p).ln_1p();
// 5. Calculate final scale translation
a * q.powf(1.0 / b)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_standard_weibull_inverse() {
// For a=1, b=1 (Exponential distribution), inverse CDF at p=0.5 is -ln(0.5)
let res = wblinv(0.5, 1.0, 1.0);
let expected = -0.5f64.ln();
assert!((res - expected).abs() < 1e-15);
}
#[test]
fn test_boundary_probabilities() {
assert_eq!(wblinv(0.0, 2.5, 1.5), 0.0);
assert_eq!(wblinv(1.0, 2.5, 1.5), f64::INFINITY);
}
#[test]
fn test_invalid_probabilities() {
assert!(wblinv(-0.01, 1.0, 1.0).is_nan());
assert!(wblinv(1.01, 1.0, 1.0).is_nan());
assert!(wblinv(f64::NAN, 1.0, 1.0).is_nan());
}
#[test]
fn test_invalid_parameters() {
assert!(wblinv(0.5, 0.0, 1.0).is_nan());
assert!(wblinv(0.5, 1.0, -2.0).is_nan());
assert!(wblinv(0.5, -1.0, 1.0).is_nan());
}
#[test]
fn test_round_trip() {
// Round trip verification from a custom parameter group
let p_initial = 0.75;
let a = 1200.0;
let b = 2.5;
// Manual verification logic mapping (wblcdf lower tail structure)
let x = wblinv(p_initial, a, b);
let z = (x / a).powf(b);
let p_recovered = -(-z).exp_m1();
assert!((p_initial - p_recovered).abs() < 1e-15);
}
}