copula_core/
model_selection.rs

1//! Model selection utilities such as cross-validation.
2use nalgebra::DMatrix;
3use rand::seq::SliceRandom;
4use rand::Rng;
5
6use crate::traits::FittableCopula;
7use crate::utils::validate_pseudo_observations;
8use crate::{CopulaError, Result};
9
10fn select_rows(matrix: &DMatrix<f64>, idx: &[usize]) -> DMatrix<f64> {
11    let ncols = matrix.ncols();
12    let mut out = DMatrix::<f64>::zeros(idx.len(), ncols);
13    for (i, &row_idx) in idx.iter().enumerate() {
14        for j in 0..ncols {
15            out[(i, j)] = matrix[(row_idx, j)];
16        }
17    }
18    out
19}
20
21/// Perform k-fold cross-validation for a copula model.
22///
23/// The function splits the pseudo-observations into `k` folds,
24/// fitting the model on k-1 folds and evaluating the log-likelihood on
25/// the remaining fold. The returned value is the average log-likelihood
26/// across all folds.
27pub fn k_fold_cv<C, R>(template: C, pseudo_obs: &DMatrix<f64>, k: usize, rng: &mut R) -> Result<f64>
28where
29    C: FittableCopula + Clone,
30    R: Rng + ?Sized,
31{
32    validate_pseudo_observations(pseudo_obs)?;
33    let n = pseudo_obs.nrows();
34    if k < 2 || k > n {
35        return Err(CopulaError::invalid_parameter("k must be between 2 and n"));
36    }
37
38    let mut indices: Vec<usize> = (0..n).collect();
39    indices.shuffle(rng);
40    let fold_size = n.div_ceil(k);
41    let mut total_ll = 0.0;
42    let mut folds_used = 0;
43
44    for fold in 0..k {
45        let start = fold * fold_size;
46        if start >= n {
47            break;
48        }
49        let end = ((fold + 1) * fold_size).min(n);
50        let test_idx = &indices[start..end];
51        if test_idx.is_empty() {
52            continue;
53        }
54        let train_idx: Vec<usize> = indices[..start]
55            .iter()
56            .chain(&indices[end..])
57            .copied()
58            .collect();
59        let train = select_rows(pseudo_obs, &train_idx);
60        let test = select_rows(pseudo_obs, test_idx);
61
62        let mut model = template.clone();
63        model.fit(&train)?;
64        let ll = model.log_likelihood(&test)?;
65        total_ll += ll;
66        folds_used += 1;
67    }
68
69    Ok(total_ll / folds_used as f64)
70}
copula_core/model_selection.rs

copula_core/
model_selection.rs
