fdars_core/scalar_on_function/
logistic.rs

1use super::*;
2
3// ---------------------------------------------------------------------------
4// Functional logistic regression
5// ---------------------------------------------------------------------------
6
7/// One IRLS step: compute working response and solve weighted least squares.
8/// Returns updated beta or None if system is singular.
9fn irls_step(design: &FdMatrix, y: &[f64], beta: &[f64]) -> Option<Vec<f64>> {
10    let (n, p) = design.shape();
11
12    // Linear predictor η = Xβ, probabilities μ = sigmoid(η)
13    let eta: Vec<f64> = (0..n)
14        .map(|i| (0..p).map(|j| design[(i, j)] * beta[j]).sum())
15        .collect();
16    let mu: Vec<f64> = eta.iter().map(|&e| sigmoid(e)).collect();
17    let w: Vec<f64> = mu.iter().map(|&p| (p * (1.0 - p)).max(1e-10)).collect();
18    let z_work: Vec<f64> = (0..n).map(|i| eta[i] + (y[i] - mu[i]) / w[i]).collect();
19
20    // Weighted normal equations: (X'WX) β = X'Wz
21    let mut xtwx = vec![0.0; p * p];
22    for k in 0..p {
23        for j in k..p {
24            let mut s = 0.0;
25            for i in 0..n {
26                s += design[(i, k)] * w[i] * design[(i, j)];
27            }
28            xtwx[k * p + j] = s;
29            xtwx[j * p + k] = s;
30        }
31    }
32
33    let mut xtwz = vec![0.0; p];
34    for k in 0..p {
35        let mut s = 0.0;
36        for i in 0..n {
37            s += design[(i, k)] * w[i] * z_work[i];
38        }
39        xtwz[k] = s;
40    }
41
42    cholesky_solve(&xtwx, &xtwz, p).ok()
43}
44
45/// Compute log-likelihood of logistic model.
46fn logistic_log_likelihood(probabilities: &[f64], y: &[f64]) -> f64 {
47    probabilities
48        .iter()
49        .zip(y)
50        .map(|(&p, &yi)| {
51            let p = p.clamp(1e-15, 1.0 - 1e-15);
52            yi * p.ln() + (1.0 - yi) * (1.0 - p).ln()
53        })
54        .sum()
55}
56
57/// Run IRLS iteration loop and return (beta, iterations).
58fn irls_loop(design: &FdMatrix, y: &[f64], max_iter: usize, tol: f64) -> (Vec<f64>, usize) {
59    let p_total = design.ncols();
60    let mut beta = vec![0.0; p_total];
61    let mut iterations = 0;
62    for iter in 0..max_iter {
63        iterations = iter + 1;
64        let beta_new = match irls_step(design, y, &beta) {
65            Some(b) => b,
66            None => break,
67        };
68        let change: f64 = beta_new
69            .iter()
70            .zip(&beta)
71            .map(|(a, b)| (a - b).abs())
72            .fold(0.0_f64, f64::max);
73        beta = beta_new;
74        if change < tol {
75            break;
76        }
77    }
78    (beta, iterations)
79}
80
81/// Build logistic result from converged beta.
82fn build_logistic_result(
83    design: &FdMatrix,
84    beta: Vec<f64>,
85    y: &[f64],
86    fpca: FpcaResult,
87    ncomp: usize,
88    m: usize,
89    iterations: usize,
90) -> FunctionalLogisticResult {
91    let (n, p) = design.shape();
92    let eta = compute_fitted(design, &beta);
93    let probabilities: Vec<f64> = eta.iter().map(|&e| sigmoid(e)).collect();
94    let predicted_classes: Vec<usize> = probabilities
95        .iter()
96        .map(|&p| if p >= 0.5 { 1 } else { 0 })
97        .collect();
98    let correct: usize = predicted_classes
99        .iter()
100        .zip(y)
101        .filter(|(&pred, &actual)| pred as f64 == actual)
102        .count();
103    let beta_t = recover_beta_t(&beta[1..1 + ncomp], &fpca.rotation, m);
104    let gamma: Vec<f64> = beta[1 + ncomp..].to_vec();
105
106    // Compute coefficient SEs from Fisher information matrix (X'WX)^{-1}
107    let w: Vec<f64> = probabilities
108        .iter()
109        .map(|&mu| (mu * (1.0 - mu)).max(1e-10))
110        .collect();
111    let mut xtwx = vec![0.0; p * p];
112    for k in 0..p {
113        for j in k..p {
114            let mut s = 0.0;
115            for i in 0..n {
116                s += design[(i, k)] * w[i] * design[(i, j)];
117            }
118            xtwx[k * p + j] = s;
119            xtwx[j * p + k] = s;
120        }
121    }
122    let std_errors = cholesky_factor(&xtwx, p).map_or_else(
123        |_| vec![f64::NAN; p],
124        |l| compute_ols_std_errors(&l, p, 1.0),
125    );
126    let beta_se = compute_beta_se(&std_errors[1..1 + ncomp], &fpca.rotation, m);
127
128    let ll = logistic_log_likelihood(&probabilities, y);
129    let deviance = -2.0 * ll;
130    let nf = n as f64;
131    let pf = p as f64;
132    let aic = deviance + 2.0 * pf;
133    let bic = deviance + nf.ln() * pf;
134
135    FunctionalLogisticResult {
136        intercept: beta[0],
137        beta_t,
138        beta_se,
139        gamma,
140        accuracy: correct as f64 / nf,
141        log_likelihood: ll,
142        probabilities,
143        predicted_classes,
144        ncomp,
145        std_errors,
146        coefficients: beta,
147        iterations,
148        fpca,
149        aic,
150        bic,
151    }
152}
153
154/// Functional logistic regression for binary outcomes.
155///
156/// Fits: `log(P(Y=1)/P(Y=0)) = α + ∫β(t)X(t)dt + γᵀz`
157/// using IRLS (iteratively reweighted least squares) on FPC scores.
158///
159/// # Arguments
160/// * `data` - Functional predictor matrix (n × m)
161/// * `y` - Binary response vector (0.0 or 1.0, length n)
162/// * `scalar_covariates` - Optional scalar covariates (n × p)
163/// * `ncomp` - Number of FPC components
164/// * `max_iter` - Maximum IRLS iterations (default: 25)
165/// * `tol` - Convergence tolerance (default: 1e-6)
166///
167/// # Errors
168///
169/// Returns [`FdarError::InvalidDimension`] if `data` has fewer than 3 rows,
170/// zero columns, or `y.len() != n`.
171/// Returns [`FdarError::InvalidParameter`] if any element of `y` is not `0.0`
172/// or `1.0`.
173/// Returns [`FdarError::ComputationFailed`] if the SVD inside FPCA fails to
174/// extract components.
175#[must_use = "expensive computation whose result should not be discarded"]
176pub fn functional_logistic(
177    data: &FdMatrix,
178    y: &[f64],
179    scalar_covariates: Option<&FdMatrix>,
180    ncomp: usize,
181    max_iter: usize,
182    tol: f64,
183) -> Result<FunctionalLogisticResult, FdarError> {
184    let (n, m) = data.shape();
185    if n < 3 {
186        return Err(FdarError::InvalidDimension {
187            parameter: "data",
188            expected: "at least 3 rows".to_string(),
189            actual: format!("{n}"),
190        });
191    }
192    if m == 0 {
193        return Err(FdarError::InvalidDimension {
194            parameter: "data",
195            expected: "at least 1 column".to_string(),
196            actual: "0".to_string(),
197        });
198    }
199    if y.len() != n {
200        return Err(FdarError::InvalidDimension {
201            parameter: "y",
202            expected: format!("{n}"),
203            actual: format!("{}", y.len()),
204        });
205    }
206    if y.iter().any(|&yi| yi != 0.0 && yi != 1.0) {
207        return Err(FdarError::InvalidParameter {
208            parameter: "y",
209            message: "all values must be 0.0 or 1.0 for binary classification".to_string(),
210        });
211    }
212
213    let ncomp = ncomp.min(n - 1).min(m);
214    let fpca = fdata_to_pc_1d(data, ncomp)?;
215    let design = build_design_matrix(&fpca.scores, ncomp, scalar_covariates, n);
216
217    let max_iter = if max_iter == 0 { 25 } else { max_iter };
218    let tol = if tol <= 0.0 { 1e-6 } else { tol };
219
220    let (beta, iterations) = irls_loop(&design, y, max_iter, tol);
221    Ok(build_logistic_result(
222        &design, beta, y, fpca, ncomp, m, iterations,
223    ))
224}
225
226/// Predict probabilities P(Y=1) for new data using a fitted functional logistic model.
227///
228/// Projects new curves through the stored FPCA, computes linear predictor,
229/// and applies sigmoid.
230///
231/// # Arguments
232/// * `fit` — A fitted [`FunctionalLogisticResult`]
233/// * `new_data` — New functional predictor matrix (n_new × m)
234/// * `new_scalar` — Optional new scalar covariates (n_new × p)
235pub fn predict_functional_logistic(
236    fit: &FunctionalLogisticResult,
237    new_data: &FdMatrix,
238    new_scalar: Option<&FdMatrix>,
239) -> Vec<f64> {
240    let (n_new, m) = new_data.shape();
241    let ncomp = fit.ncomp;
242    let p_scalar = fit.gamma.len();
243
244    (0..n_new)
245        .map(|i| {
246            let mut eta = fit.coefficients[0]; // intercept
247            for k in 0..ncomp {
248                let mut s = 0.0;
249                for j in 0..m {
250                    s += (new_data[(i, j)] - fit.fpca.mean[j]) * fit.fpca.rotation[(j, k)];
251                }
252                eta += fit.coefficients[1 + k] * s;
253            }
254            if let Some(sc) = new_scalar {
255                for j in 0..p_scalar {
256                    eta += fit.gamma[j] * sc[(i, j)];
257                }
258            }
259            sigmoid(eta)
260        })
261        .collect()
262}
fdars_core/scalar_on_function/logistic.rs

fdars_core/scalar_on_function/
logistic.rs
