1use std::collections::HashMap;
173use std::convert::Infallible;
174use std::sync::{Arc, Mutex, OnceLock};
175
176use crate::estimate::EstimationError;
177use crate::mixture_link::{
178 beta_logistic_inverse_link_jet, component_inverse_link_jet, sas_inverse_link_jet,
179};
180use gam_math::probability::{erfcx_nonnegative, normal_logcdf};
181use gam_math::special::stable_polynomial_times_exp_neg as cloglog_stable_poly_times_exp_neg;
182use gam_problem::types::{
183 GlmLikelihoodSpec, InverseLink, LinkComponent, LinkFunction, MixtureLinkState, ResponseFamily,
184 SasLinkState, StandardLink,
185};
186const N_POINTS: usize = 7;
188const SQRT_2: f64 = std::f64::consts::SQRT_2;
189const QUADRATURE_EXP_LOG_MAX: f64 = 700.0;
190
191#[inline]
194fn safe_exp(x: f64) -> f64 {
195 if x.is_nan() {
196 f64::NAN
197 } else {
198 x.min(QUADRATURE_EXP_LOG_MAX).exp()
199 }
200}
201
202#[inline]
203fn safe_expwith_saturation(x: f64) -> (f64, bool) {
204 (safe_exp(x), x > QUADRATURE_EXP_LOG_MAX)
205}
206
207#[derive(Clone, Copy, Debug, Default)]
208struct Complex {
209 re: f64,
210 im: f64,
211}
212
213pub struct QuadratureContext {
215 gh_cache: OnceLock<GaussHermiteRule>,
216 gh15_cache: OnceLock<GaussHermiteRuleDynamic>,
217 gh21_cache: OnceLock<GaussHermiteRuleDynamic>,
218 gh31_cache: OnceLock<GaussHermiteRuleDynamic>,
219 gh51_cache: OnceLock<GaussHermiteRuleDynamic>,
220 cc_cache: Mutex<HashMap<usize, Arc<ClenshawCurtisRule>>>,
224}
225
226#[derive(Clone, Copy, Debug, Eq, PartialEq)]
227pub enum IntegratedExpectationMode {
228 ExactClosedForm,
229 ExactSpecialFunction,
230 ControlledAsymptotic,
231 QuadratureFallback,
232}
233
234impl IntegratedExpectationMode {
235 #[inline]
239 pub const fn rank(self) -> u8 {
240 match self {
241 Self::ExactClosedForm => 0,
242 Self::ExactSpecialFunction => 1,
243 Self::ControlledAsymptotic => 2,
244 Self::QuadratureFallback => 3,
245 }
246 }
247}
248
249#[derive(Clone, Copy, Debug)]
250pub struct IntegratedMeanDerivative {
251 pub mean: f64,
252 pub dmean_dmu: f64,
253 pub mode: IntegratedExpectationMode,
254}
255
256#[derive(Clone, Copy, Debug)]
257pub struct IntegratedInverseLinkJet {
258 pub mean: f64,
259 pub d1: f64,
260 pub d2: f64,
261 pub d3: f64,
262 pub mode: IntegratedExpectationMode,
263}
264
265#[derive(Clone, Copy, Debug)]
266pub(crate) struct IntegratedInverseLinkJet5 {
267 pub mean: f64,
268 pub d1: f64,
269 pub d2: f64,
270 pub d3: f64,
271 pub d4: f64,
272 pub d5: f64,
273 pub mode: IntegratedExpectationMode,
274}
275
276#[inline]
277pub(crate) fn validate_latent_cloglog_inputs(eta: f64, sigma: f64) -> Result<(), EstimationError> {
278 if !eta.is_finite() || !sigma.is_finite() || sigma < 0.0 {
279 crate::bail_invalid_estim!(
280 "latent cloglog jet requires finite eta and sigma >= 0, got eta={eta}, sigma={sigma}"
281 );
282 }
283 Ok::<(), _>(())
284}
285
286#[derive(Clone, Copy, Debug)]
291pub struct IntegratedMomentsJet {
292 pub mean: f64,
293 pub variance: f64,
294 pub d1: f64,
295 pub d2: f64,
296 pub d3: f64,
297 pub mode: IntegratedExpectationMode,
298}
299
300const LOGIT_SIGMA_DEGENERATE: f64 = 1e-10;
301const LOGIT_SIGMA_TAYLOR_MAX: f64 = 2.5e-1;
302const LOGIT_TAIL_LOG_MAX: f64 = -18.0;
303const LOGIT_ERFCX_MU_MAX: f64 = 40.0;
304const LOGIT_ERFCX_SIGMA_MAX: f64 = 6.0;
305const LOGIT_JET_GHQ_SIGMA_MAX: f64 = 1.0;
329const CLOGLOG_SIGMA_DEGENERATE: f64 = 1e-10;
330const CLOGLOG_SIGMA_TAYLOR_MAX: f64 = 0.25;
331const CLOGLOG_JET_MOMENT_SIGMA_MAX: f64 = 1.0;
338const CLOGLOG_RARE_EVENT_LOG_MAX: f64 = -18.0;
339const CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN: f64 = 8.0;
340const CLOGLOG_POSITIVE_SATURATION_EDGE: f64 = 5.0;
341const CLOGLOG_POSITIVE_SATURATION_SIGMAS: f64 = 8.0;
342const CLOGLOG_GUMBEL_QUAD_ETA_LO: f64 = -40.0;
348const CLOGLOG_GUMBEL_QUAD_ETA_HI: f64 = 6.0;
349const CLOGLOG_GUMBEL_QUAD_MIN_NODES: usize = 97;
353const CLOGLOG_GUMBEL_QUAD_NODE_SCALE: f64 = 320.0;
357const CLOGLOG_GUMBEL_QUAD_MAX_NODES: usize = 513;
358const SERIES_CONSECUTIVE_SMALL_TERMS: usize = 6;
359const LOGIT_MAX_TERMS: usize = 160;
360const LOGIT_ERFCX_ACCURACY_TARGET: f64 = 1.0e-11;
376const CLOGLOG_MILES_ALPHA: f64 = 60.0;
377const CLOGLOG_MILES_MAX_TERMS: usize = 256;
378const CLOGLOG_MILES_PEAK_LOG_MAX: f64 = 0.0;
395const CLOGLOG_GAMMA_K_REF: f64 = 0.5;
396const CLOGLOG_GAMMA_T_MAX_REF: f64 = 24.0;
397const CLOGLOG_GAMMA_H_REF: f64 = 0.01;
398const CLOGLOG_CC_TOL: f64 = 1e-12;
402const CLOGLOG_CC_NODE_CAP: usize = 1025;
406const CLOGLOG_GAMMA_SAMPLE_COUNT: usize =
410 (CLOGLOG_GAMMA_T_MAX_REF / CLOGLOG_GAMMA_H_REF) as usize + 1;
411const CLOGLOG_CC_PREFER_THRESHOLD: usize = CLOGLOG_GAMMA_SAMPLE_COUNT / 3;
416const CLOGLOG_CC_MIN_N: usize = 17;
420
421impl QuadratureContext {
422 pub fn new() -> Self {
423 Self {
424 gh_cache: OnceLock::new(),
425 gh15_cache: OnceLock::new(),
426 gh21_cache: OnceLock::new(),
427 gh31_cache: OnceLock::new(),
428 gh51_cache: OnceLock::new(),
429 cc_cache: Mutex::new(HashMap::new()),
430 }
431 }
432
433 fn gauss_hermite(&self) -> &GaussHermiteRule {
434 self.gh_cache.get_or_init(compute_gauss_hermite)
435 }
436
437 fn gauss_hermite_n(&self, n: usize) -> &GaussHermiteRuleDynamic {
438 match n {
439 7 => self.gh15_cache.get_or_init(|| compute_gauss_hermite_n(15)),
442 15 => self.gh15_cache.get_or_init(|| compute_gauss_hermite_n(15)),
443 21 => self.gh21_cache.get_or_init(|| compute_gauss_hermite_n(21)),
444 31 => self.gh31_cache.get_or_init(|| compute_gauss_hermite_n(31)),
445 51 => self.gh51_cache.get_or_init(|| compute_gauss_hermite_n(51)),
446 _ => self.gh21_cache.get_or_init(|| compute_gauss_hermite_n(21)),
447 }
448 }
449
450 fn clenshaw_curtis_n(&self, n: usize) -> Arc<ClenshawCurtisRule> {
451 let mut cache = match self.cc_cache.lock() {
452 Ok(guard) => guard,
453 Err(poisoned) => poisoned.into_inner(),
454 };
455 cache
456 .entry(n)
457 .or_insert_with(|| Arc::new(compute_clenshaw_curtis_n(n)))
458 .clone()
459 }
460}
461
462impl Default for QuadratureContext {
463 fn default() -> Self {
464 Self::new()
465 }
466}
467
468struct GaussHermiteRule {
470 nodes: [f64; N_POINTS],
472 weights: [f64; N_POINTS],
474}
475
476pub(crate) struct GaussHermiteRuleDynamic {
477 pub(crate) nodes: Vec<f64>,
478 pub(crate) weights: Vec<f64>,
479}
480
481#[derive(Clone)]
482struct ClenshawCurtisRule {
483 nodes: Vec<f64>,
484 weights: Vec<f64>,
485}
486
487fn compute_clenshaw_curtis_n(n: usize) -> ClenshawCurtisRule {
488 assert!(
489 n >= 2,
490 "Clenshaw-Curtis rule requires at least two nodes: n={n}"
491 );
492 let m = n - 1;
507 let theta: Vec<f64> = (0..=m)
508 .map(|j| std::f64::consts::PI * (j as f64) / (m as f64))
509 .collect();
510 let nodes: Vec<f64> = theta.iter().map(|&th| th.cos()).collect();
511
512 if n == 2 {
513 return ClenshawCurtisRule {
514 nodes,
515 weights: vec![1.0, 1.0],
516 };
517 }
518
519 let mut weights = vec![0.0_f64; n];
520 let mut v = vec![1.0_f64; m - 1];
521
522 if m.is_multiple_of(2) {
523 let w0 = 1.0 / ((m * m - 1) as f64);
524 weights[0] = w0;
525 weights[m] = w0;
526 for k in 1..(m / 2) {
527 let denom = (4 * k * k - 1) as f64;
528 for j in 1..m {
529 v[j - 1] -= 2.0 * (2.0 * (k as f64) * theta[j]).cos() / denom;
530 }
531 }
532 for j in 1..m {
533 v[j - 1] -= ((m as f64) * theta[j]).cos() / ((m * m - 1) as f64);
534 }
535 } else {
536 let w0 = 1.0 / ((m * m) as f64);
537 weights[0] = w0;
538 weights[m] = w0;
539 for k in 1..=((m - 1) / 2) {
540 let denom = (4 * k * k - 1) as f64;
541 for j in 1..m {
542 v[j - 1] -= 2.0 * (2.0 * (k as f64) * theta[j]).cos() / denom;
543 }
544 }
545 }
546
547 for j in 1..m {
548 weights[j] = 2.0 * v[j - 1] / (m as f64);
549 }
550
551 for j in 0..=(m / 2) {
555 let jj = m - j;
556 let avg = 0.5 * (weights[j] + weights[jj]);
557 weights[j] = avg;
558 weights[jj] = avg;
559 }
560 let weight_sum: f64 = weights.iter().sum();
561 if weight_sum.is_finite() && weight_sum != 0.0 {
562 let scale = 2.0 / weight_sum;
563 for w in &mut weights {
564 *w *= scale;
565 }
566 }
567
568 ClenshawCurtisRule { nodes, weights }
569}
570
571fn cloglog_cc_required_nodes(mu: f64, sigma: f64, tol: f64) -> Result<usize, EstimationError> {
572 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0 && tol.is_finite() && tol > 0.0) {
573 crate::bail_invalid_estim!(
574 "CC cloglog backend requires finite mu, positive sigma, and positive tolerance"
575 .to_string(),
576 );
577 }
578
579 let p_tail = (tol / 8.0).clamp(1e-300, 0.25);
584 let a = gam_math::probability::standard_normal_quantile(p_tail)
585 .map(|z| -z)
586 .unwrap_or(8.0)
587 .max(1.0);
588
589 let ay = a * sigma;
590 let y = if ay > 0.0 {
591 1.0_f64.min(std::f64::consts::PI / (4.0 * ay))
592 } else {
593 1.0
594 };
595 let rho = y + (1.0 + y * y).sqrt();
596 let m_s = (0.5 * (a * y) * (a * y)).exp() / (2.0 * std::f64::consts::PI).sqrt();
597 let eps_quad = (tol / 4.0).max(1e-300);
598 let numer = ((8.0 * a * m_s) / ((rho - 1.0).max(1e-12) * eps_quad)).max(1.0);
599 let denom = rho.ln();
600 if !denom.is_finite() || denom <= 0.0 {
601 crate::bail_invalid_estim!("CC cloglog backend ellipse bound became degenerate");
602 }
603
604 let mut n = (1.0 + numer.ln() / denom).ceil() as usize;
605 n = n.max(CLOGLOG_CC_MIN_N);
606 if n.is_multiple_of(2) {
607 n += 1;
608 }
609 Ok(n)
610}
611
612#[inline]
613fn cloglog_should_prefer_cc(mu: f64, sigma: f64, tol: f64) -> bool {
614 match cloglog_cc_required_nodes(mu, sigma, tol) {
620 Ok(n) => n <= CLOGLOG_CC_PREFER_THRESHOLD,
621 Err(_) => false,
622 }
623}
624
625fn compute_gauss_hermite() -> GaussHermiteRule {
638 let mut diag = [0.0f64; N_POINTS]; let mut off_diag = [0.0f64; N_POINTS - 1];
644
645 for i in 0..(N_POINTS - 1) {
646 off_diag[i] = (((i + 1) as f64) / 2.0).sqrt();
648 }
649
650 let (eigenvalues, eigenvectors) = symmetric_tridiagonal_eigen(&mut diag, &mut off_diag);
653
654 let nodes = eigenvalues;
656 let mut weights = [0.0f64; N_POINTS];
657
658 let mu0 = std::f64::consts::PI.sqrt();
664 for i in 0..N_POINTS {
665 let v0 = eigenvectors[i][0];
666 weights[i] = mu0 * v0 * v0;
667 }
668
669 let mut indices: [usize; N_POINTS] = [0, 1, 2, 3, 4, 5, 6];
671 indices.sort_by(|&a, &b| nodes[a].total_cmp(&nodes[b]));
672
673 let sorted_nodes: [f64; N_POINTS] = std::array::from_fn(|i| nodes[indices[i]]);
674 let sortedweights: [f64; N_POINTS] = std::array::from_fn(|i| weights[indices[i]]);
675
676 GaussHermiteRule {
677 nodes: sorted_nodes,
678 weights: sortedweights,
679 }
680}
681
682pub(crate) fn compute_gauss_hermite_n(n: usize) -> GaussHermiteRuleDynamic {
683 let mut diag = vec![0.0f64; n];
684 let mut off_diag = vec![0.0f64; n.saturating_sub(1)];
685 for (i, od) in off_diag.iter_mut().enumerate() {
686 *od = (((i + 1) as f64) / 2.0).sqrt();
687 }
688 let (nodes, eigenvectors) = symmetric_tridiagonal_eigen_dynamic(&mut diag, &mut off_diag);
689 let mu0 = std::f64::consts::PI.sqrt();
690 let mut pairs = (0..n)
691 .map(|i| {
692 let v0 = eigenvectors[i][0];
693 (nodes[i], mu0 * v0 * v0)
694 })
695 .collect::<Vec<_>>();
696 pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
697 GaussHermiteRuleDynamic {
698 nodes: pairs.iter().map(|p| p.0).collect(),
699 weights: pairs.iter().map(|p| p.1).collect(),
700 }
701}
702
703fn symmetric_tridiagonal_eigen(
707 diag: &mut [f64; N_POINTS],
708 off_diag: &mut [f64; N_POINTS - 1],
709) -> ([f64; N_POINTS], [[f64; N_POINTS]; N_POINTS]) {
710 let mut diag_vec = diag.to_vec();
711 let mut off_diag_vec = off_diag.to_vec();
712 let (eigenvalues, eigenvectors) =
713 symmetric_tridiagonal_eigen_dynamic(&mut diag_vec, &mut off_diag_vec);
714
715 let mut values = [0.0; N_POINTS];
716 let mut vectors = [[0.0; N_POINTS]; N_POINTS];
717 values.copy_from_slice(&eigenvalues);
718 for i in 0..N_POINTS {
719 vectors[i].copy_from_slice(&eigenvectors[i]);
720 }
721 diag.copy_from_slice(&values);
722 off_diag.copy_from_slice(&off_diag_vec);
723 (values, vectors)
724}
725
726fn symmetric_tridiagonal_eigen_dynamic(
727 diag: &mut [f64],
728 off_diag: &mut [f64],
729) -> (Vec<f64>, Vec<Vec<f64>>) {
730 let dim = diag.len();
731 let mut z = vec![vec![0.0_f64; dim]; dim];
732 for (i, row) in z.iter_mut().enumerate().take(dim) {
733 row[i] = 1.0;
734 }
735 const DEFLATION_TOL: f64 = 1e-15;
741 const MAX_QL_SWEEPS: usize = 200;
742 let eps = DEFLATION_TOL;
743 let max_iter = MAX_QL_SWEEPS;
744 let mut t_norm = 0.0_f64;
750 for i in 0..dim {
751 let left = if i > 0 { off_diag[i - 1].abs() } else { 0.0 };
752 let right = if i + 1 < dim { off_diag[i].abs() } else { 0.0 };
753 let row_sum = diag[i].abs() + left + right;
754 if row_sum > t_norm {
755 t_norm = row_sum;
756 }
757 }
758 let mut n = dim;
759 while n > 1 {
760 let mut converged = false;
761 for _ in 0..max_iter {
762 let mut m = n - 1;
763 while m > 0 {
764 let row_scale = (diag[m - 1].abs() + diag[m].abs()).max(t_norm);
765 if off_diag[m - 1].abs() <= eps * row_scale {
766 off_diag[m - 1] = 0.0;
767 break;
768 }
769 m -= 1;
770 }
771 if m == n - 1 {
772 n -= 1;
773 converged = true;
774 break;
775 }
776 let shift = wilkinson_shift(diag[n - 2], diag[n - 1], off_diag[n - 2]);
777 let mut x = diag[m] - shift;
778 let mut y = off_diag[m];
779 for k in m..(n - 1) {
780 let (c, s) = if y.abs() > eps {
781 let r = x.hypot(y);
782 if r > 0.0 && r.is_finite() {
783 (x / r, -y / r)
784 } else {
785 (1.0, 0.0)
786 }
787 } else {
788 (1.0, 0.0)
789 };
790 if k > m {
791 off_diag[k - 1] = x.hypot(y);
792 }
793 let d1 = diag[k];
794 let d2 = diag[k + 1];
795 let e_k = off_diag[k];
796 diag[k] = c * c * d1 + s * s * d2 - 2.0 * c * s * e_k;
797 diag[k + 1] = s * s * d1 + c * c * d2 + 2.0 * c * s * e_k;
798 off_diag[k] = c * s * (d1 - d2) + (c * c - s * s) * e_k;
799 if k < n - 2 {
800 x = off_diag[k];
801 y = -s * off_diag[k + 1];
802 off_diag[k + 1] *= c;
803 }
804 for i in 0..dim {
805 let t = z[k][i];
806 z[k][i] = c * t - s * z[k + 1][i];
807 z[k + 1][i] = s * t + c * z[k + 1][i];
808 }
809 }
810 }
811 if !converged {
812 off_diag[n - 2] = 0.0;
813 n -= 1;
814 }
815 }
816 (diag.to_vec(), z)
817}
818
819#[inline]
820fn wilkinson_shift(a: f64, c: f64, b: f64) -> f64 {
821 let d = (a - c) * 0.5;
822 let t = d.hypot(b);
823 let sgn = if d >= 0.0 { 1.0 } else { -1.0 }; let denom = d + sgn * t;
825
826 if denom.abs() > f64::EPSILON * t.max(1.0) {
827 c - (b * b) / denom
828 } else {
829 c - t
831 }
832}
833
834#[inline]
845pub fn logit_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
846 match logit_posterior_meanwith_deriv_controlled(eta, se_eta) {
847 Ok(out) => out.mean,
848 Err(_) => integrate_normal_ghq_adaptive(ctx, eta, se_eta, sigmoid),
849 }
850}
851
852#[inline]
860pub fn logit_posterior_meanwith_deriv(
861 eta: f64,
862 se_eta: f64,
863) -> Result<(f64, f64), EstimationError> {
864 let out = logit_posterior_meanwith_deriv_controlled(eta, se_eta)?;
875 Ok((out.mean, out.dmean_dmu))
876}
877
878#[inline]
879pub fn probit_posterior_meanwith_deriv_exact(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
880 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 1e-12 {
908 let mean = gam_math::probability::normal_cdf(mu);
909 let dmean_dmu = gam_math::probability::normal_pdf(mu);
910 return IntegratedMeanDerivative {
911 mean,
912 dmean_dmu,
913 mode: IntegratedExpectationMode::ExactClosedForm,
914 };
915 }
916 let denom = (1.0 + sigma * sigma).sqrt();
917 let z = mu / denom;
918 IntegratedMeanDerivative {
919 mean: gam_math::probability::normal_cdf(z),
920 dmean_dmu: gam_math::probability::normal_pdf(z) / denom,
921 mode: IntegratedExpectationMode::ExactClosedForm,
922 }
923}
924
925#[inline]
926fn logistic_normal_exact_eligible(mu: f64, sigma: f64) -> bool {
927 mu.is_finite()
928 && sigma.is_finite()
929 && mu.abs() <= LOGIT_ERFCX_MU_MAX
930 && (LOGIT_SIGMA_TAYLOR_MAX..=LOGIT_ERFCX_SIGMA_MAX).contains(&sigma)
931}
932
933#[inline]
977fn logistic_normal_series_cutoff(mu: f64, sigma: f64, target_accuracy: f64) -> Option<usize> {
978 assert!(sigma > 0.0);
979 assert!(target_accuracy > 0.0);
980 let m = mu.abs();
981 let s = sigma;
982 let gauss = (-(m * m) / (2.0 * s * s)).exp();
983 let coeff_mean = m * (2.0_f64 / std::f64::consts::PI).sqrt() * gauss / (s * s * s);
984 let coeff_deriv =
985 2.0 * gauss * (m * m - s * s).abs() / ((2.0 * std::f64::consts::PI).sqrt() * s.powi(5));
986 let asymptotic_index = |coeff: f64| -> f64 {
990 if !coeff.is_finite() || coeff <= target_accuracy {
991 0.0
992 } else {
993 (coeff / target_accuracy).sqrt() - 1.0
994 }
995 };
996 let peak_floor = m / (s * s) + 1.0;
1000 let required = asymptotic_index(coeff_mean)
1001 .max(asymptotic_index(coeff_deriv))
1002 .max(peak_floor);
1003 if !required.is_finite() || required > LOGIT_MAX_TERMS as f64 {
1004 return None;
1005 }
1006 Some((required.ceil() as usize).max(4))
1009}
1010
1011#[inline]
1012fn stable_sigmoidwith_derivative(x: f64) -> (f64, f64) {
1013 let x_clamped = x.clamp(-QUADRATURE_EXP_LOG_MAX, QUADRATURE_EXP_LOG_MAX);
1014 if x_clamped != x {
1015 return (sigmoid(x), 0.0);
1016 }
1017 if x_clamped >= 0.0 {
1018 let z = (-x_clamped).exp();
1019 let denom = 1.0 + z;
1020 (1.0 / denom, z / (denom * denom))
1021 } else {
1022 let z = x_clamped.exp();
1023 let denom = 1.0 + z;
1024 (z / denom, z / (denom * denom))
1025 }
1026}
1027
1028#[inline]
1029fn logit_small_sigma_taylor(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1030 let (mean0, d1, d2, d3) = component_point_jet(LinkComponent::Logit, mu);
1038 let s2 = sigma * sigma;
1039 IntegratedMeanDerivative {
1040 mean: (mean0 + 0.5 * s2 * d2).clamp(0.0, 1.0),
1041 dmean_dmu: (d1 + 0.5 * s2 * d3).max(0.0),
1042 mode: IntegratedExpectationMode::ControlledAsymptotic,
1043 }
1044}
1045
1046#[inline]
1047fn logit_tail_asymptotic(mu: f64, sigma: f64) -> Option<IntegratedMeanDerivative> {
1048 if mu <= 0.0 {
1053 let log_mean = mu + 0.5 * sigma * sigma;
1054 if log_mean <= LOGIT_TAIL_LOG_MAX {
1055 let mean = safe_exp(log_mean);
1056 return Some(IntegratedMeanDerivative {
1057 mean,
1058 dmean_dmu: mean,
1059 mode: IntegratedExpectationMode::ControlledAsymptotic,
1060 });
1061 }
1062 } else {
1063 let log_tail = -mu + 0.5 * sigma * sigma;
1064 if log_tail <= LOGIT_TAIL_LOG_MAX {
1065 let tail = safe_exp(log_tail);
1066 return Some(IntegratedMeanDerivative {
1067 mean: 1.0 - tail,
1068 dmean_dmu: tail,
1069 mode: IntegratedExpectationMode::ControlledAsymptotic,
1070 });
1071 }
1072 }
1073 None
1074}
1075
1076#[inline]
1077fn scaled_erfcx_termwith_derivative(m: f64, s: f64, x: f64, dxdm: f64) -> (f64, f64) {
1078 let pref = 0.5 * (-(m * m) / (2.0 * s * s)).exp();
1079 if x >= 0.0 {
1080 let ex = erfcx_nonnegative(x);
1081 let term = pref * ex;
1082 let ex_prime = 2.0 * x * ex - std::f64::consts::FRAC_2_SQRT_PI;
1083 let dterm = pref * ((-m / (s * s)) * ex + ex_prime * dxdm);
1084 (term, dterm)
1085 } else {
1086 let lead = (x * x - (m * m) / (2.0 * s * s)).exp();
1087 let dlead = lead * (2.0 * x * dxdm - m / (s * s));
1088 let (rest, drest) = scaled_erfcx_termwith_derivative(m, s, -x, -dxdm);
1089 (lead - rest, dlead - drest)
1090 }
1091}
1092
1093pub(crate) fn logit_posterior_meanwith_deriv_exact(
1094 mu: f64,
1095 sigma: f64,
1096) -> Result<IntegratedMeanDerivative, EstimationError> {
1097 if !(mu.is_finite() && sigma.is_finite()) {
1115 crate::bail_invalid_estim!("logit exact expectation requires finite mu and sigma");
1116 }
1117 if sigma <= LOGIT_SIGMA_DEGENERATE {
1118 let (mean, dmean_dmu) = stable_sigmoidwith_derivative(mu);
1119 return Ok(IntegratedMeanDerivative {
1120 mean,
1121 dmean_dmu,
1122 mode: IntegratedExpectationMode::ExactClosedForm,
1123 });
1124 }
1125 if let Some(out) = logit_tail_asymptotic(mu, sigma) {
1126 return Ok(out);
1127 }
1128 if sigma < LOGIT_SIGMA_TAYLOR_MAX {
1129 return Ok(logit_small_sigma_taylor(mu, sigma));
1130 }
1131 if logistic_normal_exact_eligible(mu, sigma)
1132 && let Ok(out) = logit_posterior_meanwith_deriv_exact_erfcx(mu, sigma)
1133 {
1134 return Ok(out);
1135 }
1136 Err(EstimationError::InvalidInput(
1145 "logit analytic expectation has no certified representation in this regime".to_string(),
1146 ))
1147}
1148
1149fn logit_posterior_meanwith_deriv_exact_erfcx(
1150 mu: f64,
1151 sigma: f64,
1152) -> Result<IntegratedMeanDerivative, EstimationError> {
1153 let m = mu.abs();
1179 let s = sigma;
1180 let z = SQRT_2 * s;
1181 let phi_term = gam_math::probability::normal_cdf(m / s);
1182 let phi_prime = gam_math::probability::normal_pdf(m / s) / s;
1183 let Some(max_k) = logistic_normal_series_cutoff(mu, sigma, LOGIT_ERFCX_ACCURACY_TARGET) else {
1184 crate::bail_invalid_estim!(
1185 "logit erfcx series truncation bound exceeds LOGIT_MAX_TERMS at the required accuracy"
1186 .to_string(),
1187 );
1188 };
1189
1190 let mut sum = 0.0_f64;
1191 let mut dsum = 0.0_f64;
1192 let mut k = 1usize;
1200 while k <= max_k {
1201 for kk in [k, k + 1].into_iter().filter(|kk| *kk <= max_k) {
1202 let kf = kk as f64;
1203 let a = (kf * s * s + m) / z;
1204 let b = (kf * s * s - m) / z;
1205 let sign = if kk % 2 == 1 { 1.0 } else { -1.0 };
1206 let (va, dva) = scaled_erfcx_termwith_derivative(m, s, a, 1.0 / z);
1207 let (vb, dvb) = scaled_erfcx_termwith_derivative(m, s, b, -1.0 / z);
1208 sum += sign * (va - vb);
1209 dsum += sign * (dva - dvb);
1210 }
1211 k += 2;
1212 }
1213
1214 let mut mean = phi_term + sum;
1215 let dmean = (phi_prime + dsum).max(0.0);
1216 if mu < 0.0 {
1217 mean = 1.0 - mean;
1218 }
1219 if !(mean.is_finite() && dmean.is_finite() && dmean >= 0.0) {
1220 crate::bail_invalid_estim!("logit erfcx expectation produced non-finite values");
1221 }
1222 Ok(IntegratedMeanDerivative {
1223 mean,
1224 dmean_dmu: dmean,
1225 mode: IntegratedExpectationMode::ExactSpecialFunction,
1226 })
1227}
1228
1229#[inline]
1234fn logit_posterior_meanwith_deriv_quadrature(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1235 let mean = integrate_normal_adaptive(mu, sigma, |x| stable_sigmoidwith_derivative(x).0);
1236 let dmean_dmu =
1237 integrate_normal_adaptive(mu, sigma, |x| stable_sigmoidwith_derivative(x).1).max(0.0);
1238 IntegratedMeanDerivative {
1239 mean,
1240 dmean_dmu,
1241 mode: IntegratedExpectationMode::QuadratureFallback,
1242 }
1243}
1244
1245#[inline]
1246fn logit_posterior_meanwith_deriv_controlled(
1247 mu: f64,
1248 sigma: f64,
1249) -> Result<IntegratedMeanDerivative, EstimationError> {
1250 if !(mu.is_finite() && sigma.is_finite()) {
1251 crate::bail_invalid_estim!("logit integrated moments require finite mu and sigma");
1252 }
1253 let candidate = match logit_posterior_meanwith_deriv_exact(mu, sigma) {
1254 Ok(out) => out,
1255 Err(_) => return Ok(logit_posterior_meanwith_deriv_quadrature(mu, sigma)),
1256 };
1257 match candidate.mode {
1269 IntegratedExpectationMode::ExactSpecialFunction
1270 | IntegratedExpectationMode::ControlledAsymptotic => {
1271 let reference = logit_posterior_meanwith_deriv_quadrature(mu, sigma);
1272 if integrated_mean_derivative_drift_exceeds(
1273 &candidate, &reference, 1e-6, 1e-4, 1e-7, 1e-3,
1274 ) {
1275 Ok(reference)
1276 } else {
1277 Ok(candidate)
1278 }
1279 }
1280 _ => Ok(candidate),
1281 }
1282}
1283
1284#[inline]
1285fn log_normal_cdf_stable(x: f64) -> f64 {
1286 if !x.is_finite() {
1287 return if x.is_sign_negative() {
1288 f64::NEG_INFINITY
1289 } else {
1290 0.0
1291 };
1292 }
1293 if x < -8.0 {
1294 let u = -x / SQRT_2;
1295 -u * u + (0.5 * erfcx_nonnegative(u)).ln()
1296 } else {
1297 gam_math::probability::normal_cdf(x).ln()
1301 }
1302}
1303
1304#[inline]
1305fn cloglog_extreme_asymptotic(mu: f64, sigma: f64) -> Option<IntegratedMeanDerivative> {
1306 let rare_log = mu + 0.5 * sigma * sigma;
1316 if rare_log <= CLOGLOG_RARE_EVENT_LOG_MAX {
1317 let mean = safe_exp(rare_log);
1318 return Some(IntegratedMeanDerivative {
1319 mean,
1320 dmean_dmu: mean,
1321 mode: IntegratedExpectationMode::ControlledAsymptotic,
1322 });
1323 }
1324 if mu - CLOGLOG_POSITIVE_SATURATION_SIGMAS * sigma >= CLOGLOG_POSITIVE_SATURATION_EDGE {
1325 return Some(IntegratedMeanDerivative {
1326 mean: 1.0,
1327 dmean_dmu: 0.0,
1328 mode: IntegratedExpectationMode::ControlledAsymptotic,
1329 });
1330 }
1331 None
1337}
1338
1339#[inline]
1340fn cloglog_survival_extreme_asymptotic(
1341 mu: f64,
1342 sigma: f64,
1343) -> Option<(f64, IntegratedExpectationMode)> {
1344 let rare_log = mu + 0.5 * sigma * sigma;
1345 if rare_log <= CLOGLOG_RARE_EVENT_LOG_MAX {
1346 let mean = safe_exp(rare_log);
1347 return Some((
1348 (1.0 - mean).clamp(0.0, 1.0),
1349 IntegratedExpectationMode::ControlledAsymptotic,
1350 ));
1351 }
1352 if mu - CLOGLOG_POSITIVE_SATURATION_SIGMAS * sigma >= CLOGLOG_POSITIVE_SATURATION_EDGE {
1353 return Some((0.0, IntegratedExpectationMode::ControlledAsymptotic));
1358 }
1359 None
1363}
1364
1365#[inline]
1373fn cloglog_gumbel_quad_nodes(sigma: f64) -> usize {
1374 let target = (CLOGLOG_GUMBEL_QUAD_NODE_SCALE / sigma.min(CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN))
1375 .ceil() as usize;
1376 let n = target
1377 .max(CLOGLOG_GUMBEL_QUAD_MIN_NODES)
1378 .min(CLOGLOG_GUMBEL_QUAD_MAX_NODES);
1379 if n % 2 == 0 { n + 1 } else { n }
1380}
1381
1382fn cloglog_log_survival_gumbel_quadrature(ctx: &QuadratureContext, mu: f64, sigma: f64) -> f64 {
1405 let a = CLOGLOG_GUMBEL_QUAD_ETA_LO;
1406 let b = CLOGLOG_GUMBEL_QUAD_ETA_HI;
1407 let half = 0.5 * (b - a);
1408 let mid = 0.5 * (a + b);
1409 let rule = ctx.clenshaw_curtis_n(cloglog_gumbel_quad_nodes(sigma));
1410 let mut running_max = f64::NEG_INFINITY;
1413 let mut running_sum = 0.0_f64;
1414 for (&node, &weight) in rule.nodes.iter().zip(rule.weights.iter()) {
1415 let eta = half * node + mid;
1416 let summand = (weight * half).ln()
1417 + (eta - safe_exp(eta))
1418 + log_normal_cdf_stable((eta - mu) / sigma);
1419 if !summand.is_finite() {
1420 continue;
1421 }
1422 if summand > running_max {
1423 running_sum = running_sum * (running_max - summand).exp() + 1.0;
1424 running_max = summand;
1425 } else {
1426 running_sum += (summand - running_max).exp();
1427 }
1428 }
1429 if running_max == f64::NEG_INFINITY {
1430 f64::NEG_INFINITY
1431 } else {
1432 running_max + running_sum.ln()
1433 }
1434}
1435
1436pub(crate) fn cloglog_log_survival_term_controlled(
1445 ctx: &QuadratureContext,
1446 mu: f64,
1447 sigma: f64,
1448) -> (f64, IntegratedExpectationMode) {
1449 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
1450 return (-safe_exp(mu), IntegratedExpectationMode::ExactClosedForm);
1452 }
1453 let rare_log = mu + 0.5 * sigma * sigma;
1454 if rare_log <= CLOGLOG_RARE_EVENT_LOG_MAX {
1455 return (
1458 (-safe_exp(rare_log)).ln_1p(),
1459 IntegratedExpectationMode::ControlledAsymptotic,
1460 );
1461 }
1462 if sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN {
1463 return (
1464 cloglog_log_survival_gumbel_quadrature(ctx, mu, sigma),
1465 IntegratedExpectationMode::ControlledAsymptotic,
1466 );
1467 }
1468 let (value, mode) = cloglog_survival_term_controlled(ctx, mu, sigma);
1469 if value > 0.0 {
1470 (value.ln(), mode)
1471 } else {
1472 (
1475 cloglog_log_survival_gumbel_quadrature(ctx, mu, sigma),
1476 IntegratedExpectationMode::QuadratureFallback,
1477 )
1478 }
1479}
1480
1481#[inline]
1500fn gumbel_survival(x: f64) -> f64 {
1501 (-safe_exp(x)).exp()
1502}
1503
1504#[inline]
1509fn cloglog_mean_d1_exact(x: f64) -> f64 {
1510 let ex = safe_exp(x);
1511 if ex.is_infinite() {
1512 0.0
1513 } else {
1514 ex * (-ex).exp()
1515 }
1516}
1517
1518#[inline]
1528fn cloglog_mean_exact(x: f64) -> f64 {
1529 cloglog_negative_tail_mean(x)
1530}
1531
1532#[inline]
1548fn cloglog_negative_tail_mean(eta: f64) -> f64 {
1549 if eta < -745.0 {
1553 0.0
1555 } else {
1556 let ex = safe_exp(eta);
1559 -(-ex).exp_m1()
1560 }
1561}
1562
1563#[inline]
1568fn cloglog_small_sigma_taylor(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1569 if sigma <= CLOGLOG_SIGMA_DEGENERATE {
1593 return IntegratedMeanDerivative {
1594 mean: cloglog_mean_exact(mu),
1595 dmean_dmu: cloglog_mean_d1_exact(mu),
1596 mode: IntegratedExpectationMode::ExactClosedForm,
1597 };
1598 }
1599
1600 let ex = safe_exp(mu);
1601 if !ex.is_finite() {
1602 return IntegratedMeanDerivative {
1604 mean: 1.0,
1605 dmean_dmu: 0.0,
1606 mode: IntegratedExpectationMode::ControlledAsymptotic,
1607 };
1608 }
1609 let surv = (-ex).exp();
1610 if surv == 0.0 {
1611 return IntegratedMeanDerivative {
1613 mean: 1.0,
1614 dmean_dmu: 0.0,
1615 mode: IntegratedExpectationMode::ControlledAsymptotic,
1616 };
1617 }
1618
1619 let s2 = sigma * sigma;
1620 let s4 = s2 * s2;
1621 let s6 = s4 * s2;
1622 let s8 = s4 * s4;
1623 let e2x = ex * ex;
1624 let e3x = e2x * ex;
1625 let e4x = e3x * ex;
1626 let e5x = e4x * ex;
1627 let e6x = e5x * ex;
1628 let e7x = e6x * ex;
1629 let e8x = e7x * ex;
1630 let e9x = e8x * ex;
1631 let f0 = -(-ex).exp_m1();
1645 let f1 = ex * surv;
1646 let f2 = surv * (ex - e2x);
1647 let f3 = surv * (ex - 3.0 * e2x + e3x);
1648 let f4 = surv * (ex - 7.0 * e2x + 6.0 * e3x - e4x);
1649 let f5 = surv * (ex - 15.0 * e2x + 25.0 * e3x - 10.0 * e4x + e5x);
1650 let f6 = surv * (ex - 31.0 * e2x + 90.0 * e3x - 65.0 * e4x + 15.0 * e5x - e6x);
1651 let f7 = surv * (ex - 63.0 * e2x + 301.0 * e3x - 350.0 * e4x + 140.0 * e5x - 21.0 * e6x + e7x);
1652 let f8 = surv
1653 * (ex - 127.0 * e2x + 966.0 * e3x - 1701.0 * e4x + 1050.0 * e5x - 266.0 * e6x + 28.0 * e7x
1654 - e8x);
1655 let f9 = surv
1656 * (ex - 255.0 * e2x + 3025.0 * e3x - 7770.0 * e4x + 6951.0 * e5x - 2646.0 * e6x
1657 + 462.0 * e7x
1658 - 36.0 * e8x
1659 + e9x);
1660 IntegratedMeanDerivative {
1664 mean: f0 + 0.5 * s2 * f2 + (s4 / 8.0) * f4 + (s6 / 48.0) * f6 + (s8 / 384.0) * f8,
1665 dmean_dmu: (f1 + 0.5 * s2 * f3 + (s4 / 8.0) * f5 + (s6 / 48.0) * f7 + (s8 / 384.0) * f9)
1666 .max(0.0),
1667 mode: IntegratedExpectationMode::ControlledAsymptotic,
1668 }
1669}
1670
1671#[inline]
1672fn adaptive_simpson_refine(
1677 g: &impl Fn(f64) -> f64,
1678 a: f64,
1679 b: f64,
1680 fa: f64,
1681 fb: f64,
1682 fm: f64,
1683 whole: f64,
1684 tol: f64,
1685 depth: i32,
1686) -> f64 {
1687 let m = 0.5 * (a + b);
1688 let lm = 0.5 * (a + m);
1689 let rm = 0.5 * (m + b);
1690 let flm = g(lm);
1691 let frm = g(rm);
1692 let left = (m - a) / 6.0 * (fa + 4.0 * flm + fm);
1693 let right = (b - m) / 6.0 * (fm + 4.0 * frm + fb);
1694 let est = left + right;
1695 if depth <= 0 || (est - whole).abs() <= 15.0 * tol {
1696 return est + (est - whole) / 15.0;
1697 }
1698 adaptive_simpson_refine(g, a, m, fa, fm, flm, left, 0.5 * tol, depth - 1)
1699 + adaptive_simpson_refine(g, m, b, fm, fb, frm, right, 0.5 * tol, depth - 1)
1700}
1701
1702fn integrate_normal_adaptive(mu: f64, sigma: f64, f: impl Fn(f64) -> f64) -> f64 {
1717 if !(sigma.is_finite()) || sigma < 1e-10 {
1718 return f(mu);
1719 }
1720 const K: f64 = 15.0;
1721 const INITIAL_PANELS: usize = 24;
1722 const TOL: f64 = 1e-12;
1723 const MAX_DEPTH: i32 = 40;
1724 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
1725 let g = |u: f64| f(mu + sigma * u) * inv_sqrt_2pi * (-0.5 * u * u).exp();
1729 let panel = 2.0 * K / INITIAL_PANELS as f64;
1730 let mut total = 0.0;
1731 for p in 0..INITIAL_PANELS {
1732 let a = -K + p as f64 * panel;
1733 let b = a + panel;
1734 let fa = g(a);
1735 let fb = g(b);
1736 let fm = g(0.5 * (a + b));
1737 let whole = (b - a) / 6.0 * (fa + 4.0 * fm + fb);
1738 total += adaptive_simpson_refine(&g, a, b, fa, fb, fm, whole, TOL, MAX_DEPTH);
1739 }
1740 total
1741}
1742
1743fn cloglog_posterior_meanwith_deriv_quadrature(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1744 if sigma < 1e-10 {
1745 return IntegratedMeanDerivative {
1746 mean: cloglog_mean_exact(mu),
1747 dmean_dmu: cloglog_mean_d1_exact(mu),
1748 mode: IntegratedExpectationMode::ExactClosedForm,
1749 };
1750 }
1751 let mean = cloglog_mean_from_survival(survival_posterior_mean_quadrature(mu, sigma));
1752 let dmean_dmu = integrate_normal_adaptive(mu, sigma, cloglog_mean_d1_exact).max(0.0);
1753 IntegratedMeanDerivative {
1754 mean,
1755 dmean_dmu,
1756 mode: IntegratedExpectationMode::QuadratureFallback,
1757 }
1758}
1759
1760#[inline]
1761fn survival_posterior_mean_quadrature(eta: f64, se_eta: f64) -> f64 {
1762 integrate_normal_adaptive(eta, se_eta, gumbel_survival).clamp(0.0, 1.0)
1763}
1764
1765fn cloglog_survival_term_controlled(
1766 ctx: &QuadratureContext,
1767 mu: f64,
1768 sigma: f64,
1769) -> (f64, IntegratedExpectationMode) {
1770 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
1807 return (
1808 gumbel_survival(mu).clamp(0.0, 1.0),
1809 IntegratedExpectationMode::ExactClosedForm,
1810 );
1811 }
1812 if sigma < CLOGLOG_SIGMA_TAYLOR_MAX {
1813 let mean = cloglog_small_sigma_taylor(mu, sigma).mean;
1814 return (
1815 (1.0 - mean).clamp(0.0, 1.0),
1816 IntegratedExpectationMode::ControlledAsymptotic,
1817 );
1818 }
1819 if let Some(out) = cloglog_survival_extreme_asymptotic(mu, sigma) {
1820 return out;
1821 }
1822 if sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN {
1823 let log_s = cloglog_log_survival_gumbel_quadrature(ctx, mu, sigma);
1830 return (
1831 safe_exp(log_s).clamp(0.0, 1.0),
1832 IntegratedExpectationMode::ControlledAsymptotic,
1833 );
1834 }
1835 if cloglog_survival_miles_is_reliable(mu, sigma)
1836 && let Ok(out) = cloglog_survival_miles(mu, sigma)
1837 {
1838 return (
1839 out.clamp(0.0, 1.0),
1840 IntegratedExpectationMode::ExactSpecialFunction,
1841 );
1842 }
1843 if cloglog_should_prefer_cc(mu, sigma, CLOGLOG_CC_TOL)
1844 && let Ok(out) = cloglog_survival_cc(ctx, mu, sigma, CLOGLOG_CC_TOL)
1845 {
1846 return (
1847 out.clamp(0.0, 1.0),
1848 IntegratedExpectationMode::ExactSpecialFunction,
1849 );
1850 }
1851 if let Ok(out) = cloglog_survival_gamma_reference(mu, sigma) {
1852 return (
1853 out.clamp(0.0, 1.0),
1854 IntegratedExpectationMode::ExactSpecialFunction,
1855 );
1856 }
1857 (
1858 survival_posterior_mean_quadrature(mu, sigma),
1859 IntegratedExpectationMode::QuadratureFallback,
1860 )
1861}
1862
1863#[inline]
1864fn lognormal_laplace_term_controlled(
1865 ctx: &QuadratureContext,
1866 z: f64,
1867 mu: f64,
1868 sigma: f64,
1869) -> (f64, IntegratedExpectationMode) {
1870 if !(z.is_finite() && z > 0.0) {
1896 return (f64::NAN, IntegratedExpectationMode::QuadratureFallback);
1897 }
1898 lognormal_laplace_unit_term_shared(ctx, mu + z.ln(), sigma)
1899}
1900
1901#[inline]
1902pub(crate) fn lognormal_laplace_unit_term_shared(
1903 ctx: &QuadratureContext,
1904 shifted_mu: f64,
1905 sigma: f64,
1906) -> (f64, IntegratedExpectationMode) {
1907 cloglog_survival_term_controlled(ctx, shifted_mu, sigma)
1908}
1909
1910#[inline]
1914pub fn lognormal_laplace_unit_log_term_shared(
1915 ctx: &QuadratureContext,
1916 shifted_mu: f64,
1917 sigma: f64,
1918) -> (f64, IntegratedExpectationMode) {
1919 cloglog_log_survival_term_controlled(ctx, shifted_mu, sigma)
1920}
1921
1922#[inline]
1923fn cloglog_survivalsecond_moment_controlled(
1924 ctx: &QuadratureContext,
1925 mu: f64,
1926 sigma: f64,
1927) -> (f64, IntegratedExpectationMode) {
1928 lognormal_laplace_term_controlled(ctx, 2.0, mu, sigma)
1943}
1944
1945#[inline]
1946fn cloglog_survival_pair_controlled(
1947 ctx: &QuadratureContext,
1948 mu: f64,
1949 sigma: f64,
1950) -> (
1951 (f64, IntegratedExpectationMode),
1952 (f64, IntegratedExpectationMode),
1953) {
1954 let shiftedmu = mu + sigma * sigma;
1955
1956 if cloglog_survival_miles_is_reliable(mu, sigma)
1967 && cloglog_survival_miles_is_reliable(shiftedmu, sigma)
1968 && let (Ok(base), Ok(shifted)) = (
1969 cloglog_survival_miles(mu, sigma),
1970 cloglog_survival_miles(shiftedmu, sigma),
1971 )
1972 {
1973 return (
1974 (
1975 base.clamp(0.0, 1.0),
1976 IntegratedExpectationMode::ExactSpecialFunction,
1977 ),
1978 (
1979 shifted.clamp(0.0, 1.0),
1980 IntegratedExpectationMode::ExactSpecialFunction,
1981 ),
1982 );
1983 }
1984
1985 if cloglog_should_prefer_cc(mu, sigma, CLOGLOG_CC_TOL)
1986 && cloglog_should_prefer_cc(shiftedmu, sigma, CLOGLOG_CC_TOL)
1987 && let (Ok(base), Ok(shifted)) = (
1988 cloglog_survival_cc(ctx, mu, sigma, CLOGLOG_CC_TOL),
1989 cloglog_survival_cc(ctx, shiftedmu, sigma, CLOGLOG_CC_TOL),
1990 )
1991 {
1992 return (
1993 (
1994 base.clamp(0.0, 1.0),
1995 IntegratedExpectationMode::ExactSpecialFunction,
1996 ),
1997 (
1998 shifted.clamp(0.0, 1.0),
1999 IntegratedExpectationMode::ExactSpecialFunction,
2000 ),
2001 );
2002 }
2003
2004 if let (Ok(base), Ok(shifted)) = (
2005 cloglog_survival_gamma_reference(mu, sigma),
2006 cloglog_survival_gamma_reference(shiftedmu, sigma),
2007 ) {
2008 return (
2009 (
2010 base.clamp(0.0, 1.0),
2011 IntegratedExpectationMode::ExactSpecialFunction,
2012 ),
2013 (
2014 shifted.clamp(0.0, 1.0),
2015 IntegratedExpectationMode::ExactSpecialFunction,
2016 ),
2017 );
2018 }
2019
2020 (
2021 cloglog_survival_term_controlled(ctx, mu, sigma),
2022 cloglog_survival_term_controlled(ctx, shiftedmu, sigma),
2023 )
2024}
2025
2026#[inline]
2027fn cloglog_mean_from_survival(survival: f64) -> f64 {
2028 let survival = survival.clamp(0.0, 1.0);
2029 if survival > 0.5 {
2030 -survival.ln().exp_m1()
2040 } else {
2041 1.0 - survival
2042 }
2043}
2044
2045#[inline]
2046fn cloglog_shift_identity_derivative(mu: f64, sigma: f64, shifted_survival: f64) -> f64 {
2047 if !(mu.is_finite() && sigma.is_finite()) || shifted_survival <= 0.0 {
2061 return 0.0;
2062 }
2063 cloglog_shift_identity_derivative_log(mu, sigma, shifted_survival.ln())
2064}
2065
2066#[inline]
2076fn cloglog_shift_identity_derivative_log(mu: f64, sigma: f64, log_shifted_survival: f64) -> f64 {
2077 if !(mu.is_finite() && sigma.is_finite()) || log_shifted_survival == f64::NEG_INFINITY {
2078 return 0.0;
2079 }
2080 let log_derivative = mu + 0.5 * sigma * sigma + log_shifted_survival;
2081 let upper = 1.0 / std::f64::consts::E;
2082 if !log_derivative.is_finite() {
2083 return upper;
2086 }
2087 safe_exp(log_derivative).clamp(0.0, upper)
2088}
2089
2090#[inline]
2091fn log_half_erfc_stable(u: f64) -> f64 {
2092 if u > 0.0 {
2104 -u * u + (0.5 * erfcx_nonnegative(u)).ln()
2105 } else {
2106 normal_logcdf(-u * SQRT_2)
2107 }
2108}
2109
2110#[inline]
2150fn cloglog_survival_miles_is_reliable(mu: f64, sigma: f64) -> bool {
2151 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0) {
2152 return false;
2153 }
2154 let alpha_ln = CLOGLOG_MILES_ALPHA.ln();
2155 let shifted = mu - alpha_ln;
2156 let peak_log = CLOGLOG_MILES_ALPHA - 0.5 * shifted * shifted / (sigma * sigma);
2157 peak_log.is_finite() && peak_log <= CLOGLOG_MILES_PEAK_LOG_MAX
2158}
2159
2160fn cloglog_survival_miles(mu: f64, sigma: f64) -> Result<f64, EstimationError> {
2161 let alpha_ln = CLOGLOG_MILES_ALPHA.ln();
2191 let mut s_sum = 0.0_f64;
2192 let mut stable_pairs = 0usize;
2193
2194 for pair_start in (0..CLOGLOG_MILES_MAX_TERMS).step_by(2) {
2195 let mut pair_s = 0.0_f64;
2196 for n in pair_start..(pair_start + 2).min(CLOGLOG_MILES_MAX_TERMS) {
2197 let nf = n as f64;
2198 let sign = if n % 2 == 0 { 1.0 } else { -1.0 };
2199 let base_log = nf * mu + 0.5 * sigma * sigma * nf * nf
2200 - statrs::function::gamma::ln_gamma(nf + 1.0);
2201 let u = (mu - alpha_ln + sigma * sigma * nf) / (SQRT_2 * sigma);
2202 let log_half_erfc = log_half_erfc_stable(u);
2203 let term_log = base_log + log_half_erfc;
2204 if term_log > QUADRATURE_EXP_LOG_MAX {
2205 crate::bail_invalid_estim!("Miles cloglog series term exceeded finite exp range");
2206 }
2207 let term = sign * safe_exp(term_log);
2208 pair_s += term;
2209 }
2210 s_sum += pair_s;
2211
2212 let s_scale = s_sum.abs().max(1.0);
2213 if pair_s.abs() <= 2e-15 * s_scale {
2214 stable_pairs += 1;
2215 if stable_pairs >= SERIES_CONSECUTIVE_SMALL_TERMS {
2216 if s_sum.is_finite() && (-1e-10..=1.0 + 1e-10).contains(&s_sum) {
2217 return Ok(s_sum.clamp(0.0, 1.0));
2218 }
2219 break;
2220 }
2221 } else {
2222 stable_pairs = 0;
2223 }
2224 }
2225
2226 Err(EstimationError::InvalidInput(
2227 "Miles cloglog series did not converge safely".to_string(),
2228 ))
2229}
2230
2231fn cloglog_survival_cc(
2232 ctx: &QuadratureContext,
2233 mu: f64,
2234 sigma: f64,
2235 tol: f64,
2236) -> Result<f64, EstimationError> {
2237 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0 && tol.is_finite() && tol > 0.0) {
2238 crate::bail_invalid_estim!(
2239 "CC cloglog backend requires finite mu, positive sigma, and positive tolerance"
2240 .to_string(),
2241 );
2242 }
2243
2244 let p_tail = (tol / 8.0).clamp(1e-300, 0.25);
2282 let a = gam_math::probability::standard_normal_quantile(p_tail)
2283 .map(|z| -z)
2284 .unwrap_or(8.0)
2285 .max(1.0);
2286 let n = cloglog_cc_required_nodes(mu, sigma, tol)?;
2287 if n > CLOGLOG_CC_NODE_CAP {
2288 crate::bail_invalid_estim!("CC cloglog backend requires too many nodes");
2289 }
2290
2291 let rule = ctx.clenshaw_curtis_n(n);
2292 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
2293 let mut sum = 0.0_f64;
2294 let mut c = 0.0_f64;
2295 for (&x, &w) in rule.nodes.iter().zip(rule.weights.iter()) {
2296 let t = a * x;
2297 let u = mu + sigma * t;
2298 let e = safe_exp(u);
2299 let w0 = (-0.5 * t * t).exp() * inv_sqrt_2pi;
2300 let yk = w * w0 * (-e).exp() - c;
2301 let tk = sum + yk;
2302 c = (tk - sum) - yk;
2303 sum = tk;
2304 }
2305
2306 let survival = (a * sum).clamp(0.0, 1.0);
2307 if !survival.is_finite() {
2308 crate::bail_invalid_estim!("CC cloglog backend produced non-finite values");
2309 }
2310 Ok(survival)
2311}
2312
2313#[inline]
2314fn complex_add(a: Complex, b: Complex) -> Complex {
2315 Complex {
2316 re: a.re + b.re,
2317 im: a.im + b.im,
2318 }
2319}
2320
2321#[inline]
2322fn complex_sub(a: Complex, b: Complex) -> Complex {
2323 Complex {
2324 re: a.re - b.re,
2325 im: a.im - b.im,
2326 }
2327}
2328
2329#[inline]
2330fn complexmul(a: Complex, b: Complex) -> Complex {
2331 Complex {
2332 re: a.re * b.re - a.im * b.im,
2333 im: a.re * b.im + a.im * b.re,
2334 }
2335}
2336
2337#[inline]
2338fn complex_div(a: Complex, b: Complex) -> Complex {
2339 let den = (b.re * b.re + b.im * b.im).max(1e-300);
2340 Complex {
2341 re: (a.re * b.re + a.im * b.im) / den,
2342 im: (a.im * b.re - a.re * b.im) / den,
2343 }
2344}
2345
2346#[inline]
2347fn complex_abs(z: Complex) -> f64 {
2348 z.re.hypot(z.im)
2349}
2350
2351#[inline]
2352fn complex_ln(z: Complex) -> Complex {
2353 Complex {
2354 re: complex_abs(z).ln(),
2355 im: z.im.atan2(z.re),
2356 }
2357}
2358
2359#[inline]
2360fn complex_exp(z: Complex) -> Complex {
2361 let e = z.re.exp();
2362 Complex {
2363 re: e * z.im.cos(),
2364 im: e * z.im.sin(),
2365 }
2366}
2367
2368#[inline]
2369fn complex_sin(z: Complex) -> Complex {
2370 Complex {
2371 re: z.re.sin() * z.im.cosh(),
2372 im: z.re.cos() * z.im.sinh(),
2373 }
2374}
2375
2376fn complex_log_gamma_lanczos(z: Complex) -> Complex {
2377 const G: f64 = 7.0;
2381 const COEFFS: [f64; 9] = [
2382 0.999_999_999_999_809_9,
2383 676.520_368_121_885_1,
2384 -1_259.139_216_722_402_8,
2385 771.323_428_777_653_1,
2386 -176.615_029_162_140_6,
2387 12.507_343_278_686_905,
2388 -0.138_571_095_265_720_12,
2389 9.984_369_578_019_572e-6,
2390 1.505_632_735_149_311_6e-7,
2391 ];
2392
2393 if z.re < 0.5 {
2394 let piz = Complex {
2395 re: std::f64::consts::PI * z.re,
2396 im: std::f64::consts::PI * z.im,
2397 };
2398 let one_minusz = Complex {
2399 re: 1.0 - z.re,
2400 im: -z.im,
2401 };
2402 return complex_sub(
2403 complex_sub(
2404 Complex {
2405 re: std::f64::consts::PI.ln(),
2406 im: 0.0,
2407 },
2408 complex_ln(complex_sin(piz)),
2409 ),
2410 complex_log_gamma_lanczos(one_minusz),
2411 );
2412 }
2413
2414 let z1 = Complex {
2415 re: z.re - 1.0,
2416 im: z.im,
2417 };
2418 let mut x = Complex {
2419 re: COEFFS[0],
2420 im: 0.0,
2421 };
2422 for (i, c) in COEFFS.iter().enumerate().skip(1) {
2423 x = complex_add(
2424 x,
2425 complex_div(
2426 Complex { re: *c, im: 0.0 },
2427 Complex {
2428 re: z1.re + i as f64,
2429 im: z1.im,
2430 },
2431 ),
2432 );
2433 }
2434 let t = Complex {
2435 re: z1.re + G + 0.5,
2436 im: z1.im,
2437 };
2438 complex_add(
2439 complex_add(
2440 Complex {
2441 re: 0.5 * (2.0 * std::f64::consts::PI).ln(),
2442 im: 0.0,
2443 },
2444 complexmul(
2445 Complex {
2446 re: z1.re + 0.5,
2447 im: z1.im,
2448 },
2449 complex_ln(t),
2450 ),
2451 ),
2452 complex_sub(complex_ln(x), t),
2453 )
2454}
2455
2456fn cloglog_survival_gamma_reference(mu: f64, sigma: f64) -> Result<f64, EstimationError> {
2460 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 0.0 {
2461 crate::bail_invalid_estim!(
2462 "Gamma cloglog reference backend requires finite mu and positive sigma"
2463 );
2464 }
2465
2466 let n = (CLOGLOG_GAMMA_T_MAX_REF / CLOGLOG_GAMMA_H_REF).round() as usize;
2500 let n = if n.is_multiple_of(2) { n } else { n + 1 };
2501 let h = CLOGLOG_GAMMA_T_MAX_REF / n as f64;
2502
2503 let eval = |t: f64| -> f64 {
2504 let z = Complex {
2505 re: CLOGLOG_GAMMA_K_REF,
2506 im: t,
2507 };
2508 let log_gamma = complex_log_gamma_lanczos(z);
2509 let z_sq = complexmul(z, z);
2510 let exponent = complex_sub(
2511 complex_add(
2512 log_gamma,
2513 Complex {
2514 re: 0.5 * sigma * sigma * z_sq.re,
2515 im: 0.5 * sigma * sigma * z_sq.im,
2516 },
2517 ),
2518 Complex {
2519 re: mu * z.re,
2520 im: mu * z.im,
2521 },
2522 );
2523 complex_exp(exponent).re
2524 };
2525
2526 let f0 = eval(0.0);
2527 let fn_ = eval(CLOGLOG_GAMMA_T_MAX_REF);
2528 let mut sum_s = f0 + fn_;
2529 for i in 1..n {
2530 let t = i as f64 * h;
2531 let fi = eval(t);
2532 let w = if i % 2 == 0 { 2.0 } else { 4.0 };
2533 sum_s += w * fi;
2534 }
2535 let sval = ((h / 3.0) * sum_s / std::f64::consts::PI).clamp(0.0, 1.0);
2536 if !sval.is_finite() {
2537 crate::bail_invalid_estim!("Gamma cloglog reference backend produced non-finite values");
2538 }
2539 Ok(sval)
2540}
2541
2542pub(crate) fn cloglog_posterior_meanwith_deriv_controlled(
2543 ctx: &QuadratureContext,
2544 mu: f64,
2545 sigma: f64,
2546) -> IntegratedMeanDerivative {
2547 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
2592 return IntegratedMeanDerivative {
2593 mean: cloglog_mean_exact(mu),
2596 dmean_dmu: cloglog_mean_d1_exact(mu),
2599 mode: IntegratedExpectationMode::ExactClosedForm,
2600 };
2601 }
2602 if sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN {
2603 let (log_base, base_mode) = cloglog_log_survival_term_controlled(ctx, mu, sigma);
2609 let (log_shift, shift_mode) =
2610 cloglog_log_survival_term_controlled(ctx, mu + sigma * sigma, sigma);
2611 let mean = (-log_base.exp_m1()).clamp(0.0, 1.0);
2613 let dmean = cloglog_shift_identity_derivative_log(mu, sigma, log_shift);
2614 return IntegratedMeanDerivative {
2615 mean,
2616 dmean_dmu: dmean.max(0.0),
2617 mode: worse_integrated_expectation_mode(base_mode, shift_mode),
2618 };
2619 }
2620 let candidate = if sigma < CLOGLOG_SIGMA_TAYLOR_MAX {
2621 cloglog_small_sigma_taylor(mu, sigma)
2622 } else if let Some(out) = cloglog_extreme_asymptotic(mu, sigma) {
2623 out
2624 } else {
2625 let ((survival, mode), (shifted_survival, shifted_mode)) =
2626 cloglog_survival_pair_controlled(ctx, mu, sigma);
2627 if matches!(mode, IntegratedExpectationMode::QuadratureFallback)
2628 || matches!(shifted_mode, IntegratedExpectationMode::QuadratureFallback)
2629 {
2630 return cloglog_posterior_meanwith_deriv_quadrature(mu, sigma);
2631 }
2632 let mean = cloglog_mean_from_survival(survival);
2633 let dmean = cloglog_shift_identity_derivative(mu, sigma, shifted_survival);
2634 let mode = if matches!(mode, IntegratedExpectationMode::ControlledAsymptotic)
2635 || matches!(
2636 shifted_mode,
2637 IntegratedExpectationMode::ControlledAsymptotic
2638 ) {
2639 IntegratedExpectationMode::ControlledAsymptotic
2640 } else {
2641 mode
2642 };
2643 IntegratedMeanDerivative {
2644 mean,
2645 dmean_dmu: dmean.max(0.0),
2646 mode,
2647 }
2648 };
2649 if matches!(
2655 candidate.mode,
2656 IntegratedExpectationMode::ControlledAsymptotic
2657 ) && sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN
2658 {
2659 return candidate;
2660 }
2661 let ghq = cloglog_posterior_meanwith_deriv_quadrature(mu, sigma);
2662 if integrated_mean_derivative_drift_exceeds(&candidate, &ghq, 1e-6, 1e-4, 1e-7, 1e-3) {
2669 ghq
2670 } else {
2671 candidate
2672 }
2673}
2674
2675pub fn integrated_inverse_link_mean_and_derivative(
2676 quadctx: &QuadratureContext,
2677 link: LinkFunction,
2678 mu: f64,
2679 sigma: f64,
2680) -> Result<IntegratedMeanDerivative, EstimationError> {
2681 match link {
2710 LinkFunction::Log => {
2711 let (mean, saturated) = safe_expwith_saturation(mu + 0.5 * sigma * sigma);
2712 Ok(IntegratedMeanDerivative {
2713 mean,
2714 dmean_dmu: mean,
2715 mode: if saturated {
2716 IntegratedExpectationMode::ControlledAsymptotic
2717 } else {
2718 IntegratedExpectationMode::ExactClosedForm
2719 },
2720 })
2721 }
2722 LinkFunction::Probit => Ok(probit_posterior_meanwith_deriv_exact(mu, sigma)),
2723 LinkFunction::Logit => logit_posterior_meanwith_deriv_controlled(mu, sigma),
2724 LinkFunction::CLogLog => Ok(cloglog_posterior_meanwith_deriv_controlled(quadctx, mu, sigma)),
2725 LinkFunction::LogLog | LinkFunction::Cauchit => {
2726 let component = if matches!(link, LinkFunction::LogLog) {
2728 LinkComponent::LogLog
2729 } else {
2730 LinkComponent::Cauchit
2731 };
2732 let (mean, dmean_dmu, _, _) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
2733 component_point_jet(component, x)
2734 });
2735 Ok(IntegratedMeanDerivative {
2736 mean,
2737 dmean_dmu,
2738 mode: if sigma <= 1e-10 {
2739 IntegratedExpectationMode::ExactClosedForm
2740 } else {
2741 IntegratedExpectationMode::QuadratureFallback
2742 },
2743 })
2744 }
2745 LinkFunction::Sas => Err(EstimationError::InvalidInput(
2746 "state-less integrated SAS moments are unsupported; use SAS-aware prediction APIs with explicit (epsilon, log_delta)".to_string(),
2747 )),
2748 LinkFunction::BetaLogistic => Err(EstimationError::InvalidInput(
2749 "state-less integrated Beta-Logistic moments are unsupported; use link-aware prediction APIs with explicit (delta, epsilon)".to_string(),
2750 )),
2751 LinkFunction::Identity => Ok(IntegratedMeanDerivative {
2752 mean: mu,
2753 dmean_dmu: 1.0,
2754 mode: IntegratedExpectationMode::ExactClosedForm,
2755 }),
2756 }
2757}
2758
2759#[inline]
2760pub fn integrated_inverse_link_jet(
2761 quadctx: &QuadratureContext,
2762 link: LinkFunction,
2763 mu: f64,
2764 sigma: f64,
2765) -> Result<IntegratedInverseLinkJet, EstimationError> {
2766 match link {
2767 LinkFunction::Log => {
2768 let (mean, saturated) = safe_expwith_saturation(mu + 0.5 * sigma * sigma);
2769 Ok(IntegratedInverseLinkJet {
2770 mean,
2771 d1: mean,
2772 d2: mean,
2773 d3: mean,
2774 mode: if saturated {
2775 IntegratedExpectationMode::ControlledAsymptotic
2776 } else {
2777 IntegratedExpectationMode::ExactClosedForm
2778 },
2779 })
2780 }
2781 LinkFunction::Probit => Ok(integrated_probit_jet(mu, sigma)),
2782 LinkFunction::Logit => {
2783 if sigma > LOGIT_JET_GHQ_SIGMA_MAX {
2784 return logit_wide_sigma_jet(mu, sigma);
2788 }
2789 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
2794 component_point_jet(LinkComponent::Logit, x)
2795 });
2796 let mode = if sigma <= 1e-10 {
2797 IntegratedExpectationMode::ExactClosedForm
2798 } else {
2799 match logit_posterior_meanwith_deriv_controlled(mu, sigma) {
2803 Ok(scalar) => scalar.mode,
2804 Err(_) => IntegratedExpectationMode::QuadratureFallback,
2805 }
2806 };
2807 Ok(IntegratedInverseLinkJet {
2808 mean,
2809 d1: d1.max(0.0),
2810 d2,
2811 d3,
2812 mode,
2813 })
2814 }
2815 LinkFunction::CLogLog => {
2816 validate_latent_cloglog_inputs(mu, sigma)?;
2817 Ok(integrated_cloglog_inverse_link_jet_controlled(
2818 quadctx, mu, sigma,
2819 ))
2820 }
2821 LinkFunction::LogLog | LinkFunction::Cauchit => {
2822 let component = if matches!(link, LinkFunction::LogLog) {
2824 LinkComponent::LogLog
2825 } else {
2826 LinkComponent::Cauchit
2827 };
2828 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
2829 component_point_jet(component, x)
2830 });
2831 Ok(IntegratedInverseLinkJet {
2832 mean,
2833 d1,
2834 d2,
2835 d3,
2836 mode: if sigma <= 1e-10 {
2837 IntegratedExpectationMode::ExactClosedForm
2838 } else {
2839 IntegratedExpectationMode::QuadratureFallback
2840 },
2841 })
2842 }
2843 LinkFunction::Sas => Err(EstimationError::InvalidInput(
2844 "state-less integrated SAS jet is unsupported; use SAS-aware prediction APIs with explicit (epsilon, log_delta)".to_string(),
2845 )),
2846 LinkFunction::BetaLogistic => Err(EstimationError::InvalidInput(
2847 "state-less integrated Beta-Logistic jet is unsupported; use link-aware prediction APIs with explicit (delta, epsilon)".to_string(),
2848 )),
2849 LinkFunction::Identity => Ok(IntegratedInverseLinkJet {
2850 mean: mu,
2851 d1: 1.0,
2852 d2: 0.0,
2853 d3: 0.0,
2854 mode: IntegratedExpectationMode::ExactClosedForm,
2855 }),
2856 }
2857}
2858
2859#[inline]
2870fn logit_wide_sigma_jet(mu: f64, sigma: f64) -> Result<IntegratedInverseLinkJet, EstimationError> {
2871 let scalar = logit_posterior_meanwith_deriv_controlled(mu, sigma)?;
2872 let d2 = integrate_normal_adaptive(mu, sigma, |x| {
2873 component_point_jet(LinkComponent::Logit, x).2
2874 });
2875 let d3 = integrate_normal_adaptive(mu, sigma, |x| {
2876 component_point_jet(LinkComponent::Logit, x).3
2877 });
2878 Ok(IntegratedInverseLinkJet {
2879 mean: scalar.mean,
2880 d1: scalar.dmean_dmu.max(0.0),
2881 d2,
2882 d3,
2883 mode: scalar.mode,
2884 })
2885}
2886
2887#[inline]
2888fn sas_point_jet(x: f64, epsilon: f64, log_delta: f64) -> (f64, f64, f64, f64) {
2889 let jet = sas_inverse_link_jet(x, epsilon, log_delta)
2890 .expect("normal quadrature nodes must be finite");
2891 (jet.mu, jet.d1, jet.d2, jet.d3)
2892}
2893
2894#[inline]
2895fn beta_logistic_point_jet(x: f64, log_shape_center: f64, epsilon: f64) -> (f64, f64, f64, f64) {
2896 let jet = beta_logistic_inverse_link_jet(x, log_shape_center, epsilon);
2897 (jet.mu, jet.d1, jet.d2, jet.d3)
2898}
2899
2900#[inline]
2901fn worse_integrated_expectation_mode(
2902 lhs: IntegratedExpectationMode,
2903 rhs: IntegratedExpectationMode,
2904) -> IntegratedExpectationMode {
2905 if lhs.rank() >= rhs.rank() { lhs } else { rhs }
2906}
2907
2908#[inline]
2909fn integrated_scalar_drift_exceeds(
2910 candidate: f64,
2911 reference: f64,
2912 abs_tol: f64,
2913 rel_tol: f64,
2914) -> bool {
2915 if !(candidate.is_finite() && reference.is_finite()) {
2916 return true;
2917 }
2918 (candidate - reference).abs() > abs_tol.max(rel_tol * reference.abs().max(candidate.abs()))
2919}
2920
2921#[inline]
2922fn integrated_mean_derivative_drift_exceeds(
2923 candidate: &IntegratedMeanDerivative,
2924 reference: &IntegratedMeanDerivative,
2925 mean_abs_tol: f64,
2926 mean_rel_tol: f64,
2927 deriv_abs_tol: f64,
2928 deriv_rel_tol: f64,
2929) -> bool {
2930 integrated_scalar_drift_exceeds(candidate.mean, reference.mean, mean_abs_tol, mean_rel_tol)
2931 || integrated_scalar_drift_exceeds(
2932 candidate.dmean_dmu,
2933 reference.dmean_dmu,
2934 deriv_abs_tol,
2935 deriv_rel_tol,
2936 )
2937}
2938
2939#[inline]
2940fn component_point_jet(component: LinkComponent, x: f64) -> (f64, f64, f64, f64) {
2941 let jet = component_inverse_link_jet(component, x);
2944 (jet.mu, jet.d1, jet.d2, jet.d3)
2945}
2946
2947#[inline]
2948fn integrated_mixture_component_jet(
2949 ctx: &QuadratureContext,
2950 component: LinkComponent,
2951 mu: f64,
2952 sigma: f64,
2953) -> IntegratedInverseLinkJet {
2954 match component {
2959 LinkComponent::Logit => integrated_inverse_link_jet(ctx, LinkFunction::Logit, mu, sigma)
2960 .unwrap_or_else(|_| integrated_logit_jet_ghq(ctx, mu, sigma)),
2961 LinkComponent::Probit => integrated_probit_jet(mu, sigma),
2962 LinkComponent::CLogLog => integrated_cloglog_inverse_link_jet_controlled(ctx, mu, sigma),
2963 LinkComponent::LogLog | LinkComponent::Cauchit => {
2964 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
2965 component_point_jet(component, x)
2966 });
2967 IntegratedInverseLinkJet {
2968 mean,
2969 d1: d1.max(0.0),
2970 d2,
2971 d3,
2972 mode: if sigma <= 1e-10 {
2973 IntegratedExpectationMode::ExactClosedForm
2974 } else {
2975 IntegratedExpectationMode::QuadratureFallback
2976 },
2977 }
2978 }
2979 }
2980}
2981
2982#[inline]
2983fn integrated_mixture_jet(
2984 ctx: &QuadratureContext,
2985 mu: f64,
2986 sigma: f64,
2987 mixture_state: &MixtureLinkState,
2988) -> Result<IntegratedInverseLinkJet, EstimationError> {
2989 if mixture_state.components.is_empty() {
2994 crate::bail_invalid_estim!(
2995 "integrated mixture-link jet requires at least one blended component"
2996 );
2997 }
2998 if mixture_state.components.len() != mixture_state.pi.len() {
2999 crate::bail_invalid_estim!(
3000 "integrated mixture-link jet requires matching component and weight counts"
3001 );
3002 }
3003
3004 let mut mean = 0.0_f64;
3009 let mut d1 = 0.0_f64;
3010 let mut d2 = 0.0_f64;
3011 let mut d3 = 0.0_f64;
3012 let mut mode = IntegratedExpectationMode::ExactClosedForm;
3013 let mut saw_positive_weight = false;
3014
3015 for (&component, &weight) in mixture_state.components.iter().zip(mixture_state.pi.iter()) {
3016 if weight <= 0.0 {
3017 continue;
3018 }
3019 let jet = integrated_mixture_component_jet(ctx, component, mu, sigma);
3020 mean += weight * jet.mean;
3021 d1 += weight * jet.d1;
3022 d2 += weight * jet.d2;
3023 d3 += weight * jet.d3;
3024 if jet.mode.rank() > mode.rank() {
3025 mode = jet.mode;
3026 }
3027 saw_positive_weight = true;
3028 }
3029
3030 if !saw_positive_weight {
3031 crate::bail_invalid_estim!(
3032 "integrated mixture-link jet requires at least one positive component weight"
3033 .to_string(),
3034 );
3035 }
3036
3037 Ok(IntegratedInverseLinkJet {
3038 mean,
3039 d1: d1.max(0.0),
3040 d2,
3041 d3,
3042 mode,
3043 })
3044}
3045
3046#[inline]
3047fn integrated_sas_jet_ghq(
3048 ctx: &QuadratureContext,
3049 mu: f64,
3050 sigma: f64,
3051 sas_state: &SasLinkState,
3052) -> IntegratedInverseLinkJet {
3053 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3054 sas_point_jet(x, sas_state.epsilon, sas_state.log_delta)
3055 });
3056 IntegratedInverseLinkJet {
3057 mean,
3058 d1: d1.max(0.0),
3059 d2,
3060 d3,
3061 mode: if sigma <= 1e-10 {
3062 IntegratedExpectationMode::ExactClosedForm
3063 } else {
3064 IntegratedExpectationMode::QuadratureFallback
3065 },
3066 }
3067}
3068
3069#[inline]
3070fn integrated_beta_logistic_jet_ghq(
3071 ctx: &QuadratureContext,
3072 mu: f64,
3073 sigma: f64,
3074 beta_state: &SasLinkState,
3075) -> IntegratedInverseLinkJet {
3076 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3077 beta_logistic_point_jet(x, beta_state.log_delta, beta_state.epsilon)
3078 });
3079 IntegratedInverseLinkJet {
3080 mean,
3081 d1: d1.max(0.0),
3082 d2,
3083 d3,
3084 mode: if sigma <= 1e-10 {
3085 IntegratedExpectationMode::ExactClosedForm
3086 } else {
3087 IntegratedExpectationMode::QuadratureFallback
3088 },
3089 }
3090}
3091
3092#[inline]
3094pub fn integrated_inverse_link_jetwith_state(
3095 quadctx: &QuadratureContext,
3096 link: LinkFunction,
3097 mu: f64,
3098 sigma: f64,
3099 mixture_link_state: Option<&MixtureLinkState>,
3100 sas_link_state: Option<&SasLinkState>,
3101) -> Result<IntegratedInverseLinkJet, EstimationError> {
3102 if let Some(state) = mixture_link_state {
3103 return integrated_mixture_jet(quadctx, mu, sigma, state);
3104 }
3105 if matches!(link, LinkFunction::Sas) {
3106 let sas = sas_link_state.ok_or_else(|| {
3107 EstimationError::InvalidInput(
3108 "state-less integrated SAS jet is unsupported; explicit SasLinkState is required"
3109 .to_string(),
3110 )
3111 })?;
3112 return Ok(integrated_sas_jet_ghq(quadctx, mu, sigma, sas));
3113 }
3114 if matches!(link, LinkFunction::BetaLogistic) {
3115 let state = sas_link_state.ok_or_else(|| {
3116 EstimationError::InvalidInput(
3117 "state-less integrated Beta-Logistic jet is unsupported; explicit link state is required"
3118 .to_string(),
3119 )
3120 })?;
3121 return Ok(integrated_beta_logistic_jet_ghq(quadctx, mu, sigma, state));
3122 }
3123 integrated_inverse_link_jet(quadctx, link, mu, sigma)
3124}
3125
3126#[inline]
3136pub fn integrated_family_moments_jet(
3137 quadctx: &QuadratureContext,
3138 likelihood: &GlmLikelihoodSpec,
3139 eta: f64,
3140 se_eta: f64,
3141) -> Result<IntegratedMomentsJet, EstimationError> {
3142 const PROB_EPS: f64 = 1e-12;
3143 if !(eta.is_finite() && (-700.0..=700.0).contains(&eta)) {
3144 crate::bail_invalid_estim!(
3145 "integrated moments eta must be finite and within [-700, 700]; got {eta}"
3146 );
3147 }
3148 let e = eta;
3149 let se = se_eta.max(0.0);
3150 let resolved_scale = likelihood
3154 .resolved_scale()
3155 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3156 let spec = &likelihood.spec;
3157 let mixture_link_state: Option<&MixtureLinkState> = spec.link.mixture_state();
3158 let sas_link_state: Option<&SasLinkState> = spec.link.sas_state();
3159 match &spec.response {
3160 ResponseFamily::Binomial => match &spec.link {
3161 InverseLink::Standard(StandardLink::Logit) => {
3162 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3163 let mean = jet.mean;
3164 Ok(IntegratedMomentsJet {
3165 mean,
3166 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3167 d1: jet.d1,
3168 d2: jet.d2,
3169 d3: jet.d3,
3170 mode: jet.mode,
3171 })
3172 }
3173 InverseLink::Standard(StandardLink::Probit) => {
3174 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Probit, e, se)?;
3175 let mean = jet.mean;
3176 Ok(IntegratedMomentsJet {
3177 mean,
3178 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3179 d1: jet.d1,
3180 d2: jet.d2,
3181 d3: jet.d3,
3182 mode: jet.mode,
3183 })
3184 }
3185 InverseLink::Standard(StandardLink::CLogLog) => {
3186 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::CLogLog, e, se)?;
3187 let mean = jet.mean;
3188 Ok(IntegratedMomentsJet {
3189 mean,
3190 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3191 d1: jet.d1,
3192 d2: jet.d2,
3193 d3: jet.d3,
3194 mode: jet.mode,
3195 })
3196 }
3197 InverseLink::LatentCLogLog(_) => Err(EstimationError::InvalidInput(
3198 "Binomial+LatentCLogLog integrated moments require an explicit latent cloglog inverse-link state"
3199 .to_string(),
3200 )),
3201 InverseLink::Sas(_) => {
3202 let jet = integrated_inverse_link_jetwith_state(
3203 quadctx,
3204 LinkFunction::Sas,
3205 e,
3206 se,
3207 mixture_link_state,
3208 sas_link_state,
3209 )?;
3210 let mean = jet.mean;
3211 Ok(IntegratedMomentsJet {
3212 mean,
3213 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3214 d1: jet.d1,
3215 d2: jet.d2,
3216 d3: jet.d3,
3217 mode: jet.mode,
3218 })
3219 }
3220 InverseLink::BetaLogistic(_) => {
3221 let jet = integrated_inverse_link_jetwith_state(
3222 quadctx,
3223 LinkFunction::BetaLogistic,
3224 e,
3225 se,
3226 mixture_link_state,
3227 sas_link_state,
3228 )?;
3229 let mean = jet.mean;
3230 Ok(IntegratedMomentsJet {
3231 mean,
3232 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3233 d1: jet.d1,
3234 d2: jet.d2,
3235 d3: jet.d3,
3236 mode: jet.mode,
3237 })
3238 }
3239 InverseLink::Mixture(state) => {
3240 let jet = integrated_mixture_jet(quadctx, e, se, &state)?;
3241 let mean = jet.mean;
3242 Ok(IntegratedMomentsJet {
3243 mean,
3244 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3245 d1: jet.d1,
3246 d2: jet.d2,
3247 d3: jet.d3,
3248 mode: jet.mode,
3249 })
3250 }
3251 InverseLink::Standard(other) => Err(EstimationError::InvalidInput(format!(
3252 "Binomial response paired with unsupported standard link {other:?} for integrated moments"
3253 ))),
3254 },
3255 ResponseFamily::Gaussian => {
3256 let variance = resolved_scale
3257 .gaussian_phi()
3258 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3259 Ok(IntegratedMomentsJet {
3260 mean: e,
3261 variance,
3262 d1: 1.0,
3263 d2: 0.0,
3264 d3: 0.0,
3265 mode: IntegratedExpectationMode::ExactClosedForm,
3266 })
3267 }
3268 ResponseFamily::RoystonParmar => {
3269 let jet = integrated_inverse_link_jetwith_state(
3270 quadctx,
3271 LinkFunction::CLogLog,
3272 e,
3273 se,
3274 mixture_link_state,
3275 sas_link_state,
3276 )?;
3277 let mean = (1.0 - jet.mean).clamp(0.0, 1.0);
3278 Ok(IntegratedMomentsJet {
3279 mean,
3280 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3281 d1: -jet.d1,
3282 d2: -jet.d2,
3283 d3: -jet.d3,
3284 mode: jet.mode,
3285 })
3286 }
3287 ResponseFamily::Beta { .. } => {
3288 let precision = resolved_scale
3289 .beta_precision()
3290 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3291 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3292 let mean = jet.mean.clamp(PROB_EPS, 1.0 - PROB_EPS);
3293 Ok(IntegratedMomentsJet {
3294 mean,
3295 variance: (mean * (1.0 - mean) / (1.0 + precision)).max(PROB_EPS),
3296 d1: jet.d1,
3297 d2: jet.d2,
3298 d3: jet.d3,
3299 mode: jet.mode,
3300 })
3301 }
3302 ResponseFamily::Poisson
3303 | ResponseFamily::Tweedie { .. }
3304 | ResponseFamily::NegativeBinomial { .. }
3305 | ResponseFamily::Gamma => {
3306 let s2 = se * se;
3311 let (mean, saturated) = safe_expwith_saturation(e + 0.5 * s2);
3312 let variance = match &spec.response {
3323 ResponseFamily::Poisson => mean,
3324 ResponseFamily::Tweedie { p } => {
3325 let phi = resolved_scale
3326 .tweedie_phi()
3327 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3328 phi * mean.powf(*p)
3329 }
3330 ResponseFamily::NegativeBinomial { .. } => {
3331 let theta = resolved_scale.negative_binomial_theta().map_err(|error| {
3332 EstimationError::InvalidInput(error.to_string())
3333 })?;
3334 mean + mean * mean / theta
3335 }
3336 ResponseFamily::Gamma => {
3337 let phi = resolved_scale
3338 .gamma_phi()
3339 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3340 phi * mean * mean
3341 }
3342 other => {
3346 return Err(EstimationError::InvalidInput(format!(
3347 "integrated log-normal moments reached unexpected family {other:?}"
3348 )));
3349 }
3350 };
3351 if !(variance.is_finite() && variance >= 0.0) {
3352 return Err(EstimationError::InvalidInput(format!(
3353 "integrated {} variance is not representable: {variance:?}",
3354 spec.response.name()
3355 )));
3356 }
3357 Ok(IntegratedMomentsJet {
3358 mean,
3359 variance,
3360 d1: mean,
3361 d2: mean,
3362 d3: mean,
3363 mode: if saturated {
3364 IntegratedExpectationMode::ControlledAsymptotic
3365 } else {
3366 IntegratedExpectationMode::ExactClosedForm
3367 },
3368 })
3369 }
3370 }
3371}
3372
3373pub fn logit_posterior_meanwith_deriv_batch(
3376 ctx: &QuadratureContext,
3377 eta: &ndarray::Array1<f64>,
3378 se_eta: &ndarray::Array1<f64>,
3379) -> Result<(ndarray::Array1<f64>, ndarray::Array1<f64>), EstimationError> {
3380 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3381 let n = eta.len();
3382 let pairs: Result<Vec<(f64, f64)>, _> = (0..n)
3384 .into_par_iter()
3385 .map(|i| {
3386 let integrated = integrated_inverse_link_mean_and_derivative(
3387 ctx,
3388 LinkFunction::Logit,
3389 eta[i],
3390 se_eta[i],
3391 )?;
3392 Ok::<_, EstimationError>((integrated.mean, integrated.dmean_dmu))
3393 })
3394 .collect();
3395 let pairs = pairs?;
3396 let mut mu = ndarray::Array1::<f64>::zeros(n);
3397 let mut dmu = ndarray::Array1::<f64>::zeros(n);
3398 for (i, (m, d)) in pairs.into_iter().enumerate() {
3399 mu[i] = m;
3400 dmu[i] = d;
3401 }
3402
3403 Ok((mu, dmu))
3404}
3405
3406pub fn logit_posterior_mean_batch(
3410 ctx: &QuadratureContext,
3411 eta: &ndarray::Array1<f64>,
3412 se_eta: &ndarray::Array1<f64>,
3413) -> Result<ndarray::Array1<f64>, EstimationError> {
3414 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3415 let n = eta.len();
3416 let values: Result<Vec<f64>, EstimationError> = (0..n)
3417 .into_par_iter()
3418 .map(|i| {
3419 integrated_inverse_link_mean_and_derivative(ctx, LinkFunction::Logit, eta[i], se_eta[i])
3420 .map(|integrated| integrated.mean)
3421 })
3422 .collect();
3423 Ok(ndarray::Array1::from_vec(values?))
3424}
3425
3426pub trait GhqValue: Sized {
3427 fn zero() -> Self;
3428 fn addweighted(&mut self, weight: f64, value: Self);
3429 fn scale(self, factor: f64) -> Self;
3430}
3431
3432impl GhqValue for f64 {
3433 #[inline]
3434 fn zero() -> Self {
3435 0.0
3436 }
3437
3438 #[inline]
3439 fn addweighted(&mut self, weight: f64, value: Self) {
3440 *self += weight * value;
3441 }
3442
3443 #[inline]
3444 fn scale(self, factor: f64) -> Self {
3445 self * factor
3446 }
3447}
3448
3449impl GhqValue for (f64, f64) {
3450 #[inline]
3451 fn zero() -> Self {
3452 (0.0, 0.0)
3453 }
3454
3455 #[inline]
3456 fn addweighted(&mut self, weight: f64, value: Self) {
3457 self.0 += weight * value.0;
3458 self.1 += weight * value.1;
3459 }
3460
3461 #[inline]
3462 fn scale(self, factor: f64) -> Self {
3463 (self.0 * factor, self.1 * factor)
3464 }
3465}
3466
3467impl GhqValue for (f64, f64, f64, f64) {
3468 #[inline]
3469 fn zero() -> Self {
3470 (0.0, 0.0, 0.0, 0.0)
3471 }
3472
3473 #[inline]
3474 fn addweighted(&mut self, weight: f64, value: Self) {
3475 self.0 += weight * value.0;
3476 self.1 += weight * value.1;
3477 self.2 += weight * value.2;
3478 self.3 += weight * value.3;
3479 }
3480
3481 #[inline]
3482 fn scale(self, factor: f64) -> Self {
3483 (
3484 self.0 * factor,
3485 self.1 * factor,
3486 self.2 * factor,
3487 self.3 * factor,
3488 )
3489 }
3490}
3491
3492impl GhqValue for (f64, f64, f64, f64, f64, f64) {
3493 #[inline]
3494 fn zero() -> Self {
3495 (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
3496 }
3497
3498 #[inline]
3499 fn addweighted(&mut self, weight: f64, value: Self) {
3500 self.0 += weight * value.0;
3501 self.1 += weight * value.1;
3502 self.2 += weight * value.2;
3503 self.3 += weight * value.3;
3504 self.4 += weight * value.4;
3505 self.5 += weight * value.5;
3506 }
3507
3508 #[inline]
3509 fn scale(self, factor: f64) -> Self {
3510 (
3511 self.0 * factor,
3512 self.1 * factor,
3513 self.2 * factor,
3514 self.3 * factor,
3515 self.4 * factor,
3516 self.5 * factor,
3517 )
3518 }
3519}
3520
3521#[inline]
3522fn integrate_normal_ghq_adaptive<F, R>(ctx: &QuadratureContext, eta: f64, se_eta: f64, f: F) -> R
3523where
3524 F: Fn(f64) -> R,
3525 R: GhqValue,
3526{
3527 if se_eta < 1e-10 {
3528 return f(eta);
3529 }
3530 let n = adaptive_point_count_from_sd(se_eta.abs());
3531 with_gh_nodesweights(ctx, n, |nodes, weights| {
3532 let scale = SQRT_2 * se_eta;
3533 let mut sum = R::zero();
3534 for i in 0..n {
3535 sum.addweighted(weights[i], f(eta + scale * nodes[i]));
3536 }
3537 sum.scale(1.0 / std::f64::consts::PI.sqrt())
3538 })
3539}
3540
3541#[inline]
3542fn integrated_probit_jet(mu: f64, sigma: f64) -> IntegratedInverseLinkJet {
3543 let s = sigma.hypot(1.0);
3549 let z = mu / s;
3550 let mean = gam_math::probability::normal_cdf(z);
3551 let pdf = gam_math::probability::normal_pdf(z);
3552 if pdf == 0.0 {
3553 return IntegratedInverseLinkJet {
3554 mean,
3555 d1: 0.0,
3556 d2: 0.0,
3557 d3: 0.0,
3558 mode: IntegratedExpectationMode::ExactClosedForm,
3559 };
3560 }
3561 IntegratedInverseLinkJet {
3562 mean,
3563 d1: pdf / s,
3564 d2: -z * pdf / (s * s),
3565 d3: (z * z - 1.0) * pdf / (s * s * s),
3566 mode: IntegratedExpectationMode::ExactClosedForm,
3567 }
3568}
3569
3570#[inline]
3571fn integrated_logit_jet_ghq(
3572 ctx: &QuadratureContext,
3573 mu: f64,
3574 sigma: f64,
3575) -> IntegratedInverseLinkJet {
3576 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3577 component_point_jet(LinkComponent::Logit, x)
3578 });
3579 IntegratedInverseLinkJet {
3580 mean,
3581 d1: d1.max(0.0),
3582 d2,
3583 d3,
3584 mode: if sigma <= 1e-10 {
3585 IntegratedExpectationMode::ExactClosedForm
3586 } else {
3587 IntegratedExpectationMode::QuadratureFallback
3588 },
3589 }
3590}
3591
3592#[inline]
3593fn cloglog_inverse_link_controlled_values(
3594 ctx: &QuadratureContext,
3595 mu: f64,
3596 sigma: f64,
3597 max_order: usize,
3598) -> ([f64; 6], IntegratedExpectationMode) {
3599 assert!(max_order <= 5);
3600 if sigma <= 1e-10 {
3601 let (mean, d1, d2, d3, d4, d5) = cloglog_point_jet5(mu);
3602 return (
3603 [mean, d1, d2, d3, d4, d5],
3604 IntegratedExpectationMode::ExactClosedForm,
3605 );
3606 }
3607
3608 let (k, log_k0, mode) = latent_cloglog_kernel_terms(ctx, mu, sigma, max_order);
3609 let mut values = [0.0; 6];
3610 values[0] = if log_k0.is_finite() {
3611 -log_k0.exp_m1()
3612 } else {
3613 1.0
3614 };
3615 values[1] = k[1].max(0.0);
3616 if sigma > CLOGLOG_JET_MOMENT_SIGMA_MAX {
3617 if max_order >= 2 {
3618 values[2] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).2);
3619 }
3620 if max_order >= 3 {
3621 values[3] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).3);
3622 }
3623 if max_order >= 4 {
3624 values[4] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).4);
3625 }
3626 if max_order >= 5 {
3627 values[5] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).5);
3628 }
3629 return (
3630 values,
3631 worse_integrated_expectation_mode(mode, IntegratedExpectationMode::QuadratureFallback),
3632 );
3633 }
3634 if max_order >= 2 {
3635 values[2] = k[1] - k[2];
3636 }
3637 if max_order >= 3 {
3638 values[3] = k[1] - 3.0 * k[2] + k[3];
3639 }
3640 if max_order >= 4 {
3641 values[4] = k[1] - 7.0 * k[2] + 6.0 * k[3] - k[4];
3642 }
3643 if max_order >= 5 {
3644 values[5] = k[1] - 15.0 * k[2] + 25.0 * k[3] - 10.0 * k[4] + k[5];
3645 }
3646 (values, mode)
3647}
3648
3649#[inline]
3650pub(crate) fn latent_cloglog_inverse_link_jet5_controlled(
3651 ctx: &QuadratureContext,
3652 mu: f64,
3653 sigma: f64,
3654) -> IntegratedInverseLinkJet5 {
3655 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 5);
3656 IntegratedInverseLinkJet5 {
3657 mean: values[0],
3658 d1: values[1],
3659 d2: values[2],
3660 d3: values[3],
3661 d4: values[4],
3662 d5: values[5],
3663 mode,
3664 }
3665}
3666
3667#[derive(Clone, Copy, Debug)]
3677pub struct LatentCLogLogJet5 {
3678 pub mean: f64,
3679 pub d1: f64,
3680 pub d2: f64,
3681 pub d3: f64,
3682 pub d4: f64,
3683 pub d5: f64,
3684 pub mode: IntegratedExpectationMode,
3685}
3686
3687pub fn latent_cloglog_jet5(
3688 quadctx: &QuadratureContext,
3689 eta: f64,
3690 sigma: f64,
3691) -> Result<LatentCLogLogJet5, EstimationError> {
3692 validate_latent_cloglog_inputs(eta, sigma)?;
3693 let jet = latent_cloglog_inverse_link_jet5_controlled(quadctx, eta, sigma);
3699 Ok(LatentCLogLogJet5 {
3700 mean: jet.mean,
3701 d1: jet.d1,
3702 d2: jet.d2,
3703 d3: jet.d3,
3704 d4: jet.d4,
3705 d5: jet.d5,
3706 mode: jet.mode,
3707 })
3708}
3709
3710#[inline]
3711pub fn latent_cloglog_inverse_link_jet(
3712 quadctx: &QuadratureContext,
3713 eta: f64,
3714 sigma: f64,
3715) -> Result<IntegratedInverseLinkJet, EstimationError> {
3716 let jet = latent_cloglog_jet5(quadctx, eta, sigma)?;
3717 Ok(IntegratedInverseLinkJet {
3718 mean: jet.mean,
3719 d1: jet.d1,
3720 d2: jet.d2,
3721 d3: jet.d3,
3722 mode: jet.mode,
3723 })
3724}
3725
3726#[inline]
3727fn integrated_cloglog_inverse_link_jet_controlled(
3728 ctx: &QuadratureContext,
3729 mu: f64,
3730 sigma: f64,
3731) -> IntegratedInverseLinkJet {
3732 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 3);
3733 IntegratedInverseLinkJet {
3734 mean: values[0],
3735 d1: values[1],
3736 d2: values[2],
3737 d3: values[3],
3738 mode,
3739 }
3740}
3741
3742#[inline]
3743fn latent_cloglog_kernel_terms(
3744 ctx: &QuadratureContext,
3745 mu: f64,
3746 sigma: f64,
3747 max_order: usize,
3748) -> ([f64; 6], f64, IntegratedExpectationMode) {
3749 let sigma2 = sigma * sigma;
3750 let mut k = [0.0; 6];
3751 let mut log_k0 = f64::NEG_INFINITY;
3752 let mut mode = IntegratedExpectationMode::ExactClosedForm;
3753
3754 for (order, out) in k.iter_mut().enumerate().take(max_order + 1) {
3755 let kf = order as f64;
3756 let shifted_mu = mu + kf * sigma2;
3757 let (log_survival, term_mode) =
3765 cloglog_log_survival_term_controlled(ctx, shifted_mu, sigma);
3766 mode = worse_integrated_expectation_mode(mode, term_mode);
3767
3768 let log_value = kf * mu + 0.5 * kf * kf * sigma2 + log_survival;
3769 if order == 0 {
3770 log_k0 = log_value;
3771 }
3772 if !log_value.is_finite() {
3773 *out = 0.0;
3774 continue;
3775 }
3776 let upper = if order == 0 {
3777 1.0
3778 } else {
3779 let k_over_e = kf / std::f64::consts::E;
3780 k_over_e.powf(kf)
3781 };
3782 *out = safe_exp(log_value).clamp(0.0, upper);
3783 }
3784
3785 (k, log_k0, mode)
3786}
3787
3788#[inline]
3789pub fn normal_expectation_1d_adaptive<F>(
3790 ctx: &QuadratureContext,
3791 eta: f64,
3792 se_eta: f64,
3793 f: F,
3794) -> f64
3795where
3796 F: Fn(f64) -> f64,
3797{
3798 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
3799}
3800
3801#[inline]
3802pub fn normal_expectation_1d_adaptive_pair<F>(
3803 ctx: &QuadratureContext,
3804 eta: f64,
3805 se_eta: f64,
3806 f: F,
3807) -> (f64, f64)
3808where
3809 F: Fn(f64) -> (f64, f64),
3810{
3811 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
3812}
3813
3814fn adaptive_point_count_from_sd(max_sd: f64) -> usize {
3815 if max_sd.is_finite() && max_sd > 2.5 {
3825 51
3826 } else if max_sd.is_finite() && max_sd > 0.5 {
3827 31
3828 } else if max_sd.is_finite() && max_sd > 0.35 {
3829 21
3830 } else if max_sd.is_finite() && max_sd > 0.1 {
3831 15
3832 } else {
3833 7
3834 }
3835}
3836
3837#[inline]
3838fn with_gh_nodesweights<R>(
3839 ctx: &QuadratureContext,
3840 n: usize,
3841 f: impl FnOnce(&[f64], &[f64]) -> R,
3842) -> R {
3843 if n == 7 {
3844 let gh = ctx.gauss_hermite();
3845 f(&gh.nodes, &gh.weights)
3846 } else {
3847 let gh = ctx.gauss_hermite_n(n);
3848 f(&gh.nodes, &gh.weights)
3849 }
3850}
3851
3852#[inline]
3862fn cholesky_static<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
3863 let mut l = [[0.0_f64; D]; D];
3864 for i in 0..D {
3865 for j in 0..=i {
3866 let mut sum = cov[i][j];
3867 for k in 0..j {
3868 sum -= l[i][k] * l[j][k];
3869 }
3870 if i == j {
3871 if !sum.is_finite() || sum <= 0.0 {
3872 return None;
3873 }
3874 l[i][j] = sum.sqrt();
3875 } else {
3876 l[i][j] = sum / l[j][j];
3877 }
3878 }
3879 }
3880 Some(l)
3881}
3882
3883#[inline]
3886fn cholesky_static_with_jitter<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
3887 if D == 0 {
3888 return None;
3889 }
3890 for retry in 0..8 {
3891 let jitter = if retry == 0 {
3892 0.0
3893 } else {
3894 1e-12 * 10f64.powi(retry - 1)
3895 };
3896 if jitter == 0.0 {
3897 if let Some(l) = cholesky_static::<D>(cov) {
3898 return Some(l);
3899 }
3900 } else {
3901 let mut base = *cov;
3902 for i in 0..D {
3903 base[i][i] = cov[i][i] + jitter;
3904 }
3905 if let Some(l) = cholesky_static::<D>(&base) {
3906 return Some(l);
3907 }
3908 }
3909 }
3910 None
3911}
3912
3913#[inline]
3914fn adaptive_point_countwith_cap(max_sd: f64, max_n: usize) -> usize {
3915 adaptive_point_count_from_sd(max_sd).min(max_n)
3916}
3917
3918#[inline]
3919fn ghq_nd_integrate_try<const D: usize, F, R, E>(
3920 ctx: &QuadratureContext,
3921 mu: [f64; D],
3922 cov: [[f64; D]; D],
3923 max_n: usize,
3924 f: F,
3925) -> Result<Option<R>, E>
3926where
3927 F: Fn([f64; D]) -> Result<R, E>,
3928 R: GhqValue,
3929{
3930 let mut maxvar = 0.0_f64;
3931 for (i, row) in cov.iter().enumerate() {
3932 maxvar = maxvar.max(row[i]).max(0.0);
3933 }
3934 let n = adaptive_point_countwith_cap(maxvar.sqrt(), max_n);
3935
3936 let mut cov_arr = cov;
3941 for i in 0..D {
3942 cov_arr[i][i] = cov_arr[i][i].max(0.0);
3943 }
3944 let Some(l) = cholesky_static_with_jitter::<D>(&cov_arr) else {
3945 return Ok(None);
3946 };
3947 let norm = 1.0 / std::f64::consts::PI.powf(0.5 * D as f64);
3948
3949 with_gh_nodesweights(ctx, n, |nodes, weights| {
3950 let mut acc = R::zero();
3951 let mut idx = [0usize; D];
3952 loop {
3953 let mut z = [0.0_f64; D];
3954 let mut weight = 1.0_f64;
3955 for d in 0..D {
3956 z[d] = SQRT_2 * nodes[idx[d]];
3957 weight *= weights[idx[d]];
3958 }
3959
3960 let mut x = mu;
3961 for row in 0..D {
3962 let mut dot = 0.0_f64;
3963 for (col, zc) in z.iter().enumerate().take(row + 1) {
3964 dot += l[row][col] * *zc;
3965 }
3966 x[row] += dot;
3967 }
3968 acc.addweighted(weight, f(x)?);
3969
3970 let mut carry = true;
3971 for d in (0..D).rev() {
3972 idx[d] += 1;
3973 if idx[d] < n {
3974 carry = false;
3975 break;
3976 }
3977 idx[d] = 0;
3978 }
3979 if carry {
3980 break;
3981 }
3982 }
3983 Ok(Some(acc.scale(norm)))
3984 })
3985}
3986
3987#[inline]
3988fn ghq_nd_integrate<const D: usize, F, R>(
3989 ctx: &QuadratureContext,
3990 mu: [f64; D],
3991 cov: [[f64; D]; D],
3992 max_n: usize,
3993 f: F,
3994) -> Option<R>
3995where
3996 F: Fn([f64; D]) -> R,
3997 R: GhqValue,
3998{
3999 match ghq_nd_integrate_try::<D, _, R, Infallible>(ctx, mu, cov, max_n, |x| Ok(f(x))) {
4000 Ok(v) => v,
4001 Err(e) => match e {},
4002 }
4003}
4004
4005#[inline]
4006fn ghq_nd_integrate_result<const D: usize, F, R, E>(
4007 ctx: &QuadratureContext,
4008 mu: [f64; D],
4009 cov: [[f64; D]; D],
4010 max_n: usize,
4011 f: F,
4012) -> Result<Option<R>, E>
4013where
4014 F: Fn([f64; D]) -> Result<R, E>,
4015 R: GhqValue,
4016{
4017 ghq_nd_integrate_try::<D, _, R, E>(ctx, mu, cov, max_n, f)
4018}
4019
4020pub fn normal_expectation_nd_adaptive<const D: usize, F>(
4022 ctx: &QuadratureContext,
4023 mu: [f64; D],
4024 cov: [[f64; D]; D],
4025 max_n: usize,
4026 f: F,
4027) -> f64
4028where
4029 F: Fn([f64; D]) -> f64,
4030{
4031 match ghq_nd_integrate::<D, _, f64>(ctx, mu, cov, max_n, &f) {
4032 Some(v) => v,
4033 None => f(mu),
4034 }
4035}
4036
4037pub fn normal_expectation_nd_adaptive_result<const D: usize, F, R, E>(
4039 ctx: &QuadratureContext,
4040 mu: [f64; D],
4041 cov: [[f64; D]; D],
4042 max_n: usize,
4043 f: F,
4044) -> Result<R, E>
4045where
4046 F: Fn([f64; D]) -> Result<R, E>,
4047 R: GhqValue,
4048{
4049 match ghq_nd_integrate_result::<D, _, R, E>(ctx, mu, cov, max_n, &f)? {
4050 Some(v) => Ok(v),
4051 None => f(mu),
4052 }
4053}
4054
4055pub fn normal_expectation_2d_adaptive_result<F, E>(
4057 ctx: &QuadratureContext,
4058 mu: [f64; 2],
4059 cov: [[f64; 2]; 2],
4060 f: F,
4061) -> Result<f64, E>
4062where
4063 F: Fn(f64, f64) -> Result<f64, E>,
4064{
4065 normal_expectation_nd_adaptive_result::<2, _, _, E>(ctx, mu, cov, 21, |x| f(x[0], x[1]))
4066}
4067
4068pub fn normal_expectation_3d_adaptive<F>(
4070 ctx: &QuadratureContext,
4071 mu: [f64; 3],
4072 cov: [[f64; 3]; 3],
4073 f: F,
4074) -> f64
4075where
4076 F: Fn(f64, f64, f64) -> f64,
4077{
4078 normal_expectation_nd_adaptive::<3, _>(ctx, mu, cov, 15, |x| f(x[0], x[1], x[2]))
4080}
4081
4082#[inline]
4101pub fn probit_posterior_mean(eta: f64, se_eta: f64) -> f64 {
4102 if se_eta < 1e-10 {
4103 return gam_math::probability::normal_cdf(eta);
4104 }
4105 let denom = (1.0 + se_eta * se_eta).sqrt();
4106 gam_math::probability::normal_cdf(eta / denom)
4107}
4108
4109#[inline]
4110pub fn logit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4111 let (m1, m2) = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4112 let p = sigmoid(x);
4113 (p, p * p)
4114 });
4115 let m1 = m1.clamp(0.0, 1.0);
4116 let m2 = m2.clamp(0.0, 1.0);
4117 (m1, (m2 - m1 * m1).max(0.0))
4118}
4119
4120#[inline]
4121pub fn probit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4122 let m1 = probit_posterior_mean(eta, se_eta);
4123 let m2 = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4124 let p = gam_math::probability::normal_cdf(x);
4125 p * p
4126 })
4127 .clamp(0.0, 1.0);
4128 (m1, (m2 - m1 * m1).max(0.0))
4129}
4130
4131#[inline]
4132pub fn cloglog_posterior_meanvariance(
4133 ctx: &QuadratureContext,
4134 eta: f64,
4135 se_eta: f64,
4136) -> (f64, f64) {
4137 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4157 return (cloglog_mean_exact(eta), 0.0);
4158 }
4159 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4160 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4161 let mean = cloglog_mean_from_survival(survival);
4162 let variance = (survival_sq - survival * survival).max(0.0);
4163 (mean, variance)
4164}
4165
4166#[inline]
4200pub fn cloglog_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4201 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4205 return cloglog_mean_exact(eta);
4206 }
4207 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4208 cloglog_mean_from_survival(survival)
4209}
4210
4211#[inline]
4225pub fn survival_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4226 cloglog_survival_term_controlled(ctx, eta, se_eta)
4227 .0
4228 .clamp(0.0, 1.0)
4229}
4230
4231#[inline]
4232pub fn survival_posterior_meanvariance(
4233 ctx: &QuadratureContext,
4234 eta: f64,
4235 se_eta: f64,
4236) -> (f64, f64) {
4237 let (m1, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4238 let (m2, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4239 (m1.clamp(0.0, 1.0), (m2 - m1 * m1).max(0.0))
4240}
4241
4242pub fn logit_posterior_mean_exact(mu: f64, sigma: f64) -> f64 {
4318 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 0.0 {
4319 return sigmoid(mu);
4320 }
4321 if sigma < LOGIT_SIGMA_DEGENERATE {
4322 return sigmoid(mu);
4325 }
4326
4327 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI; let sqrt2_sigma = SQRT_2 * sigma;
4329 let coeff = (2.0_f64 * std::f64::consts::PI).sqrt() / sigma; let c = -mu / sqrt2_sigma; let beta = std::f64::consts::PI / sqrt2_sigma; let r2 = FADDEEVA_ASYMPTOTIC_RADIUS * FADDEEVA_ASYMPTOTIC_RADIUS;
4333
4334 let mut corr = 0.0_f64;
4340 let mut n = 1usize;
4341 let tail_start = loop {
4342 let b = (2.0 * (n as f64) - 1.0) * beta;
4343 let abs_xi2 = c * c + b * b;
4344 if abs_xi2 > r2 && n >= FADDEEVA_TAIL_MIN_INDEX {
4345 break n;
4346 }
4347 let xi = Complex { re: c, im: b };
4348 let d = if abs_xi2 > r2 {
4349 inv_sqrt_pi * faddeeva_asymptotic_a(xi).re
4351 } else {
4352 faddeeva_upper_halfplane(xi).im - inv_sqrt_pi * c / abs_xi2
4353 };
4354 corr += d;
4355 n += 1;
4356 };
4357
4358 corr += faddeeva_pole_series_em_tail(c, beta, tail_start, inv_sqrt_pi);
4359
4360 sigmoid(mu) - coeff * corr
4361}
4362
4363const FADDEEVA_TAIL_MIN_INDEX: usize = 48;
4367const FADDEEVA_ASYMPTOTIC_RADIUS: f64 = 7.0;
4370const FADDEEVA_ASYMPTOTIC_TERMS: usize = 14;
4373
4374fn faddeeva_asymptotic_a(xi: Complex) -> Complex {
4378 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4379 let inv2 = complexmul(inv, inv);
4380 let mut xp = complexmul(inv2, inv); let mut cm = 0.5_f64; let mut s = Complex::default();
4383 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4384 s = complex_add(
4385 s,
4386 Complex {
4387 re: cm * xp.re,
4388 im: cm * xp.im,
4389 },
4390 );
4391 cm *= (2.0 * (m as f64) + 1.0) / 2.0; xp = complexmul(xp, inv2);
4393 }
4394 s
4395}
4396
4397fn faddeeva_pole_series_em_tail(c: f64, beta: f64, tail_start: usize, inv_sqrt_pi: f64) -> f64 {
4406 let b_a = (2.0 * (tail_start as f64) - 1.0) * beta;
4407 let xi = Complex { re: c, im: b_a };
4408 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4409 let inv2 = complexmul(inv, inv);
4410 let two_i_beta = Complex {
4412 re: 0.0,
4413 im: 2.0 * beta,
4414 };
4415
4416 let mut s = Complex::default(); let mut a_acc = Complex::default(); let mut fp_inner = Complex::default(); let mut x2m = inv2; let mut x2m1 = complexmul(inv2, inv); let mut x2m2 = complexmul(inv2, inv2); let mut cm = 0.5_f64;
4424 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4425 let mf = m as f64;
4426 let inv_4ibm = Complex {
4428 re: 0.0,
4429 im: -1.0 / (4.0 * beta * mf),
4430 };
4431 s = complex_add(
4432 s,
4433 complexmul(
4434 Complex {
4435 re: cm * x2m.re,
4436 im: cm * x2m.im,
4437 },
4438 inv_4ibm,
4439 ),
4440 );
4441 a_acc = complex_add(
4442 a_acc,
4443 Complex {
4444 re: cm * x2m1.re,
4445 im: cm * x2m1.im,
4446 },
4447 );
4448 let fc = cm * (-(2.0 * mf + 1.0));
4449 fp_inner = complex_add(
4450 fp_inner,
4451 Complex {
4452 re: fc * x2m2.re,
4453 im: fc * x2m2.im,
4454 },
4455 );
4456 cm *= (2.0 * mf + 1.0) / 2.0;
4457 x2m = complexmul(x2m, inv2);
4458 x2m1 = complexmul(x2m1, inv2);
4459 x2m2 = complexmul(x2m2, inv2);
4460 }
4461
4462 s = complex_add(
4464 s,
4465 Complex {
4466 re: 0.5 * a_acc.re,
4467 im: 0.5 * a_acc.im,
4468 },
4469 );
4470 let fprime = complexmul(two_i_beta, fp_inner);
4472 s = complex_add(
4473 s,
4474 Complex {
4475 re: -fprime.re / 12.0,
4476 im: -fprime.im / 12.0,
4477 },
4478 );
4479
4480 inv_sqrt_pi * s.re
4483}
4484
4485fn faddeeva_upper_halfplane(z: Complex) -> Complex {
4498 let (l, coeffs) = faddeeva_weideman_coeffs();
4499 let iz = Complex {
4500 re: -z.im,
4501 im: z.re,
4502 }; let l_minus = Complex {
4504 re: l - iz.re,
4505 im: -iz.im,
4506 }; let l_plus = Complex {
4508 re: l + iz.re,
4509 im: iz.im,
4510 }; let zz = complex_div(l_plus, l_minus); let mut p = Complex {
4514 re: coeffs[0],
4515 im: 0.0,
4516 };
4517 for &c in &coeffs[1..] {
4518 p = complex_add(complexmul(p, zz), Complex { re: c, im: 0.0 });
4519 }
4520 let l_minus_sq = complexmul(l_minus, l_minus);
4521 let term1 = complex_div(
4522 Complex {
4523 re: 2.0 * p.re,
4524 im: 2.0 * p.im,
4525 },
4526 l_minus_sq,
4527 );
4528 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI;
4529 let term2 = complex_div(
4530 Complex {
4531 re: inv_sqrt_pi,
4532 im: 0.0,
4533 },
4534 l_minus,
4535 );
4536 complex_add(term1, term2)
4537}
4538
4539const FADDEEVA_WEIDEMAN_N: usize = 44;
4542
4543fn faddeeva_weideman_coeffs() -> &'static (f64, [f64; FADDEEVA_WEIDEMAN_N]) {
4549 static CACHE: OnceLock<(f64, [f64; FADDEEVA_WEIDEMAN_N])> = OnceLock::new();
4550 CACHE.get_or_init(|| {
4551 let n = FADDEEVA_WEIDEMAN_N;
4552 let l = (n as f64 / SQRT_2).sqrt();
4553 let m = 2 * n;
4554 let m2 = 2 * m; let mut f = vec![0.0_f64; m2];
4558 for (idx, fi) in f.iter_mut().enumerate().skip(1) {
4559 let k = (idx as isize - 1) - (m as isize - 1);
4560 let theta = (k as f64) * std::f64::consts::PI / (m as f64);
4561 let t = l * (0.5 * theta).tan();
4562 *fi = (-t * t).exp() * (l * l + t * t);
4563 }
4564 let half = m2 / 2;
4567 let mut coeffs = [0.0_f64; FADDEEVA_WEIDEMAN_N];
4568 for j in 1..=n {
4569 let mut acc = 0.0_f64;
4570 for (p, _) in f.iter().enumerate() {
4571 let fp = f[(p + half) % m2];
4572 if fp != 0.0 {
4573 acc += fp
4574 * (-2.0 * std::f64::consts::PI * (j as f64) * (p as f64) / (m2 as f64))
4575 .cos();
4576 }
4577 }
4578 coeffs[n - j] = acc / (m2 as f64);
4580 }
4581 (l, coeffs)
4582 })
4583}
4584
4585#[inline]
4587fn sigmoid(x: f64) -> f64 {
4588 let x_clamped = x.clamp(-QUADRATURE_EXP_LOG_MAX, QUADRATURE_EXP_LOG_MAX);
4589 1.0 / (1.0 + f64::exp(-x_clamped))
4590}
4591
4592#[derive(Clone, Copy, Debug)]
4608pub struct CLogLogConvolutionDerivatives {
4609 pub l: f64,
4611
4612 pub l_mu: f64,
4614 pub l_sigma: f64,
4615
4616 pub l_mumu: f64,
4618 pub l_musigma: f64,
4619 pub l_sigmasigma: f64,
4620
4621 pub l_mumumu: f64,
4623 pub l_mumusigma: f64,
4624 pub l_musigmasigma: f64,
4625 pub l_sigmasigmasigma: f64,
4626
4627 pub l_mumumumu: f64,
4629 pub l_mumumusigma: f64,
4630 pub l_mumusigmasigma: f64,
4631 pub l_musigmasigmasigma: f64,
4632 pub l_sigmasigmasigmasigma: f64,
4633}
4634
4635#[inline]
4636pub(crate) fn cloglog_point_jet5(t: f64) -> (f64, f64, f64, f64, f64, f64) {
4637 if t.is_nan() {
4638 return (f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN);
4639 }
4640 let et = safe_exp(t);
4641
4642 (
4643 -(-et).exp_m1(),
4644 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0]),
4645 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -1.0]),
4646 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -3.0, 1.0]),
4647 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -7.0, 6.0, -1.0]),
4648 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -15.0, 25.0, -10.0, 1.0]),
4649 )
4650}
4651
4652#[inline]
4664fn cloglog_g_derivatives(t: f64) -> (f64, f64, f64, f64, f64) {
4665 let (g, g1, g2, g3, g4, _) = cloglog_point_jet5(t);
4666 (g, g1, g2, g3, g4)
4667}
4668
4669pub fn cloglog_ghq_value(ctx: &QuadratureContext, mu: f64, sigma: f64, n_nodes: usize) -> f64 {
4677 if sigma.abs() < 1e-14 {
4678 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
4679 return g.clamp(0.0, 1.0);
4680 }
4681 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
4682
4683 let inv_sig2 = 1.0 / (sigma * sigma);
4715 let mut eta_hat = mu;
4716 let mut converged = false;
4717 for _ in 0..100 {
4718 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
4719 if !(g > 0.0) || !g1.is_finite() || !g2.is_finite() {
4720 break;
4721 }
4722 let r = g1 / g;
4723 let lp = r - (eta_hat - mu) * inv_sig2;
4724 let lpp = g2 / g - r * r - inv_sig2;
4725 if !lpp.is_finite() || lpp >= 0.0 {
4726 break;
4727 }
4728 let step = lp / lpp;
4729 eta_hat -= step;
4730 if step.abs() <= 1e-13 * (1.0 + eta_hat.abs()) {
4731 converged = true;
4732 break;
4733 }
4734 }
4735
4736 let tau = if converged {
4740 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
4741 if g > 0.0 {
4742 let r = g1 / g;
4743 let lpp = g2 / g - r * r - inv_sig2;
4744 let tau2 = -1.0 / lpp;
4745 if tau2.is_finite() && tau2 > 0.0 {
4746 Some(tau2.sqrt())
4747 } else {
4748 None
4749 }
4750 } else {
4751 None
4752 }
4753 } else {
4754 None
4755 };
4756
4757 let eval_at = |n: usize| -> f64 {
4759 match tau {
4760 Some(tau) => {
4761 let pref = tau * inv_sqrt_pi / sigma;
4762 with_gh_nodesweights(ctx, n, |nodes, weights| {
4763 let mut sum = 0.0_f64;
4764 for i in 0..nodes.len() {
4765 let t = nodes[i];
4766 let eta_i = eta_hat + SQRT_2 * tau * t;
4767 let (g, _, _, _, _, _) = cloglog_point_jet5(eta_i);
4768 let dev = eta_i - mu;
4769 sum += weights[i] * (t * t - 0.5 * dev * dev * inv_sig2).exp() * g;
4770 }
4771 (pref * sum).clamp(0.0, 1.0)
4772 })
4773 }
4774 None => {
4775 let scale = SQRT_2 * sigma;
4776 with_gh_nodesweights(ctx, n, |nodes, weights| {
4777 let mut sum = 0.0_f64;
4778 for i in 0..nodes.len() {
4779 let t = mu + scale * nodes[i];
4780 let (g, _, _, _, _) = cloglog_g_derivatives(t);
4781 sum += weights[i] * g;
4782 }
4783 (sum * inv_sqrt_pi).clamp(0.0, 1.0)
4784 })
4785 }
4786 }
4787 };
4788
4789 const CLOGLOG_GHQ_ORDER_LADDER: [usize; 5] = [7, 15, 21, 31, 51];
4794 const CLOGLOG_GHQ_CONV_TOL: f64 = 1e-10;
4795 let floor = n_nodes.min(*CLOGLOG_GHQ_ORDER_LADDER.last().unwrap());
4796 let mut prev: Option<f64> = None;
4797 let mut result = 0.0_f64;
4798 for &n in CLOGLOG_GHQ_ORDER_LADDER.iter().filter(|&&n| n >= floor) {
4799 let cur = eval_at(n);
4800 result = cur;
4801 if let Some(p) = prev
4802 && (cur - p).abs() < CLOGLOG_GHQ_CONV_TOL
4803 {
4804 break;
4805 }
4806 prev = Some(cur);
4807 }
4808 result
4809}
4810
4811pub fn cloglog_ghq_derivatives(
4822 ctx: &QuadratureContext,
4823 mu: f64,
4824 sigma: f64,
4825 n_nodes: usize,
4826) -> CLogLogConvolutionDerivatives {
4827 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
4828
4829 if sigma.abs() < 1e-14 {
4836 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
4837 return CLogLogConvolutionDerivatives {
4838 l: g,
4839 l_mu: g1,
4840 l_sigma: 0.0,
4841 l_mumu: g2,
4842 l_musigma: 0.0,
4843 l_sigmasigma: g2,
4844 l_mumumu: g3,
4845 l_mumusigma: 0.0,
4846 l_musigmasigma: g3,
4847 l_sigmasigmasigma: 0.0,
4848 l_mumumumu: g4,
4849 l_mumumusigma: 0.0,
4850 l_mumusigmasigma: g4,
4851 l_musigmasigmasigma: 0.0,
4852 l_sigmasigmasigmasigma: 3.0 * g4,
4853 };
4854 }
4855
4856 let scale = SQRT_2 * sigma;
4857 let sqrt2 = SQRT_2;
4858
4859 with_gh_nodesweights(ctx, n_nodes, |nodes, weights| {
4860 let mut s = [[0.0_f64; 5]; 5];
4872
4873 for i in 0..nodes.len() {
4874 let x = nodes[i];
4875 let t = mu + scale * x;
4876 let (g0, g1, g2, g3, g4) = cloglog_g_derivatives(t);
4877 let w = weights[i];
4878
4879 let x2 = x * x;
4881 let x3 = x2 * x;
4882 let x4 = x3 * x;
4883
4884 s[0][0] += w * g0;
4886
4887 s[1][0] += w * g1;
4889 s[1][1] += w * x * g1;
4890
4891 s[2][0] += w * g2;
4893 s[2][1] += w * x * g2;
4894 s[2][2] += w * x2 * g2;
4895
4896 s[3][0] += w * g3;
4898 s[3][1] += w * x * g3;
4899 s[3][2] += w * x2 * g3;
4900 s[3][3] += w * x3 * g3;
4901
4902 s[4][0] += w * g4;
4904 s[4][1] += w * x * g4;
4905 s[4][2] += w * x2 * g4;
4906 s[4][3] += w * x3 * g4;
4907 s[4][4] += w * x4 * g4;
4908 }
4909
4910 let sqrt2_1 = sqrt2;
4913 let sqrt2_2 = 2.0; let sqrt2_3 = 2.0 * sqrt2; let sqrt2_4 = 4.0; CLogLogConvolutionDerivatives {
4918 l: inv_sqrt_pi * s[0][0],
4920
4921 l_mu: inv_sqrt_pi * s[1][0],
4923 l_sigma: inv_sqrt_pi * sqrt2_1 * s[1][1],
4924
4925 l_mumu: inv_sqrt_pi * s[2][0],
4927 l_musigma: inv_sqrt_pi * sqrt2_1 * s[2][1],
4928 l_sigmasigma: inv_sqrt_pi * sqrt2_2 * s[2][2],
4929
4930 l_mumumu: inv_sqrt_pi * s[3][0],
4932 l_mumusigma: inv_sqrt_pi * sqrt2_1 * s[3][1],
4933 l_musigmasigma: inv_sqrt_pi * sqrt2_2 * s[3][2],
4934 l_sigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[3][3],
4935
4936 l_mumumumu: inv_sqrt_pi * s[4][0],
4938 l_mumumusigma: inv_sqrt_pi * sqrt2_1 * s[4][1],
4939 l_mumusigmasigma: inv_sqrt_pi * sqrt2_2 * s[4][2],
4940 l_musigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[4][3],
4941 l_sigmasigmasigmasigma: inv_sqrt_pi * sqrt2_4 * s[4][4],
4942 }
4943 })
4944}
4945
4946pub fn cloglog_ghq_derivatives_adaptive(
4952 ctx: &QuadratureContext,
4953 mu: f64,
4954 sigma: f64,
4955) -> CLogLogConvolutionDerivatives {
4956 let n = adaptive_point_count_from_sd(sigma.abs());
4957 cloglog_ghq_derivatives(ctx, mu, sigma, n)
4958}
4959
4960#[cfg(test)]
4961mod tests {
4962 use super::*;
4963 use approx::assert_relative_eq;
4964 use gam_problem::LikelihoodScaleMetadata;
4965 use gam_spec::LikelihoodSpec;
4966
4967 #[test]
4974 fn log_half_erfc_stable_matches_high_precision_reference() {
4975 let refs: &[(f64, f64)] = &[
4976 (-3.0, -1.1045309498499094e-5),
4977 (-1.5, -0.017092677825984745),
4978 (-0.5, -0.27410803278438573),
4979 (0.0, -0.69314718055994531),
4980 (0.7, -1.8257336940742865),
4981 (2.0, -6.0580884451765829),
4982 (5.0, -27.89403672609738),
4983 (12.0, -147.75386135854695),
4984 ];
4985 for &(u, reference) in refs {
4986 let got = log_half_erfc_stable(u);
4987 let rel = (got - reference).abs() / reference.abs().max(1.0e-6);
4988 assert!(
4989 rel < 1.0e-12,
4990 "log_half_erfc_stable({u}) = {got:.17e}, reference {reference:.17e}, \
4991 rel {rel:.3e} >= 1e-12"
4992 );
4993 }
4994 }
4995
4996 pub(crate) fn cloglog_posterior_meanwith_deriv_gamma_reference(
4997 mu: f64,
4998 sigma: f64,
4999 ) -> Result<IntegratedMeanDerivative, EstimationError> {
5000 let survival = cloglog_survival_gamma_reference(mu, sigma)?;
5003 let shifted_survival = cloglog_survival_gamma_reference(mu + sigma * sigma, sigma)?;
5004 let mean = cloglog_mean_from_survival(survival);
5005 let dmean = cloglog_shift_identity_derivative(mu, sigma, shifted_survival);
5006 if !(mean.is_finite() && dmean.is_finite()) {
5007 crate::bail_invalid_estim!(
5008 "Gamma cloglog reference backend produced non-finite values"
5009 );
5010 }
5011 Ok(IntegratedMeanDerivative {
5012 mean,
5013 dmean_dmu: dmean.max(0.0),
5014 mode: IntegratedExpectationMode::ExactSpecialFunction,
5015 })
5016 }
5017
5018 fn even_moment_exp_neg_x2(power: usize) -> f64 {
5019 assert!(power.is_multiple_of(2));
5020 let m = power / 2;
5021 let mut odd_double_factorial = 1.0_f64;
5022 for k in 0..m {
5023 odd_double_factorial *= (2 * k + 1) as f64;
5024 }
5025 odd_double_factorial * std::f64::consts::PI.sqrt() / 2.0_f64.powi(m as i32)
5026 }
5027
5028 fn normal_pdf(z: f64) -> f64 {
5029 (-(z * z) * 0.5).exp() / (2.0 * std::f64::consts::PI).sqrt()
5030 }
5031
5032 fn high_res_sigmoid_integral(eta: f64, se: f64) -> f64 {
5033 let a = -12.0_f64;
5035 let b = 12.0_f64;
5036 let n = 20_000usize; let h = (b - a) / n as f64;
5038
5039 let integrand = |z: f64| -> f64 { sigmoid(eta + se * z) * normal_pdf(z) };
5040
5041 let mut sum = integrand(a) + integrand(b);
5042 for i in 1..n {
5043 let x = a + (i as f64) * h;
5044 if i % 2 == 0 {
5045 sum += 2.0 * integrand(x);
5046 } else {
5047 sum += 4.0 * integrand(x);
5048 }
5049 }
5050 sum * h / 3.0
5051 }
5052
5053 #[test]
5054 fn test_computed_nodes_symmetric() {
5055 let ctx = QuadratureContext::new();
5057 let gh = ctx.gauss_hermite();
5058 for i in 0..N_POINTS / 2 {
5059 let j = N_POINTS - 1 - i;
5060 assert_relative_eq!(gh.nodes[i], -gh.nodes[j], epsilon = 1e-12);
5061 }
5062 assert_relative_eq!(gh.nodes[N_POINTS / 2], 0.0, epsilon = 1e-12);
5064 }
5065
5066 #[test]
5067 fn test_computedweights_symmetric() {
5068 let ctx = QuadratureContext::new();
5070 let gh = ctx.gauss_hermite();
5071 for i in 0..N_POINTS / 2 {
5072 let j = N_POINTS - 1 - i;
5073 assert_relative_eq!(gh.weights[i], gh.weights[j], epsilon = 1e-12);
5074 }
5075 }
5076
5077 #[test]
5078 fn testweights_sum_to_sqrt_pi() {
5079 let ctx = QuadratureContext::new();
5081 let gh = ctx.gauss_hermite();
5082 let sum: f64 = gh.weights.iter().sum();
5083 assert_relative_eq!(sum, std::f64::consts::PI.sqrt(), epsilon = 1e-10);
5084 }
5085
5086 #[test]
5087 fn test_clenshaw_curtisweights_are_symmetric_and_integrate_constants() {
5088 let rule = compute_clenshaw_curtis_n(33);
5089 let m = rule.weights.len() - 1;
5090 for j in 0..=m / 2 {
5091 assert_relative_eq!(rule.nodes[j], -rule.nodes[m - j], epsilon = 1e-14);
5092 assert_relative_eq!(rule.weights[j], rule.weights[m - j], epsilon = 1e-14);
5093 }
5094 let sum: f64 = rule.weights.iter().sum();
5095 assert_relative_eq!(sum, 2.0, epsilon = 1e-14, max_relative = 1e-14);
5096 }
5097
5098 #[test]
5099 fn test_cc_preference_prefers_moderate_central_case() {
5100 assert!(cloglog_should_prefer_cc(-0.2, 0.8, CLOGLOG_CC_TOL));
5101 }
5102
5103 #[test]
5104 fn test_cc_preference_prefers_moderately_large_case() {
5105 assert!(cloglog_should_prefer_cc(0.0, 2.0, CLOGLOG_CC_TOL));
5106 }
5107
5108 #[test]
5109 fn test_cc_preference_rejects_broad_case() {
5110 assert!(!cloglog_should_prefer_cc(0.0, 5.0, CLOGLOG_CC_TOL));
5111 }
5112
5113 #[test]
5114 fn testwilkinson_shift_finitewhen_d_iszero() {
5115 let shift = wilkinson_shift(0.0, 0.0, 1.25);
5118 assert!(shift.is_finite());
5119 assert_relative_eq!(shift, -1.25, epsilon = 1e-14);
5120 }
5121
5122 #[test]
5123 fn test_matches_abramowitz_stegun_7_point_gauss_hermite_constants() {
5124 let known_nodes = [
5128 -2.651_961_356_835_233_4,
5129 -1.673_551_628_767_471_4,
5130 -0.816_287_882_858_964_7,
5131 0.0,
5132 0.816_287_882_858_964_7,
5133 1.673_551_628_767_471_4,
5134 2.651_961_356_835_233_4,
5135 ];
5136 let knownweights = [
5137 0.000_971_781_245_099_519_1,
5138 0.054_515_582_819_127_03,
5139 0.425_607_252_610_127_8,
5140 0.810_264_617_556_807_3,
5141 0.425_607_252_610_127_8,
5142 0.054_515_582_819_127_03,
5143 0.000_971_781_245_099_519_1,
5144 ];
5145
5146 let ctx = QuadratureContext::new();
5147 let gh = ctx.gauss_hermite();
5148 for i in 0..N_POINTS {
5149 assert_relative_eq!(gh.nodes[i], known_nodes[i], epsilon = 1e-12);
5150 assert_relative_eq!(gh.weights[i], knownweights[i], epsilon = 1e-12);
5151 }
5152 }
5153
5154 #[test]
5155 fn test_gauss_hermite_weight_assembly_uses_eigenvector_rows() {
5156 let mut diag = [0.0_f64; N_POINTS];
5157 let mut off_diag = [0.0_f64; N_POINTS - 1];
5158 for (i, od) in off_diag.iter_mut().enumerate() {
5159 *od = (((i + 1) as f64) / 2.0).sqrt();
5160 }
5161 let (nodes, eigenvectors) = symmetric_tridiagonal_eigen(&mut diag, &mut off_diag);
5162 let mu0 = std::f64::consts::PI.sqrt();
5163 let mut row_pairs: Vec<(f64, f64)> = (0..N_POINTS)
5164 .map(|i| (nodes[i], mu0 * eigenvectors[i][0] * eigenvectors[i][0]))
5165 .collect();
5166 let mut column_pairs: Vec<(f64, f64)> = (0..N_POINTS)
5167 .map(|i| (nodes[i], mu0 * eigenvectors[0][i] * eigenvectors[0][i]))
5168 .collect();
5169 row_pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
5170 column_pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
5171
5172 let knownweights = [
5173 0.000_971_781_245_099_519_1,
5174 0.054_515_582_819_127_03,
5175 0.425_607_252_610_127_8,
5176 0.810_264_617_556_807_3,
5177 0.425_607_252_610_127_8,
5178 0.054_515_582_819_127_03,
5179 0.000_971_781_245_099_519_1,
5180 ];
5181
5182 for i in 0..N_POINTS {
5183 assert_relative_eq!(row_pairs[i].1, knownweights[i], epsilon = 1e-12);
5184 }
5185 let column_error: f64 = column_pairs
5186 .iter()
5187 .zip(knownweights.iter())
5188 .map(|(actual, expected)| (actual.1 - expected).abs())
5189 .sum();
5190 assert!(
5191 column_error > 1.0,
5192 "column-oriented eigenvector indexing unexpectedly matched A&S weights"
5193 );
5194 }
5195
5196 #[test]
5197 fn testzero_se_returns_mode() {
5198 let eta = 1.5;
5200 let se = 0.0;
5201 let ctx = QuadratureContext::new();
5202 let mean = logit_posterior_mean(&ctx, eta, se);
5203 let mode = sigmoid(eta);
5204 assert_relative_eq!(mean, mode, epsilon = 1e-10);
5205 }
5206
5207 #[test]
5208 fn test_symmetric_atzero() {
5209 let eta = 0.0;
5211 let se = 1.0;
5212 let ctx = QuadratureContext::new();
5213 let mean = logit_posterior_mean(&ctx, eta, se);
5214 assert_relative_eq!(mean, 0.5, epsilon = 0.01);
5216 }
5217
5218 #[test]
5219 fn test_shrinkage_at_extremes() {
5220 let eta = 3.0; let se = 1.0;
5223 let ctx = QuadratureContext::new();
5224 let mean = logit_posterior_mean(&ctx, eta, se);
5225 let mode = sigmoid(eta);
5226
5227 assert!(mean < mode, "Expected mean {} < mode {}", mean, mode);
5229 assert!(mean > 0.8, "Mean {} should still be high", mean);
5231 }
5232
5233 #[test]
5234 fn test_matches_monte_carlo() {
5235 let eta = 2.0;
5237 let se = 0.8;
5238
5239 let ctx = QuadratureContext::new();
5240 let quad_mean = logit_posterior_mean(&ctx, eta, se);
5241
5242 let n_samples = 100_000;
5244 let mut mc_sum = 0.0;
5245 let mut rng_state = 12345u64; for _ in 0..n_samples {
5247 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5249 let u1 = ((rng_state as f64) / (u64::MAX as f64)).max(1e-10); rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5251 let u2 = (rng_state as f64) / (u64::MAX as f64);
5252 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
5253 let eta_sample = eta + se * z;
5254 mc_sum += sigmoid(eta_sample);
5255 }
5256 let mc_mean = mc_sum / (n_samples as f64);
5257
5258 assert_relative_eq!(quad_mean, mc_mean, epsilon = 0.01);
5260 }
5261
5262 #[test]
5263 fn test_quadrature_integrates_x_squared() {
5264 let ctx = QuadratureContext::new();
5267 let gh = ctx.gauss_hermite();
5268 let mut sum = 0.0;
5269 for i in 0..N_POINTS {
5270 sum += gh.weights[i] * gh.nodes[i] * gh.nodes[i];
5271 }
5272 let expected = std::f64::consts::PI.sqrt() / 2.0;
5273 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5274 }
5275
5276 #[test]
5277 fn test_quadrature_integrates_x_fourth() {
5278 let ctx = QuadratureContext::new();
5281 let gh = ctx.gauss_hermite();
5282 let mut sum = 0.0;
5283 for i in 0..N_POINTS {
5284 let x = gh.nodes[i];
5285 sum += gh.weights[i] * x * x * x * x;
5286 }
5287 let expected = 3.0 * std::f64::consts::PI.sqrt() / 4.0;
5288 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5289 }
5290
5291 #[test]
5292 fn test_moment_exactness_up_to_degree_13() {
5293 let ctx = QuadratureContext::new();
5294 let gh = ctx.gauss_hermite();
5295
5296 for degree in 0..=13usize {
5297 let approx: f64 = (0..N_POINTS)
5298 .map(|i| gh.weights[i] * gh.nodes[i].powi(degree as i32))
5299 .sum();
5300
5301 let expected = if degree % 2 == 1 {
5302 0.0
5303 } else {
5304 even_moment_exp_neg_x2(degree)
5305 };
5306
5307 let err = (approx - expected).abs();
5308 let rel_scale = approx.abs().max(expected.abs()).max(1.0);
5309 assert!(
5310 err <= 1e-10 || err / rel_scale <= 1e-10,
5311 "degree={} approx={} expected={} abs_err={}",
5312 degree,
5313 approx,
5314 expected,
5315 err
5316 );
5317 }
5318 }
5319
5320 #[test]
5321 fn test_integrated_sigmoid_matches_high_res_integral_random_pairs() {
5322 let ctx = QuadratureContext::new();
5323 let mut rng_state = 0x4d595df4d0f33173u64;
5324
5325 for _ in 0..20 {
5326 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5327 let u_eta = (rng_state as f64) / (u64::MAX as f64);
5328 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5329 let u_se = (rng_state as f64) / (u64::MAX as f64);
5330
5331 let eta = -6.0 + 12.0 * u_eta;
5332 let se = 0.02 + 1.5 * u_se;
5333
5334 let ghq = logit_posterior_mean(&ctx, eta, se);
5335 let numeric = high_res_sigmoid_integral(eta, se);
5336 assert_relative_eq!(ghq, numeric, epsilon = 2e-3);
5337 }
5338 }
5339
5340 #[test]
5341 fn test_logit_posterior_derivative_remains_positive_in_positive_tail() {
5342 let eta = 20.0;
5343 let se = 0.0;
5344 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5345 .expect("logit posterior mean derivative should evaluate");
5346 assert!(dmu > 0.0);
5347 assert!(
5348 dmu < 1e-6,
5349 "positive-tail derivative should stay tiny but nonzero, got {dmu}"
5350 );
5351 }
5352
5353 #[test]
5354 fn test_logit_posterior_derivative_matches_central_difference() {
5355 let ctx = QuadratureContext::new();
5356 let eta = 1.7;
5357 let se = 0.9;
5358 let h = 1e-5;
5359
5360 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5361 .expect("logit posterior mean derivative should evaluate");
5362 let mu_plus = logit_posterior_mean(&ctx, eta + h, se);
5363 let mu_minus = logit_posterior_mean(&ctx, eta - h, se);
5364 let dmufd = (mu_plus - mu_minus) / (2.0 * h);
5365
5366 assert_eq!(dmu.signum(), dmufd.signum());
5367 assert_relative_eq!(dmu, dmufd, epsilon = 5e-6, max_relative = 2e-4);
5368 }
5369
5370 fn dense_sigmoid_normal_mean(mu: f64, sigma: f64) -> f64 {
5376 let a = -18.0_f64;
5377 let b = 18.0_f64;
5378 let n = 400_000usize; let h = (b - a) / n as f64;
5380 let integrand = |z: f64| -> f64 { sigmoid(mu + sigma * z) * normal_pdf(z) };
5381 let mut sum = integrand(a) + integrand(b);
5382 for i in 1..n {
5383 let z = a + (i as f64) * h;
5384 sum += if i % 2 == 0 { 2.0 } else { 4.0 } * integrand(z);
5385 }
5386 sum * h / 3.0
5387 }
5388
5389 #[test]
5390 fn test_logit_posterior_mean_exact_symmetry_identity() {
5391 let cases = [
5394 (-3.0, 0.5),
5395 (-1.2, 1.7),
5396 (0.0, 2.2),
5397 (2.3, 0.8),
5398 (3.0, 0.05),
5399 ];
5400 for (mu, sigma) in cases {
5401 let p = logit_posterior_mean_exact(mu, sigma);
5402 let q = logit_posterior_mean_exact(-mu, sigma);
5403 assert!(
5404 (p + q - 1.0).abs() < 1e-12,
5405 "symmetry broken at mu={mu} sigma={sigma}: p+q-1 = {:.3e}",
5406 p + q - 1.0
5407 );
5408 }
5409 }
5410
5411 #[test]
5412 fn test_logit_posterior_mean_exact_matches_high_res_integral() {
5413 let cases = [
5417 (-2.0, 0.4),
5418 (-0.7, 1.1),
5419 (0.8, 0.9),
5420 (2.4, 1.7),
5421 (3.0, 0.05),
5422 (3.0, 0.5),
5423 (-2.0, 2.0),
5424 (5.0, 3.0),
5425 ];
5426 for (mu, sigma) in cases {
5427 let exact = logit_posterior_mean_exact(mu, sigma);
5428 let numeric = dense_sigmoid_normal_mean(mu, sigma);
5429 assert!(
5430 (exact - numeric).abs() < 1e-10,
5431 "oracle ≠ dense reference at mu={mu} sigma={sigma}: \
5432 exact={exact:.13} ref={numeric:.13} err={:.3e}",
5433 (exact - numeric).abs()
5434 );
5435 }
5436 }
5437
5438 #[test]
5445 fn test_logit_posterior_mean_exact_no_truncation_bias_1459() {
5446 let table = [
5449 (1.0, 0.02),
5450 (1.0, 0.05),
5451 (1.0, 0.5),
5452 (1.0, 2.0),
5453 (3.0, 0.02),
5454 (3.0, 0.05),
5455 (3.0, 0.5),
5456 (3.0, 2.0),
5457 (-2.0, 0.02),
5458 (-2.0, 0.05),
5459 (-2.0, 0.5),
5460 (-2.0, 2.0),
5461 ];
5462 for (mu, sigma) in table {
5463 let exact = logit_posterior_mean_exact(mu, sigma);
5464 let reference = dense_sigmoid_normal_mean(mu, sigma);
5465 let err = (exact - reference).abs();
5466 assert!(
5467 err < 1e-10,
5468 "#1459 truncation bias resurfaced at mu={mu} sigma={sigma}: \
5469 err={err:.3e} (pre-fix bias here was ~{:.2e})",
5470 mu.abs() / (2.0 * std::f64::consts::PI.powi(2) * 4096.0)
5471 );
5472 }
5473
5474 let mu = 3.0;
5479 let errs: Vec<f64> = [0.05, 0.5, 2.0]
5480 .iter()
5481 .map(|&s| logit_posterior_mean_exact(mu, s) - dense_sigmoid_normal_mean(mu, s))
5482 .collect();
5483 for e in &errs {
5484 assert!(
5485 e.abs() < 1e-10,
5486 "residual {e:.3e} at mu=3 — old σ-independent plateau was 3.71e-5"
5487 );
5488 }
5489 }
5490
5491 #[test]
5499 fn test_faddeeva_weideman_matches_known_values() {
5500 let w0 = faddeeva_upper_halfplane(Complex { re: 0.0, im: 0.0 });
5502 assert!(
5503 (w0.re - 1.0).abs() < 1e-13 && w0.im.abs() < 1e-13,
5504 "w(0)={w0:?}"
5505 );
5506 let on_axis = [
5508 (0.1, 0.8964569799691268),
5509 (0.5, 0.6156903441929258),
5510 (1.0, 0.427583576155807),
5511 (2.0, 0.2553956763105058),
5512 (5.0, 0.11070463773306861),
5513 (9.0, 0.06230772403777468),
5514 ];
5515 for (y, want) in on_axis {
5516 let w = faddeeva_upper_halfplane(Complex { re: 0.0, im: y });
5517 assert!(
5518 (w.re - want).abs() < 1e-13 && w.im.abs() < 1e-13,
5519 "w(i·{y}): got {w:?}, want re={want}, err={:.2e}",
5520 (w.re - want).abs()
5521 );
5522 }
5523 let off_axis = [
5525 ((0.7, 1.3), (0.31327301971562715, 0.12443489420104513)),
5526 ((-1.5, 0.8), (0.21066359024766423, -0.27001624496296617)),
5527 ((3.0, 0.4), (0.030278754646989155, 0.1957320888774461)),
5528 ];
5529 for ((re, im), (wre, wim)) in off_axis {
5530 let w = faddeeva_upper_halfplane(Complex { re, im });
5531 assert!(
5532 (w.re - wre).abs() < 1e-13 && (w.im - wim).abs() < 1e-13,
5533 "w({re}+{im}i): got {w:?}, want ({wre},{wim})"
5534 );
5535 }
5536 let w = faddeeva_upper_halfplane(Complex { re: 3.0, im: 40.0 });
5540 assert!(
5541 (w.re - 0.01402158696172506).abs() < 1e-13
5542 && (w.im - 0.0010509664408184546).abs() < 1e-13,
5543 "tail value mismatch: w={w:?}"
5544 );
5545 }
5546
5547 #[test]
5548 fn test_integrated_logit_mean_close_to_exact_oracle() {
5549 let ctx = QuadratureContext::new();
5553 let cases = [(-3.0, 0.3), (-1.0, 0.8), (0.5, 1.2), (2.8, 1.0)];
5554 for (eta, se) in cases {
5555 let ghq = logit_posterior_mean(&ctx, eta, se);
5556 let exact = logit_posterior_mean_exact(eta, se);
5557 assert!(
5558 (ghq - exact).abs() < 1e-6,
5559 "production path drifts from oracle at eta={eta} se={se}: \
5560 ghq={ghq:.12} oracle={exact:.12} gap={:.3e}",
5561 (ghq - exact).abs()
5562 );
5563 }
5564 }
5565
5566 #[test]
5567 fn test_probit_posterior_mean_reduces_to_map_atzero_se() {
5568 let eta = 1.25;
5569 let p = probit_posterior_mean(eta, 0.0);
5570 let map = gam_math::probability::normal_cdf(eta);
5571 assert_relative_eq!(p, map, epsilon = 1e-12);
5572 }
5573
5574 #[test]
5575 fn test_probit_posterior_mean_shrinks_extremeswith_uncertainty() {
5576 let hi_eta = 3.0;
5577 let lo_eta = -3.0;
5578 let p_hi_map = probit_posterior_mean(hi_eta, 0.0);
5579 let p_hi_unc = probit_posterior_mean(hi_eta, 2.0);
5580 let p_lo_map = probit_posterior_mean(lo_eta, 0.0);
5581 let p_lo_unc = probit_posterior_mean(lo_eta, 2.0);
5582 assert!(p_hi_unc < p_hi_map);
5583 assert!(p_lo_unc > p_lo_map);
5584 }
5585
5586 #[test]
5587 fn test_survival_posterior_mean_is_bounded_and_shrinks_tail() {
5588 let ctx = QuadratureContext::new();
5589 let eta: f64 = 3.0;
5590 let map = (-(eta.exp())).exp();
5591 let pm = survival_posterior_mean(&ctx, eta, 1.5);
5592 assert!((0.0..=1.0).contains(&pm));
5593 assert!(pm > map);
5594 }
5595
5596 #[test]
5597 fn test_cloglog_and_survival_posterior_means_are_complements() {
5598 let ctx = QuadratureContext::new();
5599 let cases = [
5600 (-3.0, 0.0),
5601 (-0.2, 0.1),
5602 (0.4, 0.8),
5603 (2.0, 1.5),
5604 (10.0, 0.3),
5605 (0.0, 20.0),
5606 (10.0, 10.0),
5607 (-0.5, 100.0),
5608 ];
5609 for (eta, se) in cases {
5610 let clog = cloglog_posterior_mean(&ctx, eta, se);
5611 let surv = survival_posterior_mean(&ctx, eta, se);
5612 assert_relative_eq!(clog + surv, 1.0, epsilon = 2e-10, max_relative = 2e-10);
5613 }
5614 }
5615
5616 #[test]
5617 fn test_cloglog_and_survival_share_large_sigmaspecial_function_path() {
5618 let ctx = QuadratureContext::new();
5619 let eta = -0.2;
5620 let se = 0.8;
5621 let clog = cloglog_posterior_mean(&ctx, eta, se);
5622 let surv = survival_posterior_mean(&ctx, eta, se);
5623 let integrated =
5624 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, eta, se)
5625 .expect("cloglog integrated inverse-link moments should evaluate");
5626 assert_eq!(
5627 integrated.mode,
5628 IntegratedExpectationMode::ExactSpecialFunction
5629 );
5630 assert_relative_eq!(clog, integrated.mean, epsilon = 1e-12, max_relative = 1e-12);
5631 assert_relative_eq!(clog + surv, 1.0, epsilon = 1e-10, max_relative = 1e-10);
5632 }
5633
5634 #[test]
5635 fn test_cloglog_and_survival_posteriorvariances_match() {
5636 let ctx = QuadratureContext::new();
5637 let cases = [(-3.0, 0.0), (-0.2, 0.1), (0.4, 0.8), (2.0, 1.5)];
5638 for (eta, se) in cases {
5639 let (_, clogvar) = cloglog_posterior_meanvariance(&ctx, eta, se);
5640 let (_, survvar) = survival_posterior_meanvariance(&ctx, eta, se);
5641 assert_relative_eq!(clogvar, survvar, epsilon = 1e-12, max_relative = 1e-12);
5642 }
5643 }
5644
5645 #[test]
5646 fn test_survivalvariance_uses_exactsecond_moment_shift() {
5647 let ctx = QuadratureContext::new();
5648 let eta = -0.2;
5649 let se = 0.8;
5650 let (survival, _) = cloglog_survival_term_controlled(&ctx, eta, se);
5651 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(&ctx, eta, se);
5652 let (_, variance) = survival_posterior_meanvariance(&ctx, eta, se);
5653 assert_relative_eq!(
5654 variance,
5655 (survival_sq - survival * survival).max(0.0),
5656 epsilon = 1e-12,
5657 max_relative = 1e-12
5658 );
5659 }
5660
5661 #[test]
5662 fn test_lognormal_laplace_shift_matches_explicitmu_plus_logz() {
5663 let ctx = QuadratureContext::new();
5664 let mu = -0.2;
5665 let sigma = 0.8;
5666 let z = 2.0;
5667 let shifted = lognormal_laplace_term_controlled(&ctx, z, mu, sigma);
5668 let explicit = cloglog_survival_term_controlled(&ctx, mu + z.ln(), sigma);
5669 assert_eq!(shifted.1, explicit.1);
5670 assert_relative_eq!(shifted.0, explicit.0, epsilon = 1e-12, max_relative = 1e-12);
5671 }
5672
5673 #[test]
5674 fn test_integrated_dispatch_uses_closed_form_probit() {
5675 let ctx = QuadratureContext::new();
5676 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Probit, 0.7, 1.3)
5677 .expect("probit integrated inverse-link moments should evaluate");
5678 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
5679 let direct = probit_posterior_meanwith_deriv_exact(0.7, 1.3);
5680 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
5681 assert_relative_eq!(out.dmean_dmu, direct.dmean_dmu, epsilon = 1e-12);
5682 }
5683
5684 #[test]
5685 fn test_integrated_probit_jet_matches_closed_form_derivatives() {
5686 let ctx = QuadratureContext::new();
5687 let mu = 0.7;
5688 let sigma = 1.3;
5689 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Probit, mu, sigma)
5690 .expect("probit integrated inverse-link jet should evaluate");
5691 let s = (1.0 + sigma * sigma).sqrt();
5692 let z = mu / s;
5693 let pdf = gam_math::probability::normal_pdf(z);
5694 assert_relative_eq!(
5695 out.mean,
5696 gam_math::probability::normal_cdf(z),
5697 epsilon = 1e-12
5698 );
5699 assert_relative_eq!(out.d1, pdf / s, epsilon = 1e-12);
5700 assert_relative_eq!(out.d2, -z * pdf / (s * s), epsilon = 1e-12);
5701 assert_relative_eq!(out.d3, (z * z - 1.0) * pdf / (s * s * s), epsilon = 1e-12);
5702 }
5703
5704 #[test]
5705 fn test_integrated_logit_jet_matches_central_differences() {
5706 let ctx = QuadratureContext::new();
5719 let mu = 1.1;
5720 let sigma = 0.8;
5721 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
5722 .expect("logit integrated inverse-link jet should evaluate");
5723 assert!(matches!(
5724 out.mode,
5725 IntegratedExpectationMode::ExactSpecialFunction
5726 | IntegratedExpectationMode::QuadratureFallback
5727 ));
5728 let (ref_mean, ref_d1, ref_d2, ref_d3) = logit_reference_jet_highres_simpson(mu, sigma);
5729 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
5730 assert_relative_eq!(out.d1, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
5731 assert_relative_eq!(out.d2, ref_d2, epsilon = 1e-11, max_relative = 1e-10);
5732 assert_relative_eq!(out.d3, ref_d3, epsilon = 1e-11, max_relative = 1e-10);
5733 }
5734
5735 #[test]
5736 fn test_integrated_cloglog_jet_matches_central_differences() {
5737 let ctx = QuadratureContext::new();
5738 let mu = 0.4;
5739 let sigma = 0.6;
5740 let h = 1e-4;
5741 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
5742 .expect("cloglog integrated inverse-link jet should evaluate");
5743 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
5744 .expect("cloglog integrated inverse-link jet should evaluate");
5745 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
5746 .expect("cloglog integrated inverse-link jet should evaluate");
5747 let d1fd = (plus.mean - minus.mean) / (2.0 * h);
5748 let d2fd = (plus.d1 - minus.d1) / (2.0 * h);
5749 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
5750 assert_eq!(out.d1.signum(), d1fd.signum());
5751 assert_eq!(out.d2.signum(), d2fd.signum());
5752 assert_eq!(out.d3.signum(), d3fd.signum());
5753 assert_relative_eq!(out.d1, d1fd, epsilon = 2e-5, max_relative = 3e-4);
5754 assert_relative_eq!(out.d2, d2fd, epsilon = 4e-5, max_relative = 8e-4);
5755 assert_relative_eq!(out.d3, d3fd, epsilon = 8e-5, max_relative = 2e-3);
5756 }
5757
5758 #[test]
5759 fn test_integrated_cloglog_wide_sigma_d3_matches_simpson_and_d2_slope() {
5760 let ctx = QuadratureContext::new();
5761 let cases = [(0.0, 4.0), (-1.0, 4.0), (2.0, 3.0), (3.0, 3.0)];
5762 let h = 1e-4;
5763
5764 for (mu, sigma) in cases {
5765 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
5766 .expect("wide-sigma cloglog integrated jet should evaluate");
5767 let reference = cloglog_reference_jet_highres_simpson(mu, sigma);
5768 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
5769 .expect("wide-sigma cloglog integrated jet should evaluate");
5770 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
5771 .expect("wide-sigma cloglog integrated jet should evaluate");
5772 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
5773
5774 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
5775 assert_relative_eq!(out.mean, reference.0, epsilon = 4e-8, max_relative = 4e-8);
5776 assert_relative_eq!(out.d1, reference.1, epsilon = 4e-8, max_relative = 4e-8);
5777 assert_relative_eq!(out.d2, reference.2, epsilon = 2e-9, max_relative = 2e-7);
5778 assert_relative_eq!(out.d3, reference.3, epsilon = 2e-9, max_relative = 2e-7);
5779 assert_relative_eq!(out.d3, d3fd, epsilon = 2e-7, max_relative = 4e-5);
5780 }
5781 }
5782
5783 #[test]
5784 fn test_latent_cloglog_jet5_matches_higher_order_central_differences() {
5785 let ctx = QuadratureContext::new();
5786 let mu = 0.35;
5787 let sigma = 0.7;
5788 let h = 2e-4;
5789
5790 let out = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu, sigma);
5791 let plus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu + h, sigma);
5792 let minus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu - h, sigma);
5793
5794 let d4fd = (plus.d3 - minus.d3) / (2.0 * h);
5795 let d5fd = (plus.d4 - minus.d4) / (2.0 * h);
5796
5797 assert_eq!(out.d4.signum(), d4fd.signum());
5798 assert_eq!(out.d5.signum(), d5fd.signum());
5799 assert_relative_eq!(out.d4, d4fd, epsilon = 2e-4, max_relative = 5e-3);
5800 assert_relative_eq!(out.d5, d5fd, epsilon = 6e-4, max_relative = 2e-2);
5801 }
5802
5803 #[test]
5804 fn test_logit_exact_derivative_matches_finite_difference() {
5805 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("controlled logit");
5815 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
5816 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
5817 assert!(out.dmean_dmu > 0.0);
5818 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
5819 }
5820
5821 #[test]
5822 fn test_logit_exact_clamped_degenerate_branch_is_locally_flat() {
5823 let out = logit_posterior_meanwith_deriv_exact(-710.0, 0.0).expect("exact logit");
5824 let h = 1e-6;
5825 let plus = logit_posterior_meanwith_deriv_exact(-710.0 + h, 0.0)
5826 .expect("exact logit plus")
5827 .mean;
5828 let minus = logit_posterior_meanwith_deriv_exact(-710.0 - h, 0.0)
5829 .expect("exact logit minus")
5830 .mean;
5831 let fd = (plus - minus) / (2.0 * h);
5832 assert_eq!(fd, 0.0);
5833 assert_eq!(out.dmean_dmu, 0.0);
5834 }
5835
5836 fn simpson_integrate<F>(a: f64, b: f64, n_intervals: usize, f: F) -> f64
5837 where
5838 F: Fn(f64) -> f64,
5839 {
5840 assert_eq!(n_intervals % 2, 0, "Simpson integration requires an even n");
5841 let h = (b - a) / n_intervals as f64;
5842 let mut sum = f(a) + f(b);
5843 for i in 1..n_intervals {
5844 let x = a + i as f64 * h;
5845 let w = if i % 2 == 0 { 2.0 } else { 4.0 };
5846 sum += w * f(x);
5847 }
5848 sum * h / 3.0
5849 }
5850
5851 fn cloglog_reference_mean_and_derivative(mu: f64, sigma: f64) -> (f64, f64) {
5852 if sigma <= CLOGLOG_SIGMA_DEGENERATE {
5853 return (cloglog_mean_exact(mu), cloglog_mean_d1_exact(mu));
5854 }
5855
5856 let z_max = 12.0;
5860 let n_intervals = 4096;
5861 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5862 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5863 let eta = mu + sigma * z;
5864 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_exact(eta)
5865 });
5866 let deriv = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5867 let eta = mu + sigma * z;
5868 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_d1_exact(eta)
5869 });
5870 (mean, deriv)
5871 }
5872
5873 fn logit_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
5886 let z_max = 14.0;
5887 let n_intervals = 16384;
5888 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5889 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
5890 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5891 let eta = mu + sigma * z;
5892 let (p, _, _, _) = component_point_jet(LinkComponent::Logit, eta);
5893 phi(z) * p
5894 });
5895 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5896 let eta = mu + sigma * z;
5897 let (_, p1, _, _) = component_point_jet(LinkComponent::Logit, eta);
5898 phi(z) * p1
5899 });
5900 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5901 let eta = mu + sigma * z;
5902 let (_, _, p2, _) = component_point_jet(LinkComponent::Logit, eta);
5903 phi(z) * p2
5904 });
5905 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5906 let eta = mu + sigma * z;
5907 let (_, _, _, p3) = component_point_jet(LinkComponent::Logit, eta);
5908 phi(z) * p3
5909 });
5910 (mean, d1, d2, d3)
5911 }
5912
5913 fn cloglog_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
5914 let z_max = 14.0;
5915 let n_intervals = 16384;
5916 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5917 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
5918 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5919 let eta = mu + sigma * z;
5920 let (g, _, _, _, _, _) = cloglog_point_jet5(eta);
5921 phi(z) * g
5922 });
5923 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5924 let eta = mu + sigma * z;
5925 let (_, g1, _, _, _, _) = cloglog_point_jet5(eta);
5926 phi(z) * g1
5927 });
5928 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5929 let eta = mu + sigma * z;
5930 let (_, _, g2, _, _, _) = cloglog_point_jet5(eta);
5931 phi(z) * g2
5932 });
5933 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5934 let eta = mu + sigma * z;
5935 let (_, _, _, g3, _, _) = cloglog_point_jet5(eta);
5936 phi(z) * g3
5937 });
5938 (mean, d1, d2, d3)
5939 }
5940
5941 #[test]
5942 fn test_cloglog_taylor_negative_tail_matches_mathematical_target() {
5943 let mu = -40.0;
5944 let sigma = 0.1;
5945 let out = cloglog_small_sigma_taylor(mu, sigma);
5946 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
5947
5948 assert!(
5949 out.dmean_dmu > 0.0,
5950 "negative-tail derivative should remain positive"
5951 );
5952 assert_relative_eq!(
5953 out.mean,
5954 expected_mean,
5955 epsilon = 1e-30,
5956 max_relative = 1e-12
5957 );
5958 assert_relative_eq!(
5959 out.dmean_dmu,
5960 expected_deriv,
5961 epsilon = 1e-30,
5962 max_relative = 1e-12
5963 );
5964 }
5965
5966 #[test]
5967 fn test_cloglog_degenerate_negative_tail_matches_pointwise_target() {
5968 let ctx = QuadratureContext::new();
5969 let mu = -40.0;
5970 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, 0.0);
5971
5972 assert!(
5973 out.dmean_dmu > 0.0,
5974 "degenerate negative-tail derivative should remain positive"
5975 );
5976 assert_relative_eq!(
5977 out.mean,
5978 cloglog_mean_exact(mu),
5979 epsilon = 1e-30,
5980 max_relative = 1e-15
5981 );
5982 assert_relative_eq!(
5983 out.dmean_dmu,
5984 cloglog_mean_d1_exact(mu),
5985 epsilon = 1e-30,
5986 max_relative = 1e-15
5987 );
5988 }
5989
5990 #[test]
5991 fn test_degenerate_probit_jet_is_exact_beyond_former_clamp() {
5992 let mu = -30.1;
5993 let probit = integrated_probit_jet(mu, 0.0);
5994 let pdf = gam_math::probability::normal_pdf(mu);
5995 assert!(
5996 pdf > 0.0,
5997 "test point must have a represented Gaussian tail"
5998 );
5999 assert_eq!(probit.mean, gam_math::probability::normal_cdf(mu));
6000 assert_eq!(probit.d1, pdf);
6001 assert_eq!(probit.d2, -mu * pdf);
6002 assert_eq!(probit.d3, (mu * mu - 1.0) * pdf);
6003
6004 let tail = (-710.0_f64).exp();
6020 assert!(
6021 tail > 0.0 && tail < f64::MIN_POSITIVE,
6022 "eta=-710 must sit in the subnormal tail, not underflow"
6023 );
6024 let logit = component_point_jet(LinkComponent::Logit, -710.0);
6025 assert_eq!(logit.0, tail);
6026 assert_eq!(logit.1, tail);
6027 assert_eq!(logit.2, tail);
6028 assert_eq!(logit.3, tail);
6029
6030 assert_eq!(
6032 (-750.0_f64).exp(),
6033 0.0,
6034 "eta=-750 must underflow f64 for this arm to mean anything"
6035 );
6036 let underflowed = component_point_jet(LinkComponent::Logit, -750.0);
6037 assert_eq!(underflowed.1, 0.0);
6038 assert_eq!(underflowed.2, 0.0);
6039 assert_eq!(underflowed.3, 0.0);
6040 }
6041
6042 #[test]
6043 fn test_degenerate_cloglog_component_jet_preserves_smooth_negative_tail() {
6044 let eta: f64 = -40.0;
6045 let t = eta.exp();
6046 let s = (-t).exp();
6047 let cloglog = component_point_jet(LinkComponent::CLogLog, eta);
6048 let expected_mean = -(-t).exp_m1();
6049 let expected_d1 = t * s;
6050 let expected_d2 = (t - t * t) * s;
6051 let expected_d3 = (t - 3.0 * t * t + t * t * t) * s;
6052
6053 assert!(cloglog.1 > 0.0, "negative-tail d1 should remain positive");
6054 assert_relative_eq!(
6055 cloglog.0,
6056 expected_mean,
6057 epsilon = 1e-30,
6058 max_relative = 1e-15
6059 );
6060 assert_relative_eq!(
6061 cloglog.1,
6062 expected_d1,
6063 epsilon = 1e-30,
6064 max_relative = 1e-15
6065 );
6066 assert_relative_eq!(
6067 cloglog.2,
6068 expected_d2,
6069 epsilon = 1e-30,
6070 max_relative = 1e-15
6071 );
6072 assert_relative_eq!(
6073 cloglog.3,
6074 expected_d3,
6075 epsilon = 1e-30,
6076 max_relative = 1e-15
6077 );
6078 }
6079
6080 #[test]
6081 fn test_zero_sigma_logit_and_cloglog_share_component_tail_jets() {
6082 let ctx = QuadratureContext::new();
6083 for (link, component, eta) in [
6084 (LinkFunction::Logit, LinkComponent::Logit, 50.0),
6085 (LinkFunction::CLogLog, LinkComponent::CLogLog, -50.0),
6086 ] {
6087 let integrated = integrated_inverse_link_jet(&ctx, link, eta, 0.0)
6088 .expect("degenerate integrated jet");
6089 let point = component_inverse_link_jet(component, eta);
6090 assert_eq!(integrated.mode, IntegratedExpectationMode::ExactClosedForm);
6091 assert_eq!(integrated.mean, point.mu);
6092 assert_eq!(integrated.d1, point.d1);
6093 assert_eq!(integrated.d2, point.d2);
6094 assert_eq!(integrated.d3, point.d3);
6095 }
6096 }
6097
6098 #[test]
6099 fn test_cloglog_controlled_matches_mathematical_target_on_small_sigma_grid() {
6100 let ctx = QuadratureContext::new();
6101 let cases = [
6105 (-30.0, 1e-10),
6106 (-30.0, 0.1),
6107 (-10.0, 0.24),
6108 (-3.0, 0.2),
6109 (0.0, 0.05),
6110 (0.4, 0.1),
6111 (3.0, 0.24),
6112 (10.0, 0.1),
6113 (30.0, 0.24),
6114 ];
6115
6116 for &(mu, sigma) in &cases {
6117 let approx = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
6118 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
6119 assert_relative_eq!(
6120 approx.mean,
6121 expected_mean,
6122 epsilon = 1e-12,
6123 max_relative = 2e-3
6124 );
6125 assert_relative_eq!(
6126 approx.dmean_dmu,
6127 expected_deriv,
6128 epsilon = 1e-12,
6129 max_relative = 4e-3
6130 );
6131 }
6132 }
6133
6134 #[test]
6135 fn test_cloglog_dispatch_uses_gamma_backend_for_large_sigma_central_regime() {
6136 let ctx = QuadratureContext::new();
6137 let out =
6138 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, -0.2, 0.8)
6139 .expect("cloglog integrated inverse-link moments should evaluate");
6140 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6141 assert!(out.mean.is_finite());
6142 assert!(out.dmean_dmu.is_finite());
6143 assert!(out.dmean_dmu >= 0.0);
6144 }
6145
6146 #[test]
6147 fn test_cloglog_dispatch_uses_large_sigma_asymptotic_without_ghq() {
6148 let ctx = QuadratureContext::new();
6149 let out =
6150 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, 0.0, 20.0)
6151 .expect("cloglog integrated inverse-link moments should evaluate");
6152 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6153 assert!(out.mean.is_finite());
6154 assert!(out.dmean_dmu.is_finite());
6155 assert!(out.dmean_dmu >= 0.0);
6156 }
6157
6158 #[test]
6159 fn test_cloglog_cc_matches_gamma_reference_on_central_case() {
6160 let ctx = QuadratureContext::new();
6161 let mu = -0.2;
6162 let sigma = 0.8;
6163 let cc = cloglog_survival_cc(&ctx, mu, sigma, CLOGLOG_CC_TOL).expect("cc backend");
6164 let gamma = cloglog_survival_gamma_reference(mu, sigma).expect("gamma backend");
6165 assert_relative_eq!(cc, gamma, epsilon = 5e-6, max_relative = 5e-6);
6166 }
6167
6168 #[test]
6169 fn test_cloglog_gamma_reference_matches_seeded_monte_carlo_small_case() {
6170 let mu = -0.2;
6171 let sigma = 0.8;
6172 let gamma =
6173 cloglog_posterior_meanwith_deriv_gamma_reference(mu, sigma).expect("gamma reference");
6174 let mut rng_state = 0x9e3779b97f4a7c15u64;
6175 let mut mean_mc = 0.0f64;
6176 let mut deriv_mc = 0.0f64;
6177 let n_samples = 300_000usize;
6178 for _ in 0..n_samples {
6179 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6180 let u1 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6181 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6182 let u2 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6183 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
6184 let eta = mu + sigma * z;
6185 mean_mc += cloglog_mean_exact(eta);
6186 deriv_mc += cloglog_mean_d1_exact(eta);
6187 }
6188 mean_mc /= n_samples as f64;
6189 deriv_mc /= n_samples as f64;
6190 assert_relative_eq!(gamma.mean, mean_mc, epsilon = 2e-3, max_relative = 2e-3);
6191 assert_relative_eq!(
6192 gamma.dmean_dmu,
6193 deriv_mc,
6194 epsilon = 2e-3,
6195 max_relative = 2e-3
6196 );
6197 }
6198
6199 #[test]
6200 fn test_logit_dispatch_uses_tail_asymptotic_outside_old_guard() {
6201 let ctx = QuadratureContext::new();
6202 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 35.0, 1.0)
6203 .expect("logit integrated inverse-link moments should evaluate");
6204 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6205 assert!(out.mean.is_finite());
6206 assert!(out.dmean_dmu.is_finite());
6207 assert!(out.dmean_dmu >= 0.0);
6208 }
6209
6210 #[test]
6211 fn test_logit_dispatch_prefers_erfcx_in_moderate_regime() {
6212 let ctx = QuadratureContext::new();
6223 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 1.1, 0.8)
6224 .expect("logit integrated inverse-link moments should evaluate");
6225 assert!(matches!(
6226 out.mode,
6227 IntegratedExpectationMode::ExactSpecialFunction
6228 | IntegratedExpectationMode::QuadratureFallback
6229 ));
6230 assert!(out.mean.is_finite());
6231 assert!(out.dmean_dmu.is_finite());
6232 assert!(out.dmean_dmu >= 0.0);
6233 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6234 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6235 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6236 }
6237
6238 #[test]
6239 fn test_logit_dispatch_large_sigma_uses_accurate_quadrature_not_monahan() {
6240 let ctx = QuadratureContext::new();
6249 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 0.5, 20.0)
6250 .expect("logit integrated inverse-link moments should evaluate");
6251 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
6252 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(0.5, 20.0);
6253 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6254 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6255 let kappa = (1.0 + std::f64::consts::PI * 20.0 * 20.0 / 8.0)
6258 .sqrt()
6259 .recip();
6260 let monahan_mean = gam_math::probability::normal_cdf(0.5 * kappa);
6261 assert!(
6262 (out.mean - monahan_mean).abs() > 1e-3,
6263 "dispatcher must not return the inaccurate Monahan mean {monahan_mean}; got {}",
6264 out.mean
6265 );
6266 }
6267
6268 #[test]
6269 fn test_logit_controlled_path_keeps_exact_backend_in_moderate_regime() {
6270 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("logit controlled");
6280 assert!(matches!(
6281 out.mode,
6282 IntegratedExpectationMode::ExactSpecialFunction
6283 | IntegratedExpectationMode::QuadratureFallback
6284 ));
6285 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6286 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6287 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6288 }
6289
6290 #[test]
6291 fn test_logit_dispatch_derivative_correct_at_mu_zero_small_sigma() {
6292 let ctx = QuadratureContext::new();
6301 for &(mu, sigma) in &[(0.0, 0.3), (0.0, 0.4), (0.0, 0.5)] {
6302 let out =
6303 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6304 .expect("logit integrated inverse-link moments should evaluate");
6305 assert_relative_eq!(out.mean, 0.5, epsilon = 1e-10);
6307 assert!(
6309 out.dmean_dmu <= 0.25 + 1e-9,
6310 "E[sigmoid'] must not exceed 0.25 at (μ={mu}, σ={sigma}); got {}",
6311 out.dmean_dmu
6312 );
6313 let (_, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6314 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-6);
6315 }
6316 }
6317
6318 #[test]
6319 fn test_logit_erfcx_exact_branch_is_self_certified() {
6320 for &(mu, sigma) in &[(8.0, 1.0), (10.0, 1.0), (15.0, 2.0)] {
6327 let out = logit_posterior_meanwith_deriv_exact(mu, sigma)
6328 .expect("erfcx branch should certify");
6329 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6330 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6331 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6332 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6333 }
6334 assert!(
6338 logit_posterior_meanwith_deriv_exact(0.0, 0.3).is_err(),
6339 "erfcx branch must not claim ExactSpecialFunction when it cannot certify the derivative"
6340 );
6341 }
6342
6343 #[test]
6344 fn test_logit_integrated_derivative_is_even_in_mu() {
6345 let ctx = QuadratureContext::new();
6351 for &(mu, sigma) in &[(0.3, 0.3), (1.1, 0.8), (10.0, 1.0), (3.0, 3.0), (35.0, 1.0)] {
6352 let pos =
6353 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6354 .expect("logit moments (+μ)");
6355 let neg =
6356 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, -mu, sigma)
6357 .expect("logit moments (-μ)");
6358 assert_relative_eq!(
6359 pos.dmean_dmu,
6360 neg.dmean_dmu,
6361 epsilon = 1e-9,
6362 max_relative = 1e-7
6363 );
6364 assert_relative_eq!(
6366 neg.mean,
6367 1.0 - pos.mean,
6368 epsilon = 1e-9,
6369 max_relative = 1e-7
6370 );
6371 }
6372 }
6373
6374 #[test]
6375 fn test_logit_dmean_dmu_equals_fd_of_mean_across_regimes() {
6376 let ctx = QuadratureContext::new();
6391 let h = 1e-4;
6392 let cases = [
6393 (0.0, 0.8), (0.7, 0.8), (1.5, 1.2), (-1.1, 0.9), (8.0, 1.0), (10.0, 1.5), (-9.0, 1.0), (0.5, 0.05), (0.5, 20.0), ];
6403 for &(mu, sigma) in &cases {
6404 let at = |m: f64| {
6405 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, m, sigma)
6406 .expect("logit moments")
6407 };
6408 let out = at(mu);
6409 let fd = (at(mu + h).mean - at(mu - h).mean) / (2.0 * h);
6410 assert!(
6411 (out.dmean_dmu - fd).abs() <= 1e-5,
6412 "dmean_dmu must equal d/dμ of mean at (μ={mu}, σ={sigma}): \
6413 returned {}, FD of mean {} (mode {:?})",
6414 out.dmean_dmu,
6415 fd,
6416 out.mode
6417 );
6418 assert!(
6422 out.dmean_dmu <= 0.25 + 1e-9 && out.dmean_dmu >= 0.0,
6423 "dmean_dmu out of [0, 0.25] at (μ={mu}, σ={sigma}): {}",
6424 out.dmean_dmu
6425 );
6426 }
6427 }
6428
6429 #[test]
6430 fn test_logit_scalar_matches_jet_at_large_sigma() {
6431 let ctx = QuadratureContext::new();
6437 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0)] {
6438 let scalar =
6439 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6440 .expect("scalar logit moments");
6441 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6442 .expect("jet logit moments");
6443 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6447 assert_relative_eq!(scalar.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-8);
6448 assert_relative_eq!(
6449 scalar.dmean_dmu,
6450 ref_d1,
6451 epsilon = 1e-9,
6452 max_relative = 1e-8
6453 );
6454 assert_relative_eq!(scalar.mean, jet.mean, epsilon = 1e-12, max_relative = 1e-12);
6462 assert_relative_eq!(
6463 scalar.dmean_dmu,
6464 jet.d1,
6465 epsilon = 1e-12,
6466 max_relative = 1e-12
6467 );
6468 }
6469 }
6470
6471 #[test]
6472 fn test_logit_jet_accurate_at_wide_sigma() {
6473 let ctx = QuadratureContext::new();
6482 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0), (0.5, 20.0)] {
6483 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6484 .expect("wide-σ logit jet");
6485 let (rm, rd1, rd2, rd3) = logit_reference_jet_highres_simpson(mu, sigma);
6486 assert_relative_eq!(jet.mean, rm, epsilon = 1e-8, max_relative = 1e-7);
6487 assert_relative_eq!(jet.d1, rd1, epsilon = 1e-8, max_relative = 1e-6);
6488 assert_relative_eq!(jet.d2, rd2, epsilon = 1e-8, max_relative = 1e-6);
6489 assert_relative_eq!(jet.d3, rd3, epsilon = 1e-8, max_relative = 1e-6);
6490 let scalar =
6492 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6493 .expect("scalar logit moments");
6494 assert_relative_eq!(jet.d1, scalar.dmean_dmu, epsilon = 1e-12);
6495 assert_relative_eq!(jet.mean, scalar.mean, epsilon = 1e-12);
6496 }
6497 }
6498
6499 #[test]
6500 fn test_logit_jet_continuous_across_ghq_simpson_seam() {
6501 let ctx = QuadratureContext::new();
6509 let sigma = LOGIT_JET_GHQ_SIGMA_MAX;
6510 for mu in [-2.0, -0.5, 0.0, 0.7, 1.3, 3.0] {
6511 let ghq = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6513 .expect("jet at seam (GHQ dispatch)");
6514 let simpson = logit_wide_sigma_jet(mu, sigma).expect("jet at seam (Simpson)");
6516 assert_relative_eq!(ghq.mean, simpson.mean, epsilon = 1e-9, max_relative = 1e-8);
6519 assert_relative_eq!(ghq.d1, simpson.d1, epsilon = 1e-9, max_relative = 1e-7);
6520 assert_relative_eq!(ghq.d2, simpson.d2, epsilon = 1e-9, max_relative = 1e-7);
6521 assert_relative_eq!(ghq.d3, simpson.d3, epsilon = 1e-8, max_relative = 1e-6);
6522 }
6523 }
6524
6525 #[test]
6526 fn test_logit_batch_uses_same_dispatchvalues() {
6527 let ctx = QuadratureContext::new();
6528 let eta = ndarray::array![-2.0, 0.0, 1.25, 35.0];
6529 let se = ndarray::array![0.1, 0.5, 1.0, 1.0];
6530 let batch_mean = logit_posterior_mean_batch(&ctx, &eta, &se)
6531 .expect("logit posterior mean batch should evaluate");
6532 let (batchmu, batch_dmu) = logit_posterior_meanwith_deriv_batch(&ctx, &eta, &se)
6533 .expect("logit posterior mean derivative batch should evaluate");
6534 for i in 0..eta.len() {
6535 let direct = integrated_inverse_link_mean_and_derivative(
6536 &ctx,
6537 LinkFunction::Logit,
6538 eta[i],
6539 se[i],
6540 )
6541 .expect("logit integrated inverse-link moments should evaluate");
6542 assert_relative_eq!(batch_mean[i], direct.mean, epsilon = 1e-12);
6543 assert_relative_eq!(batchmu[i], direct.mean, epsilon = 1e-12);
6544 assert_relative_eq!(batch_dmu[i], direct.dmean_dmu, epsilon = 1e-12);
6545 }
6546 }
6547
6548 #[test]
6549 fn exact_logit_small_se_branch_loses_tail_derivative() {
6550 let eta = 50.0_f64;
6551 let stable_z = (-eta).exp();
6552 let stable_dmu = stable_z / (1.0_f64 + stable_z).powi(2);
6553 assert!(stable_dmu > 0.0);
6554 let out = logit_posterior_meanwith_deriv_exact(eta, 0.0).expect("exact branch");
6555 let dmu = out.dmean_dmu;
6556 assert!(
6557 (dmu - stable_dmu).abs() < 1e-30,
6558 "exact logit small-se branch should use the stable derivative z/(1+z)^2 at eta={eta}; got {} vs {}",
6559 dmu,
6560 stable_dmu
6561 );
6562 }
6563
6564 #[test]
6565 fn integrated_family_moments_rejects_latent_cloglog_without_concrete_handler() {
6566 let ctx = QuadratureContext::new();
6572 let latent =
6573 gam_problem::types::LatentCLogLogState::new(0.4).expect("valid latent cloglog state");
6574 let spec =
6575 LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::LatentCLogLog(latent));
6576 let likelihood = GlmLikelihoodSpec::canonical(spec);
6577 let err = integrated_family_moments_jet(
6578 &ctx,
6579 &likelihood,
6580 0.2,
6581 0.5,
6582 )
6583 .expect_err("latent cloglog moments should error in this dispatcher");
6584 assert!(format!("{err}").contains("LatentCLogLog"));
6585 }
6586
6587 #[test]
6588 fn integrated_family_moments_supports_stateful_sas() {
6589 let ctx = QuadratureContext::new();
6590 let sas = crate::mixture_link::state_from_sasspec(gam_problem::types::SasLinkSpec {
6591 initial_epsilon: 0.3,
6592 initial_log_delta: -0.2,
6593 })
6594 .expect("sas state should reconstruct from raw parameters");
6595 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Sas(sas));
6596 let likelihood = GlmLikelihoodSpec::canonical(spec);
6597 let out = integrated_family_moments_jet(
6598 &ctx,
6599 &likelihood,
6600 0.2,
6601 0.5,
6602 )
6603 .expect("stateful SAS integrated moments should evaluate");
6604 assert!(out.mean.is_finite());
6605 assert!(out.d1.is_finite());
6606 assert!(out.d2.is_finite());
6607 assert!(out.d3.is_finite());
6608 assert!(out.mean > 0.0 && out.mean < 1.0);
6609 }
6610
6611 #[test]
6612 fn integrated_family_moments_supports_pure_probit_mixture() {
6613 let ctx = QuadratureContext::new();
6614 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6615 components: vec![gam_problem::types::LinkComponent::Probit],
6616 initial_rho: ndarray::Array1::<f64>::zeros(0),
6617 })
6618 .expect("single-component probit mixture state");
6619 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
6620 let likelihood = GlmLikelihoodSpec::canonical(spec);
6621 let out = integrated_family_moments_jet(
6622 &ctx,
6623 &likelihood,
6624 0.7,
6625 1.3,
6626 )
6627 .expect("pure probit mixture integrated moments should evaluate");
6628 let exact = integrated_probit_jet(0.7, 1.3);
6629 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
6630 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
6631 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
6632 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
6633 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
6634 }
6635
6636 #[test]
6637 fn integrated_family_moments_supports_pure_logit_mixture() {
6638 let ctx = QuadratureContext::new();
6639 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6640 components: vec![gam_problem::types::LinkComponent::Logit],
6641 initial_rho: ndarray::Array1::<f64>::zeros(0),
6642 })
6643 .expect("single-component logit mixture state");
6644 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
6645 let likelihood = GlmLikelihoodSpec::canonical(spec);
6646 let out = integrated_family_moments_jet(
6647 &ctx,
6648 &likelihood,
6649 1.1,
6650 0.8,
6651 )
6652 .expect("pure logit mixture integrated moments should evaluate");
6653 let exact = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, 1.1, 0.8)
6654 .expect("canonical integrated logit jet");
6655 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
6656 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
6657 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
6658 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
6659 assert_eq!(out.mode, exact.mode);
6660 }
6661
6662 #[test]
6663 fn integrated_family_moments_supports_stateful_mixture() {
6664 let ctx = QuadratureContext::new();
6665 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6666 components: vec![
6667 gam_problem::types::LinkComponent::Logit,
6668 gam_problem::types::LinkComponent::Probit,
6669 ],
6670 initial_rho: ndarray::array![0.35],
6671 })
6672 .expect("mixture state should reconstruct from rho");
6673 let spec = LikelihoodSpec::new(
6674 ResponseFamily::Binomial,
6675 InverseLink::Mixture(state.clone()),
6676 );
6677 let likelihood = GlmLikelihoodSpec::canonical(spec);
6678 let out = integrated_family_moments_jet(
6679 &ctx,
6680 &likelihood,
6681 0.2,
6682 0.5,
6683 )
6684 .expect("stateful mixture integrated moments should evaluate");
6685 let direct = integrated_mixture_jet(&ctx, 0.2, 0.5, &state)
6686 .expect("direct integrated mixture jet should evaluate");
6687 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
6688 assert_relative_eq!(out.d1, direct.d1, epsilon = 1e-12);
6689 assert_relative_eq!(out.d2, direct.d2, epsilon = 1e-12);
6690 assert_relative_eq!(out.d3, direct.d3, epsilon = 1e-12);
6691 assert_eq!(out.mode, direct.mode);
6692 }
6693
6694 #[test]
6695 fn integrated_family_moments_use_scale_dispersion_for_tweedie_and_gamma() {
6696 let ctx = QuadratureContext::new();
6700 let e = 0.3_f64;
6702 let se = 0.5_f64;
6703 let m = (e + 0.5 * se * se).exp();
6704
6705 let p = 1.5_f64;
6707 let phi = 2.0_f64;
6708 let tweedie = LikelihoodSpec::tweedie_log(p);
6709 let tweedie_likelihood = GlmLikelihoodSpec {
6710 spec: tweedie.clone(),
6711 scale: LikelihoodScaleMetadata::EstimatedTweediePhi { phi },
6712 };
6713 let out = integrated_family_moments_jet(
6714 &ctx,
6715 &tweedie_likelihood,
6716 e,
6717 se,
6718 )
6719 .expect("tweedie integrated moments should evaluate");
6720 let expected = phi * m.powf(p);
6721 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
6722 assert_relative_eq!(out.variance / m.powf(p), phi, epsilon = 1e-12);
6724
6725 let shape = 4.0_f64;
6727 let gamma = LikelihoodSpec::gamma_log();
6728 let gamma_likelihood = GlmLikelihoodSpec {
6729 spec: gamma.clone(),
6730 scale: LikelihoodScaleMetadata::EstimatedGammaShape { shape },
6731 };
6732 let out = integrated_family_moments_jet(
6733 &ctx,
6734 &gamma_likelihood,
6735 e,
6736 se,
6737 )
6738 .expect("gamma integrated moments should evaluate");
6739 let expected = m * m / shape;
6740 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
6741 assert_relative_eq!(out.variance / (m * m), 1.0 / shape, epsilon = 1e-12);
6743
6744 let poisson = LikelihoodSpec::poisson_log();
6746 let poisson_likelihood = GlmLikelihoodSpec::canonical(poisson);
6747 let out = integrated_family_moments_jet(
6748 &ctx,
6749 &poisson_likelihood,
6750 e,
6751 se,
6752 )
6753 .expect("poisson integrated moments should evaluate");
6754 assert_relative_eq!(out.variance, m, epsilon = 1e-12);
6755
6756 let theta = 3.0_f64;
6758 let nb = LikelihoodSpec::negative_binomial_log(theta);
6759 let nb_likelihood = GlmLikelihoodSpec::canonical(nb);
6760 let out = integrated_family_moments_jet(
6761 &ctx,
6762 &nb_likelihood,
6763 e,
6764 se,
6765 )
6766 .expect("negative-binomial integrated moments should evaluate");
6767 assert_relative_eq!(out.variance, m + m * m / theta, epsilon = 1e-12);
6768
6769 let missing_gamma = GlmLikelihoodSpec {
6771 spec: gamma,
6772 scale: LikelihoodScaleMetadata::Unspecified,
6773 };
6774 let err = integrated_family_moments_jet(
6775 &ctx,
6776 &missing_gamma,
6777 e,
6778 se,
6779 )
6780 .expect_err("gamma without a shape in the scale metadata must error");
6781 assert!(
6782 format!("{err}").contains("GammaShape"),
6783 "unexpected error message: {err}"
6784 );
6785
6786 let missing_tweedie = GlmLikelihoodSpec {
6788 spec: tweedie,
6789 scale: LikelihoodScaleMetadata::Unspecified,
6790 };
6791 let err = integrated_family_moments_jet(
6792 &ctx,
6793 &missing_tweedie,
6794 e,
6795 se,
6796 )
6797 .expect_err("tweedie without a φ in the scale metadata must error");
6798 assert!(
6799 format!("{err}").contains("EstimatedTweediePhi"),
6800 "unexpected error message: {err}"
6801 );
6802 }
6803
6804 #[test]
6807 fn cloglog_g_derivatives_at_zero() {
6808 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(0.0);
6809 let expected_g = 1.0 - (-1.0_f64).exp();
6811 assert_relative_eq!(g, expected_g, epsilon = 1e-14);
6812 let e_neg1 = (-1.0_f64).exp();
6814 assert_relative_eq!(g1, e_neg1, epsilon = 1e-14);
6815 assert_relative_eq!(g2, 0.0, epsilon = 1e-14);
6817 assert_relative_eq!(g3, -e_neg1, epsilon = 1e-14);
6819 assert_relative_eq!(g4, -e_neg1, epsilon = 1e-14);
6821 }
6822
6823 #[test]
6824 fn cloglog_g_derivatives_saturation() {
6825 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(50.0);
6827 assert_relative_eq!(g, 1.0, epsilon = 1e-10);
6828 assert_eq!(g1, 0.0);
6829 assert_eq!(g2, 0.0);
6830 assert_eq!(g3, 0.0);
6831 assert_eq!(g4, 0.0);
6832
6833 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(-50.0);
6835 let expected = (-50.0_f64).exp();
6836 assert_relative_eq!(g, expected, max_relative = 1e-10);
6837 assert_relative_eq!(g1, expected, max_relative = 1e-10);
6838 assert_relative_eq!(g2, expected, max_relative = 1e-10);
6840 assert_relative_eq!(g3, expected, max_relative = 1e-10);
6841 assert_relative_eq!(g4, expected, max_relative = 1e-10);
6842 }
6843
6844 #[test]
6845 fn cloglog_ghq_value_sigma_zero_matches_pointwise() {
6846 let ctx = QuadratureContext::new();
6847 for &mu in &[-2.0, -1.0, 0.0, 0.5, 1.5] {
6849 let val = cloglog_ghq_value(&ctx, mu, 0.0, 21);
6850 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
6851 assert_relative_eq!(val, g, epsilon = 1e-14);
6852 }
6853 }
6854
6855 #[test]
6856 fn cloglog_ghq_value_bounded_zero_one() {
6857 let ctx = QuadratureContext::new();
6858 for &mu in &[-5.0, -2.0, 0.0, 1.0, 3.0, 10.0] {
6860 for &sigma in &[0.1, 0.5, 1.0, 2.0, 5.0] {
6861 let val = cloglog_ghq_value(&ctx, mu, sigma, 31);
6862 assert!((0.0..=1.0).contains(&val), "L({mu},{sigma}) = {val}");
6863 }
6864 }
6865 }
6866
6867 #[test]
6868 fn cloglog_ghq_derivatives_sigma_zero_matches_pointwise() {
6869 let ctx = QuadratureContext::new();
6870 let mu = 0.3;
6871 let d = cloglog_ghq_derivatives(&ctx, mu, 0.0, 21);
6872 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
6873 assert_relative_eq!(d.l, g, epsilon = 1e-14);
6874 assert_relative_eq!(d.l_mu, g1, epsilon = 1e-14);
6875 assert_relative_eq!(d.l_mumu, g2, epsilon = 1e-14);
6876 assert_relative_eq!(d.l_mumumu, g3, epsilon = 1e-14);
6877 assert_relative_eq!(d.l_mumumumu, g4, epsilon = 1e-14);
6878
6879 assert_eq!(d.l_sigma, 0.0);
6881 assert_eq!(d.l_musigma, 0.0);
6882 assert_eq!(d.l_mumusigma, 0.0);
6883 assert_eq!(d.l_mumumusigma, 0.0);
6884 assert_eq!(d.l_sigmasigmasigma, 0.0);
6885 assert_eq!(d.l_musigmasigmasigma, 0.0);
6886
6887 assert_relative_eq!(d.l_sigmasigma, g2, epsilon = 1e-14);
6890 assert_relative_eq!(d.l_musigmasigma, g3, epsilon = 1e-14);
6891 assert_relative_eq!(d.l_mumusigmasigma, g4, epsilon = 1e-14);
6892 assert_relative_eq!(d.l_sigmasigmasigmasigma, 3.0 * g4, epsilon = 1e-14);
6893 }
6894
6895 #[test]
6896 fn cloglog_ghq_derivatives_finite_difference_mu() {
6897 let ctx = QuadratureContext::new();
6899 let mu = 0.5;
6900 let sigma = 0.8;
6901 let h = 1e-6;
6902 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
6903 let l_plus = cloglog_ghq_value(&ctx, mu + h, sigma, 31);
6904 let l_minus = cloglog_ghq_value(&ctx, mu - h, sigma, 31);
6905 let fd_mu = (l_plus - l_minus) / (2.0 * h);
6906 assert_relative_eq!(d.l_mu, fd_mu, epsilon = 1e-5);
6907
6908 let d_plus = cloglog_ghq_derivatives(&ctx, mu + h, sigma, 31);
6910 let d_minus = cloglog_ghq_derivatives(&ctx, mu - h, sigma, 31);
6911 let fd_mumu = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
6912 assert_relative_eq!(d.l_mumu, fd_mumu, epsilon = 1e-4);
6913 }
6914
6915 #[test]
6916 fn cloglog_ghq_derivatives_finite_difference_sigma() {
6917 let ctx = QuadratureContext::new();
6919 let mu = 0.2;
6920 let sigma = 1.0;
6921 let h = 1e-6;
6922 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
6923 let l_plus = cloglog_ghq_value(&ctx, mu, sigma + h, 31);
6924 let l_minus = cloglog_ghq_value(&ctx, mu, sigma - h, 31);
6925 let fd_sigma = (l_plus - l_minus) / (2.0 * h);
6926 assert_relative_eq!(d.l_sigma, fd_sigma, epsilon = 1e-5);
6927 }
6928
6929 #[test]
6930 fn cloglog_ghq_derivatives_finite_difference_cross() {
6931 let ctx = QuadratureContext::new();
6933 let mu = -0.5;
6934 let sigma = 0.6;
6935 let h = 1e-6;
6936 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
6937 let d_plus = cloglog_ghq_derivatives(&ctx, mu, sigma + h, 31);
6938 let d_minus = cloglog_ghq_derivatives(&ctx, mu, sigma - h, 31);
6939 let fd_musigma = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
6940 assert_relative_eq!(d.l_musigma, fd_musigma, epsilon = 1e-4);
6941 }
6942
6943 #[test]
6944 fn cloglog_ghq_l_mu_nonnegative() {
6945 let ctx = QuadratureContext::new();
6947 for &mu in &[-3.0, -1.0, 0.0, 1.0, 3.0] {
6948 for &sigma in &[0.1, 0.5, 1.0, 2.0] {
6949 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 21);
6950 assert!(
6951 d.l_mu >= -1e-14,
6952 "L_mu should be non-negative at mu={mu}, sigma={sigma}: got {}",
6953 d.l_mu
6954 );
6955 }
6956 }
6957 }
6958
6959 #[test]
6960 fn cloglog_ghq_adaptive_matches_explicit() {
6961 let ctx = QuadratureContext::new();
6962 let mu = 0.7;
6963 let sigma = 1.2;
6964 let adaptive = cloglog_ghq_derivatives_adaptive(&ctx, mu, sigma);
6965 let n = adaptive_point_count_from_sd(sigma);
6966 let explicit = cloglog_ghq_derivatives(&ctx, mu, sigma, n);
6967 assert_relative_eq!(adaptive.l, explicit.l, epsilon = 1e-15);
6968 assert_relative_eq!(adaptive.l_mu, explicit.l_mu, epsilon = 1e-15);
6969 assert_relative_eq!(adaptive.l_sigma, explicit.l_sigma, epsilon = 1e-15);
6970 assert_relative_eq!(adaptive.l_mumu, explicit.l_mumu, epsilon = 1e-15);
6971 }
6972
6973 #[test]
6974 fn cloglog_ghq_value_matches_mathematical_target_in_central_regime() {
6975 let ctx = QuadratureContext::new();
6976 for &mu in &[-1.0, 0.0, 0.5, 2.0] {
6977 for &sigma in &[0.1, 0.5, 1.0] {
6978 let ghq = cloglog_ghq_value(&ctx, mu, sigma, 51);
6979 let (expected_mean, _) = cloglog_reference_mean_and_derivative(mu, sigma);
6980 assert_relative_eq!(ghq, expected_mean, epsilon = 1e-12, max_relative = 2e-8);
6981 }
6982 }
6983 }
6984
6985 #[test]
6988 fn cloglog_negative_tail_mean_matches_exact_near_transition() {
6989 let eta: f64 = -30.0;
6993 let exact = {
6994 let ex = eta.exp();
6995 -(-ex).exp_m1()
6996 };
6997 let tail = cloglog_negative_tail_mean(eta);
6998 assert!(
6999 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7000 "tail mean at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7001 );
7002 }
7003
7004 #[inline]
7005 fn cloglog_negative_tail_derivative(eta: f64) -> f64 {
7006 if eta < -745.0 {
7008 0.0
7009 } else {
7010 let ex = safe_exp(eta);
7011 (ex * (-ex).exp()).max(0.0)
7012 }
7013 }
7014
7015 #[test]
7016 fn cloglog_negative_tail_derivative_matches_exact_near_transition() {
7017 let eta: f64 = -30.0;
7019 let ex = eta.exp();
7020 let exact = ex * (-ex).exp();
7021 let tail = cloglog_negative_tail_derivative(eta);
7022 assert!(
7023 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7024 "tail derivative at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7025 );
7026 }
7027
7028 #[test]
7029 fn cloglog_negative_tail_degenerate_branch_matches_target_near_transition() {
7030 let ctx = QuadratureContext::default();
7031 let sigma = 0.0;
7032 for &mu in &[-30.001, -30.0, -29.999] {
7033 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7034 assert_relative_eq!(
7035 out.mean,
7036 cloglog_mean_exact(mu),
7037 epsilon = 1e-28,
7038 max_relative = 1e-15
7039 );
7040 assert_relative_eq!(
7041 out.dmean_dmu,
7042 cloglog_mean_d1_exact(mu),
7043 epsilon = 1e-28,
7044 max_relative = 1e-15
7045 );
7046 }
7047 }
7048
7049 #[test]
7050 fn cloglog_negative_tail_small_sigma_branch_matches_target_near_transition() {
7051 let ctx = QuadratureContext::default();
7052 let sigma = 0.1;
7053 for &mu in &[-30.001, -30.0, -29.999] {
7054 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7055 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
7056 assert_relative_eq!(
7057 out.mean,
7058 expected_mean,
7059 epsilon = 1e-24,
7060 max_relative = 1e-10
7061 );
7062 assert_relative_eq!(
7063 out.dmean_dmu,
7064 expected_deriv,
7065 epsilon = 1e-24,
7066 max_relative = 1e-10
7067 );
7068 }
7069 }
7070
7071 fn ref_cholesky_heap(cov: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
7075 let n = cov.len();
7076 if n == 0 || cov.iter().any(|r| r.len() != n) {
7077 return None;
7078 }
7079 let mut base = cov.to_vec();
7080 for retry in 0..8 {
7081 let jitter = if retry == 0 {
7082 0.0
7083 } else {
7084 1e-12 * 10f64.powi(retry - 1)
7085 };
7086 if jitter > 0.0 {
7087 for i in 0..n {
7088 base[i][i] = cov[i][i] + jitter;
7089 }
7090 }
7091 let mut l = vec![vec![0.0_f64; n]; n];
7092 let mut ok = true;
7093 for i in 0..n {
7094 for j in 0..=i {
7095 let mut sum = base[i][j];
7096 for k in 0..j {
7097 sum -= l[i][k] * l[j][k];
7098 }
7099 if i == j {
7100 if !sum.is_finite() || sum <= 0.0 {
7101 ok = false;
7102 break;
7103 }
7104 l[i][j] = sum.sqrt();
7105 } else {
7106 l[i][j] = sum / l[j][j];
7107 }
7108 }
7109 if !ok {
7110 break;
7111 }
7112 }
7113 if ok {
7114 return Some(l);
7115 }
7116 }
7117 None
7118 }
7119
7120 #[test]
7121 fn cholesky_static_matches_heap_d2() {
7122 let cases: &[[[f64; 2]; 2]] = &[
7125 [[1.0, 0.0], [0.0, 1.0]],
7126 [[2.5, 0.3], [0.3, 0.75]],
7127 [[1.0, 0.9999], [0.9999, 1.0]],
7128 [[1e-10, 0.0], [0.0, 1e-10]],
7129 [[4.0, -1.5], [-1.5, 2.25]],
7130 ];
7131 for cov in cases {
7132 let stack = cholesky_static_with_jitter::<2>(cov).expect("stack cholesky");
7133 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7134 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7135 for i in 0..2 {
7136 for j in 0..2 {
7137 assert_eq!(
7138 stack[i][j].to_bits(),
7139 heap[i][j].to_bits(),
7140 "mismatch at ({i},{j}) for cov={cov:?}"
7141 );
7142 }
7143 }
7144 }
7145 }
7146
7147 #[test]
7148 fn cholesky_static_matches_heap_d3() {
7149 let cases: &[[[f64; 3]; 3]] = &[
7150 [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
7151 [[2.0, 0.5, 0.1], [0.5, 1.5, -0.2], [0.1, -0.2, 0.8]],
7152 [[4.0, 1.0, 0.5], [1.0, 3.0, 0.25], [0.5, 0.25, 2.0]],
7153 ];
7154 for cov in cases {
7155 let stack = cholesky_static_with_jitter::<3>(cov).expect("stack cholesky");
7156 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7157 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7158 for i in 0..3 {
7159 for j in 0..3 {
7160 assert_eq!(
7161 stack[i][j].to_bits(),
7162 heap[i][j].to_bits(),
7163 "mismatch at ({i},{j}) for cov={cov:?}"
7164 );
7165 }
7166 }
7167 }
7168 }
7169
7170 #[test]
7171 fn cholesky_static_d1() {
7172 let l = cholesky_static_with_jitter::<1>(&[[2.25]]).expect("d=1");
7173 assert_eq!(l[0][0], 1.5);
7174 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e-13]]).is_some());
7184 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e3]]).is_none());
7187 }
7188}