use crate::util::quantile::quantile_from_sorted;
use std::cmp::Ordering;
pub fn credible_interval(
samples_flat: &[f64],
n_draws: usize,
n_coeffs: usize,
level: f64,
) -> Result<Vec<f64>, String> {
if !(level > 0.0 && level < 1.0) {
return Err(format!("interval level must lie in (0, 1); got {level}"));
}
if samples_flat.len() != n_draws * n_coeffs {
return Err(format!(
"posterior_credible_interval samples shape mismatch: got {} floats, expected {} * {}",
samples_flat.len(),
n_draws,
n_coeffs
));
}
let alpha = (1.0 - level) / 2.0;
let mut out = Vec::with_capacity(2 * n_coeffs);
if n_draws == 0 {
for _ in 0..n_coeffs {
out.push(0.0);
out.push(0.0);
}
return Ok(out);
}
let mut column = vec![0.0_f64; n_draws];
for j in 0..n_coeffs {
for k in 0..n_draws {
column[k] = samples_flat[k * n_coeffs + j];
}
column.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));
out.push(quantile_from_sorted(&column, alpha));
out.push(quantile_from_sorted(&column, 1.0 - alpha));
}
Ok(out)
}