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
//! Intercept of an additive Gamma GLM by Newton iteration.
//!
//! Pure-Rust port of limma's `fitGammaIntercept.R` ([`fit_gamma_intercept`]),
//! used by the gene-set machinery to estimate a common additive offset.
/// `fitGammaIntercept(y, offset, maxit)`. Estimates the intercept `x` solving
/// `sum(y / (offset + x)) = length(y)` by Newton's method. A constant `offset`
/// (range within `1e-14`, including a scalar passed as a length-1 slice)
/// short-circuits to `mean(y) - offset`. Panics on negative `y` or `offset`.
/// limma's default `maxit` is 1000.
pub fn fit_gamma_intercept(y: &[f64], offset: &[f64], maxit: usize) -> f64 {
assert!(y.iter().all(|&v| v >= 0.0), "negative y not permitted");
assert!(offset.iter().all(|&v| v >= 0.0), "offsets must be positive");
assert!(
offset.len() == y.len() || offset.len() == 1,
"offset must be a scalar or match the length of y"
);
let (omin, omax) = offset
.iter()
.fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), &v| {
(lo.min(v), hi.max(v))
});
if omin + 1e-14 > omax {
let mean = y.iter().sum::<f64>() / y.len() as f64;
return mean - omin;
}
// Non-constant offset implies offset.len() == y.len() here.
let n = y.len() as f64;
let mut x = 0.0_f64;
let mut iter = 0usize;
loop {
iter += 1;
let mut q = 0.0;
let mut dq = 0.0;
for (i, &yi) in y.iter().enumerate() {
let d = offset[i] + x;
q += yi / d;
dq += yi / (d * d);
}
let dif = (q - n) / dq;
x += dif;
if dif.abs() < 1e-8 || iter > maxit {
break;
}
}
x
}
#[cfg(test)]
mod tests {
use super::*;
fn close(a: f64, b: f64, tol: f64) -> bool {
(a - b).abs() <= tol + tol * b.abs()
}
#[test]
fn nonconstant_offset_matches_r() {
// Reference: fitGammaIntercept(1:5, c(0.5,1,1.5,2,2.5)) in limma 3.68.3.
let y = [1.0, 2.0, 3.0, 4.0, 5.0];
let offset = [0.5, 1.0, 1.5, 2.0, 2.5];
assert!(close(
fit_gamma_intercept(&y, &offset, 1000),
1.31379054922341,
1e-9
));
}
#[test]
fn constant_offset_is_mean_shift() {
// Reference: fitGammaIntercept(c(2,4,6), c(1,1,1)) == mean - 1 == 3.
assert!(close(
fit_gamma_intercept(&[2.0, 4.0, 6.0], &[1.0, 1.0, 1.0], 1000),
3.0,
1e-12
));
}
#[test]
fn scalar_zero_offset_is_mean() {
// Reference: fitGammaIntercept(c(2,4,6), 0) == mean == 4.
assert!(close(
fit_gamma_intercept(&[2.0, 4.0, 6.0], &[0.0], 1000),
4.0,
1e-12
));
}
}