fdars_core/irreg_fdata/
kernels.rs

1//! Kernel functions and mean estimation for irregular functional data.
2
3use crate::slice_maybe_parallel;
4#[cfg(feature = "parallel")]
5use rayon::iter::ParallelIterator;
6
7use super::IrregFdata;
8
9/// Kernel function type for smoothing operations.
10#[derive(Clone, Copy, Debug, PartialEq, Eq)]
11pub enum KernelType {
12    /// Epanechnikov kernel: K(u) = 0.75(1 - u^2) for |u| <= 1
13    Epanechnikov,
14    /// Gaussian kernel: K(u) = exp(-u^2/2) / sqrt(2*pi)
15    Gaussian,
16}
17
18/// Epanechnikov kernel function.
19#[inline]
20pub(crate) fn kernel_epanechnikov(u: f64) -> f64 {
21    if u.abs() <= 1.0 {
22        0.75 * (1.0 - u * u)
23    } else {
24        0.0
25    }
26}
27
28/// Gaussian kernel function.
29#[inline]
30pub(crate) fn kernel_gaussian(u: f64) -> f64 {
31    (-0.5 * u * u).exp() / (2.0 * std::f64::consts::PI).sqrt()
32}
33
34impl KernelType {
35    #[inline]
36    pub(crate) fn as_fn(self) -> fn(f64) -> f64 {
37        match self {
38            KernelType::Epanechnikov => kernel_epanechnikov,
39            KernelType::Gaussian => kernel_gaussian,
40        }
41    }
42}
43
44/// Estimate mean function at specified target points using kernel smoothing.
45///
46/// Uses local weighted averaging (Nadaraya-Watson estimator) at each target point:
47/// mu_hat(t) = sum_{i,j} K_h(t - t_{ij}) x_{ij} / sum_{i,j} K_h(t - t_{ij})
48///
49/// # Arguments
50/// * `ifd` - Irregular functional data
51/// * `target_argvals` - Points at which to estimate the mean
52/// * `bandwidth` - Kernel bandwidth
53/// * `kernel_type` - Kernel function to use
54///
55/// # Returns
56/// Estimated mean function values at target points
57pub fn mean_irreg(
58    ifd: &IrregFdata,
59    target_argvals: &[f64],
60    bandwidth: f64,
61    kernel_type: KernelType,
62) -> Vec<f64> {
63    let n = ifd.n_obs();
64    let kernel = kernel_type.as_fn();
65
66    slice_maybe_parallel!(target_argvals)
67        .map(|&t| {
68            let mut sum_weights = 0.0;
69            let mut sum_values = 0.0;
70
71            for i in 0..n {
72                let (obs_t, obs_x) = ifd.get_obs(i);
73
74                for (&ti, &xi) in obs_t.iter().zip(obs_x.iter()) {
75                    let u = (ti - t) / bandwidth;
76                    let w = kernel(u);
77                    sum_weights += w;
78                    sum_values += w * xi;
79                }
80            }
81
82            if sum_weights > 0.0 {
83                sum_values / sum_weights
84            } else {
85                f64::NAN
86            }
87        })
88        .collect()
89}
fdars_core/irreg_fdata/kernels.rs

fdars_core/irreg_fdata/
kernels.rs
