fdars_core/scalar_on_function/
logistic.rs

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

fdars_core/scalar_on_function/
logistic.rs
