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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
use crate::{nctcdf, nctpdf, norminv, tinv};
/// Inverse of the noncentral T cumulative distribution function (CDF).
///
/// Given a probability `p`, degrees of freedom `nu`, and noncentrality parameter `delta`,
/// this function returns the value `x` such that the probability of a noncentral T
/// random variable being less than or equal to `x` is `p`.
///
/// # Mathematical Definition
/// The function finds $x$ such that:
/// <math display="block">
/// <msub><mi>F</mi><mrow><mi>NC</mi><mi>T</mi></mrow></msub>
/// <mo>(</mo><mi>x</mi><mo>;</mo><mi>ν</mi><mo>,</mo><mi>δ</mi><mo>)</mo>
/// <mo>=</mo>
/// <mi>p</mi>
/// </math>
/// where <math><msub><mi>F</mi><mrow><mi>NC</mi><mi>T</mi></mrow></msub></math> is the noncentral T CDF.
///
/// # Implementation
/// The inversion is performed using Newton's method with the standard normal
/// quantile as an initial guess. It incorporates step-size limiting and
/// back-tracking to ensure convergence.
///
/// # Domain
/// - <math><mn>0</mn><mo>≤</mo><mi>p</mi><mo>≤</mo><mn>1</mn></math>
/// - <math><mi>ν</mi><mo>></mo><mn>0</mn></math>
/// - <math><mi>δ</mi></math> must be finite.
/// - Returns `NaN` if `p` is out of range, `nu <= 0`, or `delta` is non-finite.
///
/// # Examples
/// ```
/// use abax::nctinv;
///
/// // For delta = 0, it reduces to the central Student's T distribution inverse.
/// let x = nctinv(0.5, 5.0, 0.0);
/// assert!((x - 0.0).abs() < 1e-12);
///
/// // Noncentral case
/// let x_nc = nctinv(0.5, 10.0, 2.0);
/// assert!((x_nc - 2.0536911511184894041).abs() < 1e-10);
/// ```
#[allow(non_snake_case)]
pub fn nctinv(p: f64, nu: f64, delta: f64) -> f64 {
if delta == 0.0 {
let x = tinv(p, nu);
return x;
}
let ok_params = (nu > 0.0 && !delta.is_nan()) && !delta.is_infinite();
if !ok_params {
return f64::NAN;
}
match p {
0.0 => return -f64::INFINITY,
1.0 => return f64::INFINITY,
_ if p < 0.0 || p > 1.0 => return f64::NAN,
_ => (),
}
let crit = f64::sqrt(f64::EPSILON);
let pk = p;
let vk = nu;
let dk = delta;
// Newton's Method
// Permit no more than count_limit iterations.
const COUNT_LIMIT: usize = 100;
let mut count = 0;
// Use norminv as a starting guess for x.
let mut xk = norminv(pk, dk, 1.0);
let mut h = 1.0;
// Break out of the iteration loop for the following:
// 1) The last update is very small (compared to x or in abs.value).
// 2) There are more than 100 iterations.
let mut F = nctcdf(xk, vk, dk, false);
while f64::abs(h) > crit * f64::abs(xk) && f64::abs(h) > crit && count < COUNT_LIMIT {
count += 1;
let f = nctpdf(xk, vk, dk);
h = (F - pk) / f;
// If h is inf or NaN, step closer to the original instead
let blowup = h.is_infinite() || h.is_nan();
if blowup {
h = xk / 10.0;
}
// Avoid stepping too far
let mut xnew = f64::max(-5.0 * f64::abs(xk), f64::min(5.0 * f64::abs(xk), xk - h));
// Back off if the step gives a worse result
let mut Fnew = nctcdf(xnew, vk, dk, false);
loop {
let worse = (f64::abs(Fnew - pk) > f64::abs(F - pk) * (1.0 + crit))
&& (f64::abs(xk - xnew) > crit * f64::abs(xk));
if !worse {
break;
}
xnew = 0.5 * (xnew + xk);
Fnew = nctcdf(xnew, vk, dk, false);
}
xk = xnew;
F = Fnew;
}
// Return the converged value(s).
if cfg!(debug_assertions) {
assert!(count < COUNT_LIMIT);
}
return xk;
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_nctinv_central() {
let nu = 10.0;
// Median should be 0.0
assert!((nctinv(0.5, nu, 0.0) - 0.0).abs() < 1e-14);
// Compare with tinv for non-median
assert!((nctinv(0.975, nu, 0.0) - tinv(0.975, nu)).abs() < 1e-14);
}
#[test]
fn test_nctinv_known_values() {
let tol = 1e-10;
assert!((nctinv(0.5, 5.0, 1.0) - 1.052851040947396e+00).abs() < tol);
assert!((nctinv(0.9, 10.0, 2.0) - 3.746633570263350e+00).abs() < tol);
}
#[test]
fn test_nctinv_symmetry() {
let p = 0.3;
let nu = 8.0;
let delta = 1.5;
// Inverse identity: nctinv(p, v, d) = -nctinv(1-p, v, -d)
let x1 = nctinv(p, nu, delta);
let x2 = nctinv(1.0 - p, nu, -delta);
assert!((x1 + x2).abs() < 1e-12);
}
#[test]
fn test_nctinv_boundaries() {
assert_eq!(nctinv(0.0, 5.0, 1.0), f64::NEG_INFINITY);
assert_eq!(nctinv(1.0, 5.0, 1.0), f64::INFINITY);
}
#[test]
fn test_nctinv_invalid() {
assert!(nctinv(-0.1, 5.0, 1.0).is_nan());
assert!(nctinv(1.1, 5.0, 1.0).is_nan());
assert!(nctinv(0.5, 0.0, 1.0).is_nan());
assert!(nctinv(0.5, -1.0, 1.0).is_nan());
assert!(nctinv(0.5, 5.0, f64::INFINITY).is_nan());
assert!(nctinv(0.5, 5.0, f64::NAN).is_nan());
}
}