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]
2888pub fn integrated_logit_inverse_link_jet_pirls(
2889 quadctx: &QuadratureContext,
2890 mu: f64,
2891 sigma: f64,
2892) -> Result<IntegratedInverseLinkJet, EstimationError> {
2893 if sigma <= 1e-10 {
2898 let (mean, d1, d2, d3) = component_point_jet(LinkComponent::Logit, mu);
2899 return Ok(IntegratedInverseLinkJet {
2900 mean,
2901 d1,
2902 d2,
2903 d3,
2904 mode: IntegratedExpectationMode::ExactClosedForm,
2905 });
2906 }
2907 if sigma > LOGIT_JET_GHQ_SIGMA_MAX {
2908 return logit_wide_sigma_jet(mu, sigma);
2909 }
2910 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
2911 component_point_jet(LinkComponent::Logit, x)
2912 });
2913 let mode = match logit_posterior_meanwith_deriv_controlled(mu, sigma) {
2914 Ok(scalar) => scalar.mode,
2915 Err(_) => IntegratedExpectationMode::QuadratureFallback,
2916 };
2917 Ok(IntegratedInverseLinkJet {
2918 mean,
2919 d1: d1.max(0.0),
2920 d2,
2921 d3,
2922 mode,
2923 })
2924}
2925
2926#[inline]
2927fn sas_point_jet(x: f64, epsilon: f64, log_delta: f64) -> (f64, f64, f64, f64) {
2928 let jet = sas_inverse_link_jet(x, epsilon, log_delta)
2929 .expect("normal quadrature nodes must be finite");
2930 (jet.mu, jet.d1, jet.d2, jet.d3)
2931}
2932
2933#[inline]
2934fn beta_logistic_point_jet(x: f64, log_shape_center: f64, epsilon: f64) -> (f64, f64, f64, f64) {
2935 let jet = beta_logistic_inverse_link_jet(x, log_shape_center, epsilon);
2936 (jet.mu, jet.d1, jet.d2, jet.d3)
2937}
2938
2939#[inline]
2940fn worse_integrated_expectation_mode(
2941 lhs: IntegratedExpectationMode,
2942 rhs: IntegratedExpectationMode,
2943) -> IntegratedExpectationMode {
2944 if lhs.rank() >= rhs.rank() { lhs } else { rhs }
2945}
2946
2947#[inline]
2948fn integrated_scalar_drift_exceeds(
2949 candidate: f64,
2950 reference: f64,
2951 abs_tol: f64,
2952 rel_tol: f64,
2953) -> bool {
2954 if !(candidate.is_finite() && reference.is_finite()) {
2955 return true;
2956 }
2957 (candidate - reference).abs() > abs_tol.max(rel_tol * reference.abs().max(candidate.abs()))
2958}
2959
2960#[inline]
2961fn integrated_mean_derivative_drift_exceeds(
2962 candidate: &IntegratedMeanDerivative,
2963 reference: &IntegratedMeanDerivative,
2964 mean_abs_tol: f64,
2965 mean_rel_tol: f64,
2966 deriv_abs_tol: f64,
2967 deriv_rel_tol: f64,
2968) -> bool {
2969 integrated_scalar_drift_exceeds(candidate.mean, reference.mean, mean_abs_tol, mean_rel_tol)
2970 || integrated_scalar_drift_exceeds(
2971 candidate.dmean_dmu,
2972 reference.dmean_dmu,
2973 deriv_abs_tol,
2974 deriv_rel_tol,
2975 )
2976}
2977
2978#[inline]
2979fn component_point_jet(component: LinkComponent, x: f64) -> (f64, f64, f64, f64) {
2980 let jet = component_inverse_link_jet(component, x);
2983 (jet.mu, jet.d1, jet.d2, jet.d3)
2984}
2985
2986#[inline]
2987fn integrated_mixture_component_jet(
2988 ctx: &QuadratureContext,
2989 component: LinkComponent,
2990 mu: f64,
2991 sigma: f64,
2992) -> IntegratedInverseLinkJet {
2993 match component {
2998 LinkComponent::Logit => integrated_inverse_link_jet(ctx, LinkFunction::Logit, mu, sigma)
2999 .unwrap_or_else(|_| integrated_logit_jet_ghq(ctx, mu, sigma)),
3000 LinkComponent::Probit => integrated_probit_jet(mu, sigma),
3001 LinkComponent::CLogLog => integrated_cloglog_inverse_link_jet_controlled(ctx, mu, sigma),
3002 LinkComponent::LogLog | LinkComponent::Cauchit => {
3003 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3004 component_point_jet(component, x)
3005 });
3006 IntegratedInverseLinkJet {
3007 mean,
3008 d1: d1.max(0.0),
3009 d2,
3010 d3,
3011 mode: if sigma <= 1e-10 {
3012 IntegratedExpectationMode::ExactClosedForm
3013 } else {
3014 IntegratedExpectationMode::QuadratureFallback
3015 },
3016 }
3017 }
3018 }
3019}
3020
3021#[inline]
3022fn integrated_mixture_jet(
3023 ctx: &QuadratureContext,
3024 mu: f64,
3025 sigma: f64,
3026 mixture_state: &MixtureLinkState,
3027) -> Result<IntegratedInverseLinkJet, EstimationError> {
3028 if mixture_state.components.is_empty() {
3033 crate::bail_invalid_estim!(
3034 "integrated mixture-link jet requires at least one blended component"
3035 );
3036 }
3037 if mixture_state.components.len() != mixture_state.pi.len() {
3038 crate::bail_invalid_estim!(
3039 "integrated mixture-link jet requires matching component and weight counts"
3040 );
3041 }
3042
3043 let mut mean = 0.0_f64;
3048 let mut d1 = 0.0_f64;
3049 let mut d2 = 0.0_f64;
3050 let mut d3 = 0.0_f64;
3051 let mut mode = IntegratedExpectationMode::ExactClosedForm;
3052 let mut saw_positive_weight = false;
3053
3054 for (&component, &weight) in mixture_state.components.iter().zip(mixture_state.pi.iter()) {
3055 if weight <= 0.0 {
3056 continue;
3057 }
3058 let jet = integrated_mixture_component_jet(ctx, component, mu, sigma);
3059 mean += weight * jet.mean;
3060 d1 += weight * jet.d1;
3061 d2 += weight * jet.d2;
3062 d3 += weight * jet.d3;
3063 if jet.mode.rank() > mode.rank() {
3064 mode = jet.mode;
3065 }
3066 saw_positive_weight = true;
3067 }
3068
3069 if !saw_positive_weight {
3070 crate::bail_invalid_estim!(
3071 "integrated mixture-link jet requires at least one positive component weight"
3072 .to_string(),
3073 );
3074 }
3075
3076 Ok(IntegratedInverseLinkJet {
3077 mean,
3078 d1: d1.max(0.0),
3079 d2,
3080 d3,
3081 mode,
3082 })
3083}
3084
3085#[inline]
3086fn integrated_sas_jet_ghq(
3087 ctx: &QuadratureContext,
3088 mu: f64,
3089 sigma: f64,
3090 sas_state: &SasLinkState,
3091) -> IntegratedInverseLinkJet {
3092 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3093 sas_point_jet(x, sas_state.epsilon, sas_state.log_delta)
3094 });
3095 IntegratedInverseLinkJet {
3096 mean,
3097 d1: d1.max(0.0),
3098 d2,
3099 d3,
3100 mode: if sigma <= 1e-10 {
3101 IntegratedExpectationMode::ExactClosedForm
3102 } else {
3103 IntegratedExpectationMode::QuadratureFallback
3104 },
3105 }
3106}
3107
3108#[inline]
3109fn integrated_beta_logistic_jet_ghq(
3110 ctx: &QuadratureContext,
3111 mu: f64,
3112 sigma: f64,
3113 beta_state: &SasLinkState,
3114) -> IntegratedInverseLinkJet {
3115 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3116 beta_logistic_point_jet(x, beta_state.log_delta, beta_state.epsilon)
3117 });
3118 IntegratedInverseLinkJet {
3119 mean,
3120 d1: d1.max(0.0),
3121 d2,
3122 d3,
3123 mode: if sigma <= 1e-10 {
3124 IntegratedExpectationMode::ExactClosedForm
3125 } else {
3126 IntegratedExpectationMode::QuadratureFallback
3127 },
3128 }
3129}
3130
3131#[inline]
3133pub fn integrated_inverse_link_jetwith_state(
3134 quadctx: &QuadratureContext,
3135 link: LinkFunction,
3136 mu: f64,
3137 sigma: f64,
3138 mixture_link_state: Option<&MixtureLinkState>,
3139 sas_link_state: Option<&SasLinkState>,
3140) -> Result<IntegratedInverseLinkJet, EstimationError> {
3141 if let Some(state) = mixture_link_state {
3142 return integrated_mixture_jet(quadctx, mu, sigma, state);
3143 }
3144 if matches!(link, LinkFunction::Sas) {
3145 let sas = sas_link_state.ok_or_else(|| {
3146 EstimationError::InvalidInput(
3147 "state-less integrated SAS jet is unsupported; explicit SasLinkState is required"
3148 .to_string(),
3149 )
3150 })?;
3151 return Ok(integrated_sas_jet_ghq(quadctx, mu, sigma, sas));
3152 }
3153 if matches!(link, LinkFunction::BetaLogistic) {
3154 let state = sas_link_state.ok_or_else(|| {
3155 EstimationError::InvalidInput(
3156 "state-less integrated Beta-Logistic jet is unsupported; explicit link state is required"
3157 .to_string(),
3158 )
3159 })?;
3160 return Ok(integrated_beta_logistic_jet_ghq(quadctx, mu, sigma, state));
3161 }
3162 integrated_inverse_link_jet(quadctx, link, mu, sigma)
3163}
3164
3165#[inline]
3175pub fn integrated_family_moments_jet(
3176 quadctx: &QuadratureContext,
3177 likelihood: &GlmLikelihoodSpec,
3178 eta: f64,
3179 se_eta: f64,
3180) -> Result<IntegratedMomentsJet, EstimationError> {
3181 const PROB_EPS: f64 = 1e-12;
3182 if !(eta.is_finite() && (-700.0..=700.0).contains(&eta)) {
3183 crate::bail_invalid_estim!(
3184 "integrated moments eta must be finite and within [-700, 700]; got {eta}"
3185 );
3186 }
3187 let e = eta;
3188 let se = se_eta.max(0.0);
3189 let resolved_scale = likelihood
3193 .resolved_scale()
3194 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3195 let spec = &likelihood.spec;
3196 let mixture_link_state: Option<&MixtureLinkState> = spec.link.mixture_state();
3197 let sas_link_state: Option<&SasLinkState> = spec.link.sas_state();
3198 match &spec.response {
3199 ResponseFamily::Binomial => match &spec.link {
3200 InverseLink::Standard(StandardLink::Logit) => {
3201 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3202 let mean = jet.mean;
3203 Ok(IntegratedMomentsJet {
3204 mean,
3205 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3206 d1: jet.d1,
3207 d2: jet.d2,
3208 d3: jet.d3,
3209 mode: jet.mode,
3210 })
3211 }
3212 InverseLink::Standard(StandardLink::Probit) => {
3213 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Probit, e, se)?;
3214 let mean = jet.mean;
3215 Ok(IntegratedMomentsJet {
3216 mean,
3217 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3218 d1: jet.d1,
3219 d2: jet.d2,
3220 d3: jet.d3,
3221 mode: jet.mode,
3222 })
3223 }
3224 InverseLink::Standard(StandardLink::CLogLog) => {
3225 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::CLogLog, e, se)?;
3226 let mean = jet.mean;
3227 Ok(IntegratedMomentsJet {
3228 mean,
3229 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3230 d1: jet.d1,
3231 d2: jet.d2,
3232 d3: jet.d3,
3233 mode: jet.mode,
3234 })
3235 }
3236 InverseLink::LatentCLogLog(_) => Err(EstimationError::InvalidInput(
3237 "Binomial+LatentCLogLog integrated moments require an explicit latent cloglog inverse-link state"
3238 .to_string(),
3239 )),
3240 InverseLink::Sas(_) => {
3241 let jet = integrated_inverse_link_jetwith_state(
3242 quadctx,
3243 LinkFunction::Sas,
3244 e,
3245 se,
3246 mixture_link_state,
3247 sas_link_state,
3248 )?;
3249 let mean = jet.mean;
3250 Ok(IntegratedMomentsJet {
3251 mean,
3252 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3253 d1: jet.d1,
3254 d2: jet.d2,
3255 d3: jet.d3,
3256 mode: jet.mode,
3257 })
3258 }
3259 InverseLink::BetaLogistic(_) => {
3260 let jet = integrated_inverse_link_jetwith_state(
3261 quadctx,
3262 LinkFunction::BetaLogistic,
3263 e,
3264 se,
3265 mixture_link_state,
3266 sas_link_state,
3267 )?;
3268 let mean = jet.mean;
3269 Ok(IntegratedMomentsJet {
3270 mean,
3271 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3272 d1: jet.d1,
3273 d2: jet.d2,
3274 d3: jet.d3,
3275 mode: jet.mode,
3276 })
3277 }
3278 InverseLink::Mixture(state) => {
3279 let jet = integrated_mixture_jet(quadctx, e, se, &state)?;
3280 let mean = jet.mean;
3281 Ok(IntegratedMomentsJet {
3282 mean,
3283 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3284 d1: jet.d1,
3285 d2: jet.d2,
3286 d3: jet.d3,
3287 mode: jet.mode,
3288 })
3289 }
3290 InverseLink::Standard(other) => Err(EstimationError::InvalidInput(format!(
3291 "Binomial response paired with unsupported standard link {other:?} for integrated moments"
3292 ))),
3293 },
3294 ResponseFamily::Gaussian => {
3295 let variance = resolved_scale
3296 .gaussian_phi()
3297 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3298 Ok(IntegratedMomentsJet {
3299 mean: e,
3300 variance,
3301 d1: 1.0,
3302 d2: 0.0,
3303 d3: 0.0,
3304 mode: IntegratedExpectationMode::ExactClosedForm,
3305 })
3306 }
3307 ResponseFamily::RoystonParmar => {
3308 let jet = integrated_inverse_link_jetwith_state(
3309 quadctx,
3310 LinkFunction::CLogLog,
3311 e,
3312 se,
3313 mixture_link_state,
3314 sas_link_state,
3315 )?;
3316 let mean = (1.0 - jet.mean).clamp(0.0, 1.0);
3317 Ok(IntegratedMomentsJet {
3318 mean,
3319 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3320 d1: -jet.d1,
3321 d2: -jet.d2,
3322 d3: -jet.d3,
3323 mode: jet.mode,
3324 })
3325 }
3326 ResponseFamily::Beta { .. } => {
3327 let precision = resolved_scale
3328 .beta_precision()
3329 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3330 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3331 let mean = jet.mean.clamp(PROB_EPS, 1.0 - PROB_EPS);
3332 Ok(IntegratedMomentsJet {
3333 mean,
3334 variance: (mean * (1.0 - mean) / (1.0 + precision)).max(PROB_EPS),
3335 d1: jet.d1,
3336 d2: jet.d2,
3337 d3: jet.d3,
3338 mode: jet.mode,
3339 })
3340 }
3341 ResponseFamily::Poisson
3342 | ResponseFamily::Tweedie { .. }
3343 | ResponseFamily::NegativeBinomial { .. }
3344 | ResponseFamily::Gamma => {
3345 let s2 = se * se;
3350 let (mean, saturated) = safe_expwith_saturation(e + 0.5 * s2);
3351 let variance = match &spec.response {
3362 ResponseFamily::Poisson => mean,
3363 ResponseFamily::Tweedie { p } => {
3364 let phi = resolved_scale
3365 .tweedie_phi()
3366 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3367 phi * mean.powf(*p)
3368 }
3369 ResponseFamily::NegativeBinomial { .. } => {
3370 let theta = resolved_scale.negative_binomial_theta().map_err(|error| {
3371 EstimationError::InvalidInput(error.to_string())
3372 })?;
3373 mean + mean * mean / theta
3374 }
3375 ResponseFamily::Gamma => {
3376 let phi = resolved_scale
3377 .gamma_phi()
3378 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3379 phi * mean * mean
3380 }
3381 other => {
3385 return Err(EstimationError::InvalidInput(format!(
3386 "integrated log-normal moments reached unexpected family {other:?}"
3387 )));
3388 }
3389 };
3390 if !(variance.is_finite() && variance >= 0.0) {
3391 return Err(EstimationError::InvalidInput(format!(
3392 "integrated {} variance is not representable: {variance:?}",
3393 spec.response.name()
3394 )));
3395 }
3396 Ok(IntegratedMomentsJet {
3397 mean,
3398 variance,
3399 d1: mean,
3400 d2: mean,
3401 d3: mean,
3402 mode: if saturated {
3403 IntegratedExpectationMode::ControlledAsymptotic
3404 } else {
3405 IntegratedExpectationMode::ExactClosedForm
3406 },
3407 })
3408 }
3409 }
3410}
3411
3412pub fn logit_posterior_meanwith_deriv_batch(
3415 ctx: &QuadratureContext,
3416 eta: &ndarray::Array1<f64>,
3417 se_eta: &ndarray::Array1<f64>,
3418) -> Result<(ndarray::Array1<f64>, ndarray::Array1<f64>), EstimationError> {
3419 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3420 let n = eta.len();
3421 let pairs: Result<Vec<(f64, f64)>, _> = (0..n)
3423 .into_par_iter()
3424 .map(|i| {
3425 let integrated = integrated_inverse_link_mean_and_derivative(
3426 ctx,
3427 LinkFunction::Logit,
3428 eta[i],
3429 se_eta[i],
3430 )?;
3431 Ok::<_, EstimationError>((integrated.mean, integrated.dmean_dmu))
3432 })
3433 .collect();
3434 let pairs = pairs?;
3435 let mut mu = ndarray::Array1::<f64>::zeros(n);
3436 let mut dmu = ndarray::Array1::<f64>::zeros(n);
3437 for (i, (m, d)) in pairs.into_iter().enumerate() {
3438 mu[i] = m;
3439 dmu[i] = d;
3440 }
3441
3442 Ok((mu, dmu))
3443}
3444
3445pub fn logit_posterior_mean_batch(
3449 ctx: &QuadratureContext,
3450 eta: &ndarray::Array1<f64>,
3451 se_eta: &ndarray::Array1<f64>,
3452) -> Result<ndarray::Array1<f64>, EstimationError> {
3453 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3454 let n = eta.len();
3455 let values: Result<Vec<f64>, EstimationError> = (0..n)
3456 .into_par_iter()
3457 .map(|i| {
3458 integrated_inverse_link_mean_and_derivative(ctx, LinkFunction::Logit, eta[i], se_eta[i])
3459 .map(|integrated| integrated.mean)
3460 })
3461 .collect();
3462 Ok(ndarray::Array1::from_vec(values?))
3463}
3464
3465pub trait GhqValue: Sized {
3466 fn zero() -> Self;
3467 fn addweighted(&mut self, weight: f64, value: Self);
3468 fn scale(self, factor: f64) -> Self;
3469}
3470
3471impl GhqValue for f64 {
3472 #[inline]
3473 fn zero() -> Self {
3474 0.0
3475 }
3476
3477 #[inline]
3478 fn addweighted(&mut self, weight: f64, value: Self) {
3479 *self += weight * value;
3480 }
3481
3482 #[inline]
3483 fn scale(self, factor: f64) -> Self {
3484 self * factor
3485 }
3486}
3487
3488impl GhqValue for (f64, f64) {
3489 #[inline]
3490 fn zero() -> Self {
3491 (0.0, 0.0)
3492 }
3493
3494 #[inline]
3495 fn addweighted(&mut self, weight: f64, value: Self) {
3496 self.0 += weight * value.0;
3497 self.1 += weight * value.1;
3498 }
3499
3500 #[inline]
3501 fn scale(self, factor: f64) -> Self {
3502 (self.0 * factor, self.1 * factor)
3503 }
3504}
3505
3506impl GhqValue for (f64, f64, f64, f64) {
3507 #[inline]
3508 fn zero() -> Self {
3509 (0.0, 0.0, 0.0, 0.0)
3510 }
3511
3512 #[inline]
3513 fn addweighted(&mut self, weight: f64, value: Self) {
3514 self.0 += weight * value.0;
3515 self.1 += weight * value.1;
3516 self.2 += weight * value.2;
3517 self.3 += weight * value.3;
3518 }
3519
3520 #[inline]
3521 fn scale(self, factor: f64) -> Self {
3522 (
3523 self.0 * factor,
3524 self.1 * factor,
3525 self.2 * factor,
3526 self.3 * factor,
3527 )
3528 }
3529}
3530
3531impl GhqValue for (f64, f64, f64, f64, f64, f64) {
3532 #[inline]
3533 fn zero() -> Self {
3534 (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
3535 }
3536
3537 #[inline]
3538 fn addweighted(&mut self, weight: f64, value: Self) {
3539 self.0 += weight * value.0;
3540 self.1 += weight * value.1;
3541 self.2 += weight * value.2;
3542 self.3 += weight * value.3;
3543 self.4 += weight * value.4;
3544 self.5 += weight * value.5;
3545 }
3546
3547 #[inline]
3548 fn scale(self, factor: f64) -> Self {
3549 (
3550 self.0 * factor,
3551 self.1 * factor,
3552 self.2 * factor,
3553 self.3 * factor,
3554 self.4 * factor,
3555 self.5 * factor,
3556 )
3557 }
3558}
3559
3560#[inline]
3561fn integrate_normal_ghq_adaptive<F, R>(ctx: &QuadratureContext, eta: f64, se_eta: f64, f: F) -> R
3562where
3563 F: Fn(f64) -> R,
3564 R: GhqValue,
3565{
3566 if se_eta < 1e-10 {
3567 return f(eta);
3568 }
3569 let n = adaptive_point_count_from_sd(se_eta.abs());
3570 with_gh_nodesweights(ctx, n, |nodes, weights| {
3571 let scale = SQRT_2 * se_eta;
3572 let mut sum = R::zero();
3573 for i in 0..n {
3574 sum.addweighted(weights[i], f(eta + scale * nodes[i]));
3575 }
3576 sum.scale(1.0 / std::f64::consts::PI.sqrt())
3577 })
3578}
3579
3580#[inline]
3581fn integrated_probit_jet(mu: f64, sigma: f64) -> IntegratedInverseLinkJet {
3582 let s = sigma.hypot(1.0);
3588 let z = mu / s;
3589 let mean = gam_math::probability::normal_cdf(z);
3590 let pdf = gam_math::probability::normal_pdf(z);
3591 if pdf == 0.0 {
3592 return IntegratedInverseLinkJet {
3593 mean,
3594 d1: 0.0,
3595 d2: 0.0,
3596 d3: 0.0,
3597 mode: IntegratedExpectationMode::ExactClosedForm,
3598 };
3599 }
3600 IntegratedInverseLinkJet {
3601 mean,
3602 d1: pdf / s,
3603 d2: -z * pdf / (s * s),
3604 d3: (z * z - 1.0) * pdf / (s * s * s),
3605 mode: IntegratedExpectationMode::ExactClosedForm,
3606 }
3607}
3608
3609#[inline]
3610fn integrated_logit_jet_ghq(
3611 ctx: &QuadratureContext,
3612 mu: f64,
3613 sigma: f64,
3614) -> IntegratedInverseLinkJet {
3615 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3616 component_point_jet(LinkComponent::Logit, x)
3617 });
3618 IntegratedInverseLinkJet {
3619 mean,
3620 d1: d1.max(0.0),
3621 d2,
3622 d3,
3623 mode: if sigma <= 1e-10 {
3624 IntegratedExpectationMode::ExactClosedForm
3625 } else {
3626 IntegratedExpectationMode::QuadratureFallback
3627 },
3628 }
3629}
3630
3631#[inline]
3632fn cloglog_inverse_link_controlled_values(
3633 ctx: &QuadratureContext,
3634 mu: f64,
3635 sigma: f64,
3636 max_order: usize,
3637) -> ([f64; 6], IntegratedExpectationMode) {
3638 assert!(max_order <= 5);
3639 if sigma <= 1e-10 {
3640 let (mean, d1, d2, d3, d4, d5) = cloglog_point_jet5(mu);
3641 return (
3642 [mean, d1, d2, d3, d4, d5],
3643 IntegratedExpectationMode::ExactClosedForm,
3644 );
3645 }
3646
3647 let (k, log_k0, mode) = latent_cloglog_kernel_terms(ctx, mu, sigma, max_order);
3648 let mut values = [0.0; 6];
3649 values[0] = if log_k0.is_finite() {
3650 -log_k0.exp_m1()
3651 } else {
3652 1.0
3653 };
3654 values[1] = k[1].max(0.0);
3655 if sigma > CLOGLOG_JET_MOMENT_SIGMA_MAX {
3656 if max_order >= 2 {
3657 values[2] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).2);
3658 }
3659 if max_order >= 3 {
3660 values[3] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).3);
3661 }
3662 if max_order >= 4 {
3663 values[4] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).4);
3664 }
3665 if max_order >= 5 {
3666 values[5] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).5);
3667 }
3668 return (
3669 values,
3670 worse_integrated_expectation_mode(mode, IntegratedExpectationMode::QuadratureFallback),
3671 );
3672 }
3673 if max_order >= 2 {
3674 values[2] = k[1] - k[2];
3675 }
3676 if max_order >= 3 {
3677 values[3] = k[1] - 3.0 * k[2] + k[3];
3678 }
3679 if max_order >= 4 {
3680 values[4] = k[1] - 7.0 * k[2] + 6.0 * k[3] - k[4];
3681 }
3682 if max_order >= 5 {
3683 values[5] = k[1] - 15.0 * k[2] + 25.0 * k[3] - 10.0 * k[4] + k[5];
3684 }
3685 (values, mode)
3686}
3687
3688#[inline]
3689pub(crate) fn latent_cloglog_inverse_link_jet5_controlled(
3690 ctx: &QuadratureContext,
3691 mu: f64,
3692 sigma: f64,
3693) -> IntegratedInverseLinkJet5 {
3694 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 5);
3695 IntegratedInverseLinkJet5 {
3696 mean: values[0],
3697 d1: values[1],
3698 d2: values[2],
3699 d3: values[3],
3700 d4: values[4],
3701 d5: values[5],
3702 mode,
3703 }
3704}
3705
3706#[derive(Clone, Copy, Debug)]
3716pub struct LatentCLogLogJet5 {
3717 pub mean: f64,
3718 pub d1: f64,
3719 pub d2: f64,
3720 pub d3: f64,
3721 pub d4: f64,
3722 pub d5: f64,
3723 pub mode: IntegratedExpectationMode,
3724}
3725
3726pub fn latent_cloglog_jet5(
3727 quadctx: &QuadratureContext,
3728 eta: f64,
3729 sigma: f64,
3730) -> Result<LatentCLogLogJet5, EstimationError> {
3731 validate_latent_cloglog_inputs(eta, sigma)?;
3732 let jet = latent_cloglog_inverse_link_jet5_controlled(quadctx, eta, sigma);
3738 Ok(LatentCLogLogJet5 {
3739 mean: jet.mean,
3740 d1: jet.d1,
3741 d2: jet.d2,
3742 d3: jet.d3,
3743 d4: jet.d4,
3744 d5: jet.d5,
3745 mode: jet.mode,
3746 })
3747}
3748
3749#[inline]
3750pub fn latent_cloglog_inverse_link_jet(
3751 quadctx: &QuadratureContext,
3752 eta: f64,
3753 sigma: f64,
3754) -> Result<IntegratedInverseLinkJet, EstimationError> {
3755 let jet = latent_cloglog_jet5(quadctx, eta, sigma)?;
3756 Ok(IntegratedInverseLinkJet {
3757 mean: jet.mean,
3758 d1: jet.d1,
3759 d2: jet.d2,
3760 d3: jet.d3,
3761 mode: jet.mode,
3762 })
3763}
3764
3765#[inline]
3766fn integrated_cloglog_inverse_link_jet_controlled(
3767 ctx: &QuadratureContext,
3768 mu: f64,
3769 sigma: f64,
3770) -> IntegratedInverseLinkJet {
3771 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 3);
3772 IntegratedInverseLinkJet {
3773 mean: values[0],
3774 d1: values[1],
3775 d2: values[2],
3776 d3: values[3],
3777 mode,
3778 }
3779}
3780
3781#[inline]
3782fn latent_cloglog_kernel_terms(
3783 ctx: &QuadratureContext,
3784 mu: f64,
3785 sigma: f64,
3786 max_order: usize,
3787) -> ([f64; 6], f64, IntegratedExpectationMode) {
3788 let sigma2 = sigma * sigma;
3789 let mut k = [0.0; 6];
3790 let mut log_k0 = f64::NEG_INFINITY;
3791 let mut mode = IntegratedExpectationMode::ExactClosedForm;
3792
3793 for (order, out) in k.iter_mut().enumerate().take(max_order + 1) {
3794 let kf = order as f64;
3795 let shifted_mu = mu + kf * sigma2;
3796 let (log_survival, term_mode) =
3804 cloglog_log_survival_term_controlled(ctx, shifted_mu, sigma);
3805 mode = worse_integrated_expectation_mode(mode, term_mode);
3806
3807 let log_value = kf * mu + 0.5 * kf * kf * sigma2 + log_survival;
3808 if order == 0 {
3809 log_k0 = log_value;
3810 }
3811 if !log_value.is_finite() {
3812 *out = 0.0;
3813 continue;
3814 }
3815 let upper = if order == 0 {
3816 1.0
3817 } else {
3818 let k_over_e = kf / std::f64::consts::E;
3819 k_over_e.powf(kf)
3820 };
3821 *out = safe_exp(log_value).clamp(0.0, upper);
3822 }
3823
3824 (k, log_k0, mode)
3825}
3826
3827#[inline]
3828pub fn normal_expectation_1d_adaptive<F>(
3829 ctx: &QuadratureContext,
3830 eta: f64,
3831 se_eta: f64,
3832 f: F,
3833) -> f64
3834where
3835 F: Fn(f64) -> f64,
3836{
3837 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
3838}
3839
3840#[inline]
3841pub fn normal_expectation_1d_adaptive_pair<F>(
3842 ctx: &QuadratureContext,
3843 eta: f64,
3844 se_eta: f64,
3845 f: F,
3846) -> (f64, f64)
3847where
3848 F: Fn(f64) -> (f64, f64),
3849{
3850 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
3851}
3852
3853fn adaptive_point_count_from_sd(max_sd: f64) -> usize {
3854 if max_sd.is_finite() && max_sd > 2.5 {
3864 51
3865 } else if max_sd.is_finite() && max_sd > 0.5 {
3866 31
3867 } else if max_sd.is_finite() && max_sd > 0.35 {
3868 21
3869 } else if max_sd.is_finite() && max_sd > 0.1 {
3870 15
3871 } else {
3872 7
3873 }
3874}
3875
3876#[inline]
3877fn with_gh_nodesweights<R>(
3878 ctx: &QuadratureContext,
3879 n: usize,
3880 f: impl FnOnce(&[f64], &[f64]) -> R,
3881) -> R {
3882 if n == 7 {
3883 let gh = ctx.gauss_hermite();
3884 f(&gh.nodes, &gh.weights)
3885 } else {
3886 let gh = ctx.gauss_hermite_n(n);
3887 f(&gh.nodes, &gh.weights)
3888 }
3889}
3890
3891#[inline]
3901fn cholesky_static<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
3902 let mut l = [[0.0_f64; D]; D];
3903 for i in 0..D {
3904 for j in 0..=i {
3905 let mut sum = cov[i][j];
3906 for k in 0..j {
3907 sum -= l[i][k] * l[j][k];
3908 }
3909 if i == j {
3910 if !sum.is_finite() || sum <= 0.0 {
3911 return None;
3912 }
3913 l[i][j] = sum.sqrt();
3914 } else {
3915 l[i][j] = sum / l[j][j];
3916 }
3917 }
3918 }
3919 Some(l)
3920}
3921
3922#[inline]
3925fn cholesky_static_with_jitter<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
3926 if D == 0 {
3927 return None;
3928 }
3929 for retry in 0..8 {
3930 let jitter = if retry == 0 {
3931 0.0
3932 } else {
3933 1e-12 * 10f64.powi(retry - 1)
3934 };
3935 if jitter == 0.0 {
3936 if let Some(l) = cholesky_static::<D>(cov) {
3937 return Some(l);
3938 }
3939 } else {
3940 let mut base = *cov;
3941 for i in 0..D {
3942 base[i][i] = cov[i][i] + jitter;
3943 }
3944 if let Some(l) = cholesky_static::<D>(&base) {
3945 return Some(l);
3946 }
3947 }
3948 }
3949 None
3950}
3951
3952#[inline]
3953fn adaptive_point_countwith_cap(max_sd: f64, max_n: usize) -> usize {
3954 adaptive_point_count_from_sd(max_sd).min(max_n)
3955}
3956
3957#[inline]
3958fn ghq_nd_integrate_try<const D: usize, F, R, E>(
3959 ctx: &QuadratureContext,
3960 mu: [f64; D],
3961 cov: [[f64; D]; D],
3962 max_n: usize,
3963 f: F,
3964) -> Result<Option<R>, E>
3965where
3966 F: Fn([f64; D]) -> Result<R, E>,
3967 R: GhqValue,
3968{
3969 let mut maxvar = 0.0_f64;
3970 for (i, row) in cov.iter().enumerate() {
3971 maxvar = maxvar.max(row[i]).max(0.0);
3972 }
3973 let n = adaptive_point_countwith_cap(maxvar.sqrt(), max_n);
3974
3975 let mut cov_arr = cov;
3980 for i in 0..D {
3981 cov_arr[i][i] = cov_arr[i][i].max(0.0);
3982 }
3983 let Some(l) = cholesky_static_with_jitter::<D>(&cov_arr) else {
3984 return Ok(None);
3985 };
3986 let norm = 1.0 / std::f64::consts::PI.powf(0.5 * D as f64);
3987
3988 with_gh_nodesweights(ctx, n, |nodes, weights| {
3989 let mut acc = R::zero();
3990 let mut idx = [0usize; D];
3991 loop {
3992 let mut z = [0.0_f64; D];
3993 let mut weight = 1.0_f64;
3994 for d in 0..D {
3995 z[d] = SQRT_2 * nodes[idx[d]];
3996 weight *= weights[idx[d]];
3997 }
3998
3999 let mut x = mu;
4000 for row in 0..D {
4001 let mut dot = 0.0_f64;
4002 for (col, zc) in z.iter().enumerate().take(row + 1) {
4003 dot += l[row][col] * *zc;
4004 }
4005 x[row] += dot;
4006 }
4007 acc.addweighted(weight, f(x)?);
4008
4009 let mut carry = true;
4010 for d in (0..D).rev() {
4011 idx[d] += 1;
4012 if idx[d] < n {
4013 carry = false;
4014 break;
4015 }
4016 idx[d] = 0;
4017 }
4018 if carry {
4019 break;
4020 }
4021 }
4022 Ok(Some(acc.scale(norm)))
4023 })
4024}
4025
4026#[inline]
4027fn ghq_nd_integrate<const D: usize, F, R>(
4028 ctx: &QuadratureContext,
4029 mu: [f64; D],
4030 cov: [[f64; D]; D],
4031 max_n: usize,
4032 f: F,
4033) -> Option<R>
4034where
4035 F: Fn([f64; D]) -> R,
4036 R: GhqValue,
4037{
4038 match ghq_nd_integrate_try::<D, _, R, Infallible>(ctx, mu, cov, max_n, |x| Ok(f(x))) {
4039 Ok(v) => v,
4040 Err(e) => match e {},
4041 }
4042}
4043
4044#[inline]
4045fn ghq_nd_integrate_result<const D: usize, F, R, E>(
4046 ctx: &QuadratureContext,
4047 mu: [f64; D],
4048 cov: [[f64; D]; D],
4049 max_n: usize,
4050 f: F,
4051) -> Result<Option<R>, E>
4052where
4053 F: Fn([f64; D]) -> Result<R, E>,
4054 R: GhqValue,
4055{
4056 ghq_nd_integrate_try::<D, _, R, E>(ctx, mu, cov, max_n, f)
4057}
4058
4059pub fn normal_expectation_nd_adaptive<const D: usize, F>(
4061 ctx: &QuadratureContext,
4062 mu: [f64; D],
4063 cov: [[f64; D]; D],
4064 max_n: usize,
4065 f: F,
4066) -> f64
4067where
4068 F: Fn([f64; D]) -> f64,
4069{
4070 match ghq_nd_integrate::<D, _, f64>(ctx, mu, cov, max_n, &f) {
4071 Some(v) => v,
4072 None => f(mu),
4073 }
4074}
4075
4076pub fn normal_expectation_nd_adaptive_result<const D: usize, F, R, E>(
4078 ctx: &QuadratureContext,
4079 mu: [f64; D],
4080 cov: [[f64; D]; D],
4081 max_n: usize,
4082 f: F,
4083) -> Result<R, E>
4084where
4085 F: Fn([f64; D]) -> Result<R, E>,
4086 R: GhqValue,
4087{
4088 match ghq_nd_integrate_result::<D, _, R, E>(ctx, mu, cov, max_n, &f)? {
4089 Some(v) => Ok(v),
4090 None => f(mu),
4091 }
4092}
4093
4094pub fn normal_expectation_2d_adaptive_result<F, E>(
4096 ctx: &QuadratureContext,
4097 mu: [f64; 2],
4098 cov: [[f64; 2]; 2],
4099 f: F,
4100) -> Result<f64, E>
4101where
4102 F: Fn(f64, f64) -> Result<f64, E>,
4103{
4104 normal_expectation_nd_adaptive_result::<2, _, _, E>(ctx, mu, cov, 21, |x| f(x[0], x[1]))
4105}
4106
4107pub fn normal_expectation_3d_adaptive<F>(
4109 ctx: &QuadratureContext,
4110 mu: [f64; 3],
4111 cov: [[f64; 3]; 3],
4112 f: F,
4113) -> f64
4114where
4115 F: Fn(f64, f64, f64) -> f64,
4116{
4117 normal_expectation_nd_adaptive::<3, _>(ctx, mu, cov, 15, |x| f(x[0], x[1], x[2]))
4119}
4120
4121#[inline]
4140pub fn probit_posterior_mean(eta: f64, se_eta: f64) -> f64 {
4141 if se_eta < 1e-10 {
4142 return gam_math::probability::normal_cdf(eta);
4143 }
4144 let denom = (1.0 + se_eta * se_eta).sqrt();
4145 gam_math::probability::normal_cdf(eta / denom)
4146}
4147
4148#[inline]
4149pub fn logit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4150 let (m1, m2) = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4151 let p = sigmoid(x);
4152 (p, p * p)
4153 });
4154 let m1 = m1.clamp(0.0, 1.0);
4155 let m2 = m2.clamp(0.0, 1.0);
4156 (m1, (m2 - m1 * m1).max(0.0))
4157}
4158
4159#[inline]
4160pub fn probit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4161 let m1 = probit_posterior_mean(eta, se_eta);
4162 let m2 = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4163 let p = gam_math::probability::normal_cdf(x);
4164 p * p
4165 })
4166 .clamp(0.0, 1.0);
4167 (m1, (m2 - m1 * m1).max(0.0))
4168}
4169
4170#[inline]
4171pub fn cloglog_posterior_meanvariance(
4172 ctx: &QuadratureContext,
4173 eta: f64,
4174 se_eta: f64,
4175) -> (f64, f64) {
4176 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4196 return (cloglog_mean_exact(eta), 0.0);
4197 }
4198 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4199 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4200 let mean = cloglog_mean_from_survival(survival);
4201 let variance = (survival_sq - survival * survival).max(0.0);
4202 (mean, variance)
4203}
4204
4205#[inline]
4239pub fn cloglog_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4240 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4244 return cloglog_mean_exact(eta);
4245 }
4246 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4247 cloglog_mean_from_survival(survival)
4248}
4249
4250#[inline]
4264pub fn survival_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4265 cloglog_survival_term_controlled(ctx, eta, se_eta)
4266 .0
4267 .clamp(0.0, 1.0)
4268}
4269
4270#[inline]
4271pub fn survival_posterior_meanvariance(
4272 ctx: &QuadratureContext,
4273 eta: f64,
4274 se_eta: f64,
4275) -> (f64, f64) {
4276 let (m1, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4277 let (m2, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4278 (m1.clamp(0.0, 1.0), (m2 - m1 * m1).max(0.0))
4279}
4280
4281pub fn logit_posterior_mean_exact(mu: f64, sigma: f64) -> f64 {
4357 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 0.0 {
4358 return sigmoid(mu);
4359 }
4360 if sigma < LOGIT_SIGMA_DEGENERATE {
4361 return sigmoid(mu);
4364 }
4365
4366 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI; let sqrt2_sigma = SQRT_2 * sigma;
4368 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;
4372
4373 let mut corr = 0.0_f64;
4379 let mut n = 1usize;
4380 let tail_start = loop {
4381 let b = (2.0 * (n as f64) - 1.0) * beta;
4382 let abs_xi2 = c * c + b * b;
4383 if abs_xi2 > r2 && n >= FADDEEVA_TAIL_MIN_INDEX {
4384 break n;
4385 }
4386 let xi = Complex { re: c, im: b };
4387 let d = if abs_xi2 > r2 {
4388 inv_sqrt_pi * faddeeva_asymptotic_a(xi).re
4390 } else {
4391 faddeeva_upper_halfplane(xi).im - inv_sqrt_pi * c / abs_xi2
4392 };
4393 corr += d;
4394 n += 1;
4395 };
4396
4397 corr += faddeeva_pole_series_em_tail(c, beta, tail_start, inv_sqrt_pi);
4398
4399 sigmoid(mu) - coeff * corr
4400}
4401
4402const FADDEEVA_TAIL_MIN_INDEX: usize = 48;
4406const FADDEEVA_ASYMPTOTIC_RADIUS: f64 = 7.0;
4409const FADDEEVA_ASYMPTOTIC_TERMS: usize = 14;
4412
4413fn faddeeva_asymptotic_a(xi: Complex) -> Complex {
4417 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4418 let inv2 = complexmul(inv, inv);
4419 let mut xp = complexmul(inv2, inv); let mut cm = 0.5_f64; let mut s = Complex::default();
4422 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4423 s = complex_add(
4424 s,
4425 Complex {
4426 re: cm * xp.re,
4427 im: cm * xp.im,
4428 },
4429 );
4430 cm *= (2.0 * (m as f64) + 1.0) / 2.0; xp = complexmul(xp, inv2);
4432 }
4433 s
4434}
4435
4436fn faddeeva_pole_series_em_tail(c: f64, beta: f64, tail_start: usize, inv_sqrt_pi: f64) -> f64 {
4445 let b_a = (2.0 * (tail_start as f64) - 1.0) * beta;
4446 let xi = Complex { re: c, im: b_a };
4447 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4448 let inv2 = complexmul(inv, inv);
4449 let two_i_beta = Complex {
4451 re: 0.0,
4452 im: 2.0 * beta,
4453 };
4454
4455 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;
4463 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4464 let mf = m as f64;
4465 let inv_4ibm = Complex {
4467 re: 0.0,
4468 im: -1.0 / (4.0 * beta * mf),
4469 };
4470 s = complex_add(
4471 s,
4472 complexmul(
4473 Complex {
4474 re: cm * x2m.re,
4475 im: cm * x2m.im,
4476 },
4477 inv_4ibm,
4478 ),
4479 );
4480 a_acc = complex_add(
4481 a_acc,
4482 Complex {
4483 re: cm * x2m1.re,
4484 im: cm * x2m1.im,
4485 },
4486 );
4487 let fc = cm * (-(2.0 * mf + 1.0));
4488 fp_inner = complex_add(
4489 fp_inner,
4490 Complex {
4491 re: fc * x2m2.re,
4492 im: fc * x2m2.im,
4493 },
4494 );
4495 cm *= (2.0 * mf + 1.0) / 2.0;
4496 x2m = complexmul(x2m, inv2);
4497 x2m1 = complexmul(x2m1, inv2);
4498 x2m2 = complexmul(x2m2, inv2);
4499 }
4500
4501 s = complex_add(
4503 s,
4504 Complex {
4505 re: 0.5 * a_acc.re,
4506 im: 0.5 * a_acc.im,
4507 },
4508 );
4509 let fprime = complexmul(two_i_beta, fp_inner);
4511 s = complex_add(
4512 s,
4513 Complex {
4514 re: -fprime.re / 12.0,
4515 im: -fprime.im / 12.0,
4516 },
4517 );
4518
4519 inv_sqrt_pi * s.re
4522}
4523
4524fn faddeeva_upper_halfplane(z: Complex) -> Complex {
4537 let (l, coeffs) = faddeeva_weideman_coeffs();
4538 let iz = Complex {
4539 re: -z.im,
4540 im: z.re,
4541 }; let l_minus = Complex {
4543 re: l - iz.re,
4544 im: -iz.im,
4545 }; let l_plus = Complex {
4547 re: l + iz.re,
4548 im: iz.im,
4549 }; let zz = complex_div(l_plus, l_minus); let mut p = Complex {
4553 re: coeffs[0],
4554 im: 0.0,
4555 };
4556 for &c in &coeffs[1..] {
4557 p = complex_add(complexmul(p, zz), Complex { re: c, im: 0.0 });
4558 }
4559 let l_minus_sq = complexmul(l_minus, l_minus);
4560 let term1 = complex_div(
4561 Complex {
4562 re: 2.0 * p.re,
4563 im: 2.0 * p.im,
4564 },
4565 l_minus_sq,
4566 );
4567 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI;
4568 let term2 = complex_div(
4569 Complex {
4570 re: inv_sqrt_pi,
4571 im: 0.0,
4572 },
4573 l_minus,
4574 );
4575 complex_add(term1, term2)
4576}
4577
4578const FADDEEVA_WEIDEMAN_N: usize = 44;
4581
4582fn faddeeva_weideman_coeffs() -> &'static (f64, [f64; FADDEEVA_WEIDEMAN_N]) {
4588 static CACHE: OnceLock<(f64, [f64; FADDEEVA_WEIDEMAN_N])> = OnceLock::new();
4589 CACHE.get_or_init(|| {
4590 let n = FADDEEVA_WEIDEMAN_N;
4591 let l = (n as f64 / SQRT_2).sqrt();
4592 let m = 2 * n;
4593 let m2 = 2 * m; let mut f = vec![0.0_f64; m2];
4597 for (idx, fi) in f.iter_mut().enumerate().skip(1) {
4598 let k = (idx as isize - 1) - (m as isize - 1);
4599 let theta = (k as f64) * std::f64::consts::PI / (m as f64);
4600 let t = l * (0.5 * theta).tan();
4601 *fi = (-t * t).exp() * (l * l + t * t);
4602 }
4603 let half = m2 / 2;
4606 let mut coeffs = [0.0_f64; FADDEEVA_WEIDEMAN_N];
4607 for j in 1..=n {
4608 let mut acc = 0.0_f64;
4609 for (p, _) in f.iter().enumerate() {
4610 let fp = f[(p + half) % m2];
4611 if fp != 0.0 {
4612 acc += fp
4613 * (-2.0 * std::f64::consts::PI * (j as f64) * (p as f64) / (m2 as f64))
4614 .cos();
4615 }
4616 }
4617 coeffs[n - j] = acc / (m2 as f64);
4619 }
4620 (l, coeffs)
4621 })
4622}
4623
4624#[inline]
4626fn sigmoid(x: f64) -> f64 {
4627 let x_clamped = x.clamp(-QUADRATURE_EXP_LOG_MAX, QUADRATURE_EXP_LOG_MAX);
4628 1.0 / (1.0 + f64::exp(-x_clamped))
4629}
4630
4631#[derive(Clone, Copy, Debug)]
4647pub struct CLogLogConvolutionDerivatives {
4648 pub l: f64,
4650
4651 pub l_mu: f64,
4653 pub l_sigma: f64,
4654
4655 pub l_mumu: f64,
4657 pub l_musigma: f64,
4658 pub l_sigmasigma: f64,
4659
4660 pub l_mumumu: f64,
4662 pub l_mumusigma: f64,
4663 pub l_musigmasigma: f64,
4664 pub l_sigmasigmasigma: f64,
4665
4666 pub l_mumumumu: f64,
4668 pub l_mumumusigma: f64,
4669 pub l_mumusigmasigma: f64,
4670 pub l_musigmasigmasigma: f64,
4671 pub l_sigmasigmasigmasigma: f64,
4672}
4673
4674#[inline]
4675pub(crate) fn cloglog_point_jet5(t: f64) -> (f64, f64, f64, f64, f64, f64) {
4676 if t.is_nan() {
4677 return (f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN);
4678 }
4679 let et = safe_exp(t);
4680
4681 (
4682 -(-et).exp_m1(),
4683 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0]),
4684 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -1.0]),
4685 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -3.0, 1.0]),
4686 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -7.0, 6.0, -1.0]),
4687 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -15.0, 25.0, -10.0, 1.0]),
4688 )
4689}
4690
4691#[inline]
4703fn cloglog_g_derivatives(t: f64) -> (f64, f64, f64, f64, f64) {
4704 let (g, g1, g2, g3, g4, _) = cloglog_point_jet5(t);
4705 (g, g1, g2, g3, g4)
4706}
4707
4708pub fn cloglog_ghq_value(ctx: &QuadratureContext, mu: f64, sigma: f64, n_nodes: usize) -> f64 {
4716 if sigma.abs() < 1e-14 {
4717 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
4718 return g.clamp(0.0, 1.0);
4719 }
4720 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
4721
4722 let inv_sig2 = 1.0 / (sigma * sigma);
4754 let mut eta_hat = mu;
4755 let mut converged = false;
4756 for _ in 0..100 {
4757 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
4758 if !(g > 0.0) || !g1.is_finite() || !g2.is_finite() {
4759 break;
4760 }
4761 let r = g1 / g;
4762 let lp = r - (eta_hat - mu) * inv_sig2;
4763 let lpp = g2 / g - r * r - inv_sig2;
4764 if !lpp.is_finite() || lpp >= 0.0 {
4765 break;
4766 }
4767 let step = lp / lpp;
4768 eta_hat -= step;
4769 if step.abs() <= 1e-13 * (1.0 + eta_hat.abs()) {
4770 converged = true;
4771 break;
4772 }
4773 }
4774
4775 let tau = if converged {
4779 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
4780 if g > 0.0 {
4781 let r = g1 / g;
4782 let lpp = g2 / g - r * r - inv_sig2;
4783 let tau2 = -1.0 / lpp;
4784 if tau2.is_finite() && tau2 > 0.0 {
4785 Some(tau2.sqrt())
4786 } else {
4787 None
4788 }
4789 } else {
4790 None
4791 }
4792 } else {
4793 None
4794 };
4795
4796 let eval_at = |n: usize| -> f64 {
4798 match tau {
4799 Some(tau) => {
4800 let pref = tau * inv_sqrt_pi / sigma;
4801 with_gh_nodesweights(ctx, n, |nodes, weights| {
4802 let mut sum = 0.0_f64;
4803 for i in 0..nodes.len() {
4804 let t = nodes[i];
4805 let eta_i = eta_hat + SQRT_2 * tau * t;
4806 let (g, _, _, _, _, _) = cloglog_point_jet5(eta_i);
4807 let dev = eta_i - mu;
4808 sum += weights[i] * (t * t - 0.5 * dev * dev * inv_sig2).exp() * g;
4809 }
4810 (pref * sum).clamp(0.0, 1.0)
4811 })
4812 }
4813 None => {
4814 let scale = SQRT_2 * sigma;
4815 with_gh_nodesweights(ctx, n, |nodes, weights| {
4816 let mut sum = 0.0_f64;
4817 for i in 0..nodes.len() {
4818 let t = mu + scale * nodes[i];
4819 let (g, _, _, _, _) = cloglog_g_derivatives(t);
4820 sum += weights[i] * g;
4821 }
4822 (sum * inv_sqrt_pi).clamp(0.0, 1.0)
4823 })
4824 }
4825 }
4826 };
4827
4828 const CLOGLOG_GHQ_ORDER_LADDER: [usize; 5] = [7, 15, 21, 31, 51];
4833 const CLOGLOG_GHQ_CONV_TOL: f64 = 1e-10;
4834 let floor = n_nodes.min(*CLOGLOG_GHQ_ORDER_LADDER.last().unwrap());
4835 let mut prev: Option<f64> = None;
4836 let mut result = 0.0_f64;
4837 for &n in CLOGLOG_GHQ_ORDER_LADDER.iter().filter(|&&n| n >= floor) {
4838 let cur = eval_at(n);
4839 result = cur;
4840 if let Some(p) = prev
4841 && (cur - p).abs() < CLOGLOG_GHQ_CONV_TOL
4842 {
4843 break;
4844 }
4845 prev = Some(cur);
4846 }
4847 result
4848}
4849
4850pub fn cloglog_ghq_derivatives(
4861 ctx: &QuadratureContext,
4862 mu: f64,
4863 sigma: f64,
4864 n_nodes: usize,
4865) -> CLogLogConvolutionDerivatives {
4866 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
4867
4868 if sigma.abs() < 1e-14 {
4875 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
4876 return CLogLogConvolutionDerivatives {
4877 l: g,
4878 l_mu: g1,
4879 l_sigma: 0.0,
4880 l_mumu: g2,
4881 l_musigma: 0.0,
4882 l_sigmasigma: g2,
4883 l_mumumu: g3,
4884 l_mumusigma: 0.0,
4885 l_musigmasigma: g3,
4886 l_sigmasigmasigma: 0.0,
4887 l_mumumumu: g4,
4888 l_mumumusigma: 0.0,
4889 l_mumusigmasigma: g4,
4890 l_musigmasigmasigma: 0.0,
4891 l_sigmasigmasigmasigma: 3.0 * g4,
4892 };
4893 }
4894
4895 let scale = SQRT_2 * sigma;
4896 let sqrt2 = SQRT_2;
4897
4898 with_gh_nodesweights(ctx, n_nodes, |nodes, weights| {
4899 let mut s = [[0.0_f64; 5]; 5];
4911
4912 for i in 0..nodes.len() {
4913 let x = nodes[i];
4914 let t = mu + scale * x;
4915 let (g0, g1, g2, g3, g4) = cloglog_g_derivatives(t);
4916 let w = weights[i];
4917
4918 let x2 = x * x;
4920 let x3 = x2 * x;
4921 let x4 = x3 * x;
4922
4923 s[0][0] += w * g0;
4925
4926 s[1][0] += w * g1;
4928 s[1][1] += w * x * g1;
4929
4930 s[2][0] += w * g2;
4932 s[2][1] += w * x * g2;
4933 s[2][2] += w * x2 * g2;
4934
4935 s[3][0] += w * g3;
4937 s[3][1] += w * x * g3;
4938 s[3][2] += w * x2 * g3;
4939 s[3][3] += w * x3 * g3;
4940
4941 s[4][0] += w * g4;
4943 s[4][1] += w * x * g4;
4944 s[4][2] += w * x2 * g4;
4945 s[4][3] += w * x3 * g4;
4946 s[4][4] += w * x4 * g4;
4947 }
4948
4949 let sqrt2_1 = sqrt2;
4952 let sqrt2_2 = 2.0; let sqrt2_3 = 2.0 * sqrt2; let sqrt2_4 = 4.0; CLogLogConvolutionDerivatives {
4957 l: inv_sqrt_pi * s[0][0],
4959
4960 l_mu: inv_sqrt_pi * s[1][0],
4962 l_sigma: inv_sqrt_pi * sqrt2_1 * s[1][1],
4963
4964 l_mumu: inv_sqrt_pi * s[2][0],
4966 l_musigma: inv_sqrt_pi * sqrt2_1 * s[2][1],
4967 l_sigmasigma: inv_sqrt_pi * sqrt2_2 * s[2][2],
4968
4969 l_mumumu: inv_sqrt_pi * s[3][0],
4971 l_mumusigma: inv_sqrt_pi * sqrt2_1 * s[3][1],
4972 l_musigmasigma: inv_sqrt_pi * sqrt2_2 * s[3][2],
4973 l_sigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[3][3],
4974
4975 l_mumumumu: inv_sqrt_pi * s[4][0],
4977 l_mumumusigma: inv_sqrt_pi * sqrt2_1 * s[4][1],
4978 l_mumusigmasigma: inv_sqrt_pi * sqrt2_2 * s[4][2],
4979 l_musigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[4][3],
4980 l_sigmasigmasigmasigma: inv_sqrt_pi * sqrt2_4 * s[4][4],
4981 }
4982 })
4983}
4984
4985pub fn cloglog_ghq_derivatives_adaptive(
4991 ctx: &QuadratureContext,
4992 mu: f64,
4993 sigma: f64,
4994) -> CLogLogConvolutionDerivatives {
4995 let n = adaptive_point_count_from_sd(sigma.abs());
4996 cloglog_ghq_derivatives(ctx, mu, sigma, n)
4997}
4998
4999#[cfg(test)]
5000mod tests {
5001 use super::*;
5002 use approx::assert_relative_eq;
5003 use gam_problem::LikelihoodScaleMetadata;
5004 use gam_spec::LikelihoodSpec;
5005
5006 #[test]
5013 fn log_half_erfc_stable_matches_high_precision_reference() {
5014 let refs: &[(f64, f64)] = &[
5015 (-3.0, -1.1045309498499094e-5),
5016 (-1.5, -0.017092677825984745),
5017 (-0.5, -0.27410803278438573),
5018 (0.0, -0.69314718055994531),
5019 (0.7, -1.8257336940742865),
5020 (2.0, -6.0580884451765829),
5021 (5.0, -27.89403672609738),
5022 (12.0, -147.75386135854695),
5023 ];
5024 for &(u, reference) in refs {
5025 let got = log_half_erfc_stable(u);
5026 let rel = (got - reference).abs() / reference.abs().max(1.0e-6);
5027 assert!(
5028 rel < 1.0e-12,
5029 "log_half_erfc_stable({u}) = {got:.17e}, reference {reference:.17e}, \
5030 rel {rel:.3e} >= 1e-12"
5031 );
5032 }
5033 }
5034
5035 pub(crate) fn cloglog_posterior_meanwith_deriv_gamma_reference(
5036 mu: f64,
5037 sigma: f64,
5038 ) -> Result<IntegratedMeanDerivative, EstimationError> {
5039 let survival = cloglog_survival_gamma_reference(mu, sigma)?;
5042 let shifted_survival = cloglog_survival_gamma_reference(mu + sigma * sigma, sigma)?;
5043 let mean = cloglog_mean_from_survival(survival);
5044 let dmean = cloglog_shift_identity_derivative(mu, sigma, shifted_survival);
5045 if !(mean.is_finite() && dmean.is_finite()) {
5046 crate::bail_invalid_estim!(
5047 "Gamma cloglog reference backend produced non-finite values"
5048 );
5049 }
5050 Ok(IntegratedMeanDerivative {
5051 mean,
5052 dmean_dmu: dmean.max(0.0),
5053 mode: IntegratedExpectationMode::ExactSpecialFunction,
5054 })
5055 }
5056
5057 fn even_moment_exp_neg_x2(power: usize) -> f64 {
5058 assert!(power.is_multiple_of(2));
5059 let m = power / 2;
5060 let mut odd_double_factorial = 1.0_f64;
5061 for k in 0..m {
5062 odd_double_factorial *= (2 * k + 1) as f64;
5063 }
5064 odd_double_factorial * std::f64::consts::PI.sqrt() / 2.0_f64.powi(m as i32)
5065 }
5066
5067 fn normal_pdf(z: f64) -> f64 {
5068 (-(z * z) * 0.5).exp() / (2.0 * std::f64::consts::PI).sqrt()
5069 }
5070
5071 fn high_res_sigmoid_integral(eta: f64, se: f64) -> f64 {
5072 let a = -12.0_f64;
5074 let b = 12.0_f64;
5075 let n = 20_000usize; let h = (b - a) / n as f64;
5077
5078 let integrand = |z: f64| -> f64 { sigmoid(eta + se * z) * normal_pdf(z) };
5079
5080 let mut sum = integrand(a) + integrand(b);
5081 for i in 1..n {
5082 let x = a + (i as f64) * h;
5083 if i % 2 == 0 {
5084 sum += 2.0 * integrand(x);
5085 } else {
5086 sum += 4.0 * integrand(x);
5087 }
5088 }
5089 sum * h / 3.0
5090 }
5091
5092 #[test]
5093 fn test_computed_nodes_symmetric() {
5094 let ctx = QuadratureContext::new();
5096 let gh = ctx.gauss_hermite();
5097 for i in 0..N_POINTS / 2 {
5098 let j = N_POINTS - 1 - i;
5099 assert_relative_eq!(gh.nodes[i], -gh.nodes[j], epsilon = 1e-12);
5100 }
5101 assert_relative_eq!(gh.nodes[N_POINTS / 2], 0.0, epsilon = 1e-12);
5103 }
5104
5105 #[test]
5106 fn test_computedweights_symmetric() {
5107 let ctx = QuadratureContext::new();
5109 let gh = ctx.gauss_hermite();
5110 for i in 0..N_POINTS / 2 {
5111 let j = N_POINTS - 1 - i;
5112 assert_relative_eq!(gh.weights[i], gh.weights[j], epsilon = 1e-12);
5113 }
5114 }
5115
5116 #[test]
5117 fn testweights_sum_to_sqrt_pi() {
5118 let ctx = QuadratureContext::new();
5120 let gh = ctx.gauss_hermite();
5121 let sum: f64 = gh.weights.iter().sum();
5122 assert_relative_eq!(sum, std::f64::consts::PI.sqrt(), epsilon = 1e-10);
5123 }
5124
5125 #[test]
5126 fn test_clenshaw_curtisweights_are_symmetric_and_integrate_constants() {
5127 let rule = compute_clenshaw_curtis_n(33);
5128 let m = rule.weights.len() - 1;
5129 for j in 0..=m / 2 {
5130 assert_relative_eq!(rule.nodes[j], -rule.nodes[m - j], epsilon = 1e-14);
5131 assert_relative_eq!(rule.weights[j], rule.weights[m - j], epsilon = 1e-14);
5132 }
5133 let sum: f64 = rule.weights.iter().sum();
5134 assert_relative_eq!(sum, 2.0, epsilon = 1e-14, max_relative = 1e-14);
5135 }
5136
5137 #[test]
5138 fn test_cc_preference_prefers_moderate_central_case() {
5139 assert!(cloglog_should_prefer_cc(-0.2, 0.8, CLOGLOG_CC_TOL));
5140 }
5141
5142 #[test]
5143 fn test_cc_preference_prefers_moderately_large_case() {
5144 assert!(cloglog_should_prefer_cc(0.0, 2.0, CLOGLOG_CC_TOL));
5145 }
5146
5147 #[test]
5148 fn test_cc_preference_rejects_broad_case() {
5149 assert!(!cloglog_should_prefer_cc(0.0, 5.0, CLOGLOG_CC_TOL));
5150 }
5151
5152 #[test]
5153 fn testwilkinson_shift_finitewhen_d_iszero() {
5154 let shift = wilkinson_shift(0.0, 0.0, 1.25);
5157 assert!(shift.is_finite());
5158 assert_relative_eq!(shift, -1.25, epsilon = 1e-14);
5159 }
5160
5161 #[test]
5162 fn test_matches_abramowitz_stegun_7_point_gauss_hermite_constants() {
5163 let known_nodes = [
5167 -2.651_961_356_835_233_4,
5168 -1.673_551_628_767_471_4,
5169 -0.816_287_882_858_964_7,
5170 0.0,
5171 0.816_287_882_858_964_7,
5172 1.673_551_628_767_471_4,
5173 2.651_961_356_835_233_4,
5174 ];
5175 let knownweights = [
5176 0.000_971_781_245_099_519_1,
5177 0.054_515_582_819_127_03,
5178 0.425_607_252_610_127_8,
5179 0.810_264_617_556_807_3,
5180 0.425_607_252_610_127_8,
5181 0.054_515_582_819_127_03,
5182 0.000_971_781_245_099_519_1,
5183 ];
5184
5185 let ctx = QuadratureContext::new();
5186 let gh = ctx.gauss_hermite();
5187 for i in 0..N_POINTS {
5188 assert_relative_eq!(gh.nodes[i], known_nodes[i], epsilon = 1e-12);
5189 assert_relative_eq!(gh.weights[i], knownweights[i], epsilon = 1e-12);
5190 }
5191 }
5192
5193 #[test]
5194 fn test_gauss_hermite_weight_assembly_uses_eigenvector_rows() {
5195 let mut diag = [0.0_f64; N_POINTS];
5196 let mut off_diag = [0.0_f64; N_POINTS - 1];
5197 for (i, od) in off_diag.iter_mut().enumerate() {
5198 *od = (((i + 1) as f64) / 2.0).sqrt();
5199 }
5200 let (nodes, eigenvectors) = symmetric_tridiagonal_eigen(&mut diag, &mut off_diag);
5201 let mu0 = std::f64::consts::PI.sqrt();
5202 let mut row_pairs: Vec<(f64, f64)> = (0..N_POINTS)
5203 .map(|i| (nodes[i], mu0 * eigenvectors[i][0] * eigenvectors[i][0]))
5204 .collect();
5205 let mut column_pairs: Vec<(f64, f64)> = (0..N_POINTS)
5206 .map(|i| (nodes[i], mu0 * eigenvectors[0][i] * eigenvectors[0][i]))
5207 .collect();
5208 row_pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
5209 column_pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
5210
5211 let knownweights = [
5212 0.000_971_781_245_099_519_1,
5213 0.054_515_582_819_127_03,
5214 0.425_607_252_610_127_8,
5215 0.810_264_617_556_807_3,
5216 0.425_607_252_610_127_8,
5217 0.054_515_582_819_127_03,
5218 0.000_971_781_245_099_519_1,
5219 ];
5220
5221 for i in 0..N_POINTS {
5222 assert_relative_eq!(row_pairs[i].1, knownweights[i], epsilon = 1e-12);
5223 }
5224 let column_error: f64 = column_pairs
5225 .iter()
5226 .zip(knownweights.iter())
5227 .map(|(actual, expected)| (actual.1 - expected).abs())
5228 .sum();
5229 assert!(
5230 column_error > 1.0,
5231 "column-oriented eigenvector indexing unexpectedly matched A&S weights"
5232 );
5233 }
5234
5235 #[test]
5236 fn testzero_se_returns_mode() {
5237 let eta = 1.5;
5239 let se = 0.0;
5240 let ctx = QuadratureContext::new();
5241 let mean = logit_posterior_mean(&ctx, eta, se);
5242 let mode = sigmoid(eta);
5243 assert_relative_eq!(mean, mode, epsilon = 1e-10);
5244 }
5245
5246 #[test]
5247 fn test_symmetric_atzero() {
5248 let eta = 0.0;
5250 let se = 1.0;
5251 let ctx = QuadratureContext::new();
5252 let mean = logit_posterior_mean(&ctx, eta, se);
5253 assert_relative_eq!(mean, 0.5, epsilon = 0.01);
5255 }
5256
5257 #[test]
5258 fn test_shrinkage_at_extremes() {
5259 let eta = 3.0; let se = 1.0;
5262 let ctx = QuadratureContext::new();
5263 let mean = logit_posterior_mean(&ctx, eta, se);
5264 let mode = sigmoid(eta);
5265
5266 assert!(mean < mode, "Expected mean {} < mode {}", mean, mode);
5268 assert!(mean > 0.8, "Mean {} should still be high", mean);
5270 }
5271
5272 #[test]
5273 fn test_matches_monte_carlo() {
5274 let eta = 2.0;
5276 let se = 0.8;
5277
5278 let ctx = QuadratureContext::new();
5279 let quad_mean = logit_posterior_mean(&ctx, eta, se);
5280
5281 let n_samples = 100_000;
5283 let mut mc_sum = 0.0;
5284 let mut rng_state = 12345u64; for _ in 0..n_samples {
5286 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5288 let u1 = ((rng_state as f64) / (u64::MAX as f64)).max(1e-10); rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5290 let u2 = (rng_state as f64) / (u64::MAX as f64);
5291 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
5292 let eta_sample = eta + se * z;
5293 mc_sum += sigmoid(eta_sample);
5294 }
5295 let mc_mean = mc_sum / (n_samples as f64);
5296
5297 assert_relative_eq!(quad_mean, mc_mean, epsilon = 0.01);
5299 }
5300
5301 #[test]
5302 fn test_quadrature_integrates_x_squared() {
5303 let ctx = QuadratureContext::new();
5306 let gh = ctx.gauss_hermite();
5307 let mut sum = 0.0;
5308 for i in 0..N_POINTS {
5309 sum += gh.weights[i] * gh.nodes[i] * gh.nodes[i];
5310 }
5311 let expected = std::f64::consts::PI.sqrt() / 2.0;
5312 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5313 }
5314
5315 #[test]
5316 fn test_quadrature_integrates_x_fourth() {
5317 let ctx = QuadratureContext::new();
5320 let gh = ctx.gauss_hermite();
5321 let mut sum = 0.0;
5322 for i in 0..N_POINTS {
5323 let x = gh.nodes[i];
5324 sum += gh.weights[i] * x * x * x * x;
5325 }
5326 let expected = 3.0 * std::f64::consts::PI.sqrt() / 4.0;
5327 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5328 }
5329
5330 #[test]
5331 fn test_moment_exactness_up_to_degree_13() {
5332 let ctx = QuadratureContext::new();
5333 let gh = ctx.gauss_hermite();
5334
5335 for degree in 0..=13usize {
5336 let approx: f64 = (0..N_POINTS)
5337 .map(|i| gh.weights[i] * gh.nodes[i].powi(degree as i32))
5338 .sum();
5339
5340 let expected = if degree % 2 == 1 {
5341 0.0
5342 } else {
5343 even_moment_exp_neg_x2(degree)
5344 };
5345
5346 let err = (approx - expected).abs();
5347 let rel_scale = approx.abs().max(expected.abs()).max(1.0);
5348 assert!(
5349 err <= 1e-10 || err / rel_scale <= 1e-10,
5350 "degree={} approx={} expected={} abs_err={}",
5351 degree,
5352 approx,
5353 expected,
5354 err
5355 );
5356 }
5357 }
5358
5359 #[test]
5360 fn test_integrated_sigmoid_matches_high_res_integral_random_pairs() {
5361 let ctx = QuadratureContext::new();
5362 let mut rng_state = 0x4d595df4d0f33173u64;
5363
5364 for _ in 0..20 {
5365 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5366 let u_eta = (rng_state as f64) / (u64::MAX as f64);
5367 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5368 let u_se = (rng_state as f64) / (u64::MAX as f64);
5369
5370 let eta = -6.0 + 12.0 * u_eta;
5371 let se = 0.02 + 1.5 * u_se;
5372
5373 let ghq = logit_posterior_mean(&ctx, eta, se);
5374 let numeric = high_res_sigmoid_integral(eta, se);
5375 assert_relative_eq!(ghq, numeric, epsilon = 2e-3);
5376 }
5377 }
5378
5379 #[test]
5380 fn test_logit_posterior_derivative_remains_positive_in_positive_tail() {
5381 let eta = 20.0;
5382 let se = 0.0;
5383 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5384 .expect("logit posterior mean derivative should evaluate");
5385 assert!(dmu > 0.0);
5386 assert!(
5387 dmu < 1e-6,
5388 "positive-tail derivative should stay tiny but nonzero, got {dmu}"
5389 );
5390 }
5391
5392 #[test]
5393 fn test_logit_posterior_derivative_matches_central_difference() {
5394 let ctx = QuadratureContext::new();
5395 let eta = 1.7;
5396 let se = 0.9;
5397 let h = 1e-5;
5398
5399 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5400 .expect("logit posterior mean derivative should evaluate");
5401 let mu_plus = logit_posterior_mean(&ctx, eta + h, se);
5402 let mu_minus = logit_posterior_mean(&ctx, eta - h, se);
5403 let dmufd = (mu_plus - mu_minus) / (2.0 * h);
5404
5405 assert_eq!(dmu.signum(), dmufd.signum());
5406 assert_relative_eq!(dmu, dmufd, epsilon = 5e-6, max_relative = 2e-4);
5407 }
5408
5409 fn dense_sigmoid_normal_mean(mu: f64, sigma: f64) -> f64 {
5415 let a = -18.0_f64;
5416 let b = 18.0_f64;
5417 let n = 400_000usize; let h = (b - a) / n as f64;
5419 let integrand = |z: f64| -> f64 { sigmoid(mu + sigma * z) * normal_pdf(z) };
5420 let mut sum = integrand(a) + integrand(b);
5421 for i in 1..n {
5422 let z = a + (i as f64) * h;
5423 sum += if i % 2 == 0 { 2.0 } else { 4.0 } * integrand(z);
5424 }
5425 sum * h / 3.0
5426 }
5427
5428 #[test]
5429 fn test_logit_posterior_mean_exact_symmetry_identity() {
5430 let cases = [
5433 (-3.0, 0.5),
5434 (-1.2, 1.7),
5435 (0.0, 2.2),
5436 (2.3, 0.8),
5437 (3.0, 0.05),
5438 ];
5439 for (mu, sigma) in cases {
5440 let p = logit_posterior_mean_exact(mu, sigma);
5441 let q = logit_posterior_mean_exact(-mu, sigma);
5442 assert!(
5443 (p + q - 1.0).abs() < 1e-12,
5444 "symmetry broken at mu={mu} sigma={sigma}: p+q-1 = {:.3e}",
5445 p + q - 1.0
5446 );
5447 }
5448 }
5449
5450 #[test]
5451 fn test_logit_posterior_mean_exact_matches_high_res_integral() {
5452 let cases = [
5456 (-2.0, 0.4),
5457 (-0.7, 1.1),
5458 (0.8, 0.9),
5459 (2.4, 1.7),
5460 (3.0, 0.05),
5461 (3.0, 0.5),
5462 (-2.0, 2.0),
5463 (5.0, 3.0),
5464 ];
5465 for (mu, sigma) in cases {
5466 let exact = logit_posterior_mean_exact(mu, sigma);
5467 let numeric = dense_sigmoid_normal_mean(mu, sigma);
5468 assert!(
5469 (exact - numeric).abs() < 1e-10,
5470 "oracle ≠ dense reference at mu={mu} sigma={sigma}: \
5471 exact={exact:.13} ref={numeric:.13} err={:.3e}",
5472 (exact - numeric).abs()
5473 );
5474 }
5475 }
5476
5477 #[test]
5484 fn test_logit_posterior_mean_exact_no_truncation_bias_1459() {
5485 let table = [
5488 (1.0, 0.02),
5489 (1.0, 0.05),
5490 (1.0, 0.5),
5491 (1.0, 2.0),
5492 (3.0, 0.02),
5493 (3.0, 0.05),
5494 (3.0, 0.5),
5495 (3.0, 2.0),
5496 (-2.0, 0.02),
5497 (-2.0, 0.05),
5498 (-2.0, 0.5),
5499 (-2.0, 2.0),
5500 ];
5501 for (mu, sigma) in table {
5502 let exact = logit_posterior_mean_exact(mu, sigma);
5503 let reference = dense_sigmoid_normal_mean(mu, sigma);
5504 let err = (exact - reference).abs();
5505 assert!(
5506 err < 1e-10,
5507 "#1459 truncation bias resurfaced at mu={mu} sigma={sigma}: \
5508 err={err:.3e} (pre-fix bias here was ~{:.2e})",
5509 mu.abs() / (2.0 * std::f64::consts::PI.powi(2) * 4096.0)
5510 );
5511 }
5512
5513 let mu = 3.0;
5518 let errs: Vec<f64> = [0.05, 0.5, 2.0]
5519 .iter()
5520 .map(|&s| logit_posterior_mean_exact(mu, s) - dense_sigmoid_normal_mean(mu, s))
5521 .collect();
5522 for e in &errs {
5523 assert!(
5524 e.abs() < 1e-10,
5525 "residual {e:.3e} at mu=3 — old σ-independent plateau was 3.71e-5"
5526 );
5527 }
5528 }
5529
5530 #[test]
5538 fn test_faddeeva_weideman_matches_known_values() {
5539 let w0 = faddeeva_upper_halfplane(Complex { re: 0.0, im: 0.0 });
5541 assert!(
5542 (w0.re - 1.0).abs() < 1e-13 && w0.im.abs() < 1e-13,
5543 "w(0)={w0:?}"
5544 );
5545 let on_axis = [
5547 (0.1, 0.8964569799691268),
5548 (0.5, 0.6156903441929258),
5549 (1.0, 0.427583576155807),
5550 (2.0, 0.2553956763105058),
5551 (5.0, 0.11070463773306861),
5552 (9.0, 0.06230772403777468),
5553 ];
5554 for (y, want) in on_axis {
5555 let w = faddeeva_upper_halfplane(Complex { re: 0.0, im: y });
5556 assert!(
5557 (w.re - want).abs() < 1e-13 && w.im.abs() < 1e-13,
5558 "w(i·{y}): got {w:?}, want re={want}, err={:.2e}",
5559 (w.re - want).abs()
5560 );
5561 }
5562 let off_axis = [
5564 ((0.7, 1.3), (0.31327301971562715, 0.12443489420104513)),
5565 ((-1.5, 0.8), (0.21066359024766423, -0.27001624496296617)),
5566 ((3.0, 0.4), (0.030278754646989155, 0.1957320888774461)),
5567 ];
5568 for ((re, im), (wre, wim)) in off_axis {
5569 let w = faddeeva_upper_halfplane(Complex { re, im });
5570 assert!(
5571 (w.re - wre).abs() < 1e-13 && (w.im - wim).abs() < 1e-13,
5572 "w({re}+{im}i): got {w:?}, want ({wre},{wim})"
5573 );
5574 }
5575 let w = faddeeva_upper_halfplane(Complex { re: 3.0, im: 40.0 });
5579 assert!(
5580 (w.re - 0.01402158696172506).abs() < 1e-13
5581 && (w.im - 0.0010509664408184546).abs() < 1e-13,
5582 "tail value mismatch: w={w:?}"
5583 );
5584 }
5585
5586 #[test]
5587 fn test_integrated_logit_mean_close_to_exact_oracle() {
5588 let ctx = QuadratureContext::new();
5592 let cases = [(-3.0, 0.3), (-1.0, 0.8), (0.5, 1.2), (2.8, 1.0)];
5593 for (eta, se) in cases {
5594 let ghq = logit_posterior_mean(&ctx, eta, se);
5595 let exact = logit_posterior_mean_exact(eta, se);
5596 assert!(
5597 (ghq - exact).abs() < 1e-6,
5598 "production path drifts from oracle at eta={eta} se={se}: \
5599 ghq={ghq:.12} oracle={exact:.12} gap={:.3e}",
5600 (ghq - exact).abs()
5601 );
5602 }
5603 }
5604
5605 #[test]
5606 fn test_probit_posterior_mean_reduces_to_map_atzero_se() {
5607 let eta = 1.25;
5608 let p = probit_posterior_mean(eta, 0.0);
5609 let map = gam_math::probability::normal_cdf(eta);
5610 assert_relative_eq!(p, map, epsilon = 1e-12);
5611 }
5612
5613 #[test]
5614 fn test_probit_posterior_mean_shrinks_extremeswith_uncertainty() {
5615 let hi_eta = 3.0;
5616 let lo_eta = -3.0;
5617 let p_hi_map = probit_posterior_mean(hi_eta, 0.0);
5618 let p_hi_unc = probit_posterior_mean(hi_eta, 2.0);
5619 let p_lo_map = probit_posterior_mean(lo_eta, 0.0);
5620 let p_lo_unc = probit_posterior_mean(lo_eta, 2.0);
5621 assert!(p_hi_unc < p_hi_map);
5622 assert!(p_lo_unc > p_lo_map);
5623 }
5624
5625 #[test]
5626 fn test_survival_posterior_mean_is_bounded_and_shrinks_tail() {
5627 let ctx = QuadratureContext::new();
5628 let eta: f64 = 3.0;
5629 let map = (-(eta.exp())).exp();
5630 let pm = survival_posterior_mean(&ctx, eta, 1.5);
5631 assert!((0.0..=1.0).contains(&pm));
5632 assert!(pm > map);
5633 }
5634
5635 #[test]
5636 fn test_cloglog_and_survival_posterior_means_are_complements() {
5637 let ctx = QuadratureContext::new();
5638 let cases = [
5639 (-3.0, 0.0),
5640 (-0.2, 0.1),
5641 (0.4, 0.8),
5642 (2.0, 1.5),
5643 (10.0, 0.3),
5644 (0.0, 20.0),
5645 (10.0, 10.0),
5646 (-0.5, 100.0),
5647 ];
5648 for (eta, se) in cases {
5649 let clog = cloglog_posterior_mean(&ctx, eta, se);
5650 let surv = survival_posterior_mean(&ctx, eta, se);
5651 assert_relative_eq!(clog + surv, 1.0, epsilon = 2e-10, max_relative = 2e-10);
5652 }
5653 }
5654
5655 #[test]
5656 fn test_cloglog_and_survival_share_large_sigmaspecial_function_path() {
5657 let ctx = QuadratureContext::new();
5658 let eta = -0.2;
5659 let se = 0.8;
5660 let clog = cloglog_posterior_mean(&ctx, eta, se);
5661 let surv = survival_posterior_mean(&ctx, eta, se);
5662 let integrated =
5663 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, eta, se)
5664 .expect("cloglog integrated inverse-link moments should evaluate");
5665 assert_eq!(
5666 integrated.mode,
5667 IntegratedExpectationMode::ExactSpecialFunction
5668 );
5669 assert_relative_eq!(clog, integrated.mean, epsilon = 1e-12, max_relative = 1e-12);
5670 assert_relative_eq!(clog + surv, 1.0, epsilon = 1e-10, max_relative = 1e-10);
5671 }
5672
5673 #[test]
5674 fn test_cloglog_and_survival_posteriorvariances_match() {
5675 let ctx = QuadratureContext::new();
5676 let cases = [(-3.0, 0.0), (-0.2, 0.1), (0.4, 0.8), (2.0, 1.5)];
5677 for (eta, se) in cases {
5678 let (_, clogvar) = cloglog_posterior_meanvariance(&ctx, eta, se);
5679 let (_, survvar) = survival_posterior_meanvariance(&ctx, eta, se);
5680 assert_relative_eq!(clogvar, survvar, epsilon = 1e-12, max_relative = 1e-12);
5681 }
5682 }
5683
5684 #[test]
5685 fn test_survivalvariance_uses_exactsecond_moment_shift() {
5686 let ctx = QuadratureContext::new();
5687 let eta = -0.2;
5688 let se = 0.8;
5689 let (survival, _) = cloglog_survival_term_controlled(&ctx, eta, se);
5690 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(&ctx, eta, se);
5691 let (_, variance) = survival_posterior_meanvariance(&ctx, eta, se);
5692 assert_relative_eq!(
5693 variance,
5694 (survival_sq - survival * survival).max(0.0),
5695 epsilon = 1e-12,
5696 max_relative = 1e-12
5697 );
5698 }
5699
5700 #[test]
5701 fn test_lognormal_laplace_shift_matches_explicitmu_plus_logz() {
5702 let ctx = QuadratureContext::new();
5703 let mu = -0.2;
5704 let sigma = 0.8;
5705 let z = 2.0;
5706 let shifted = lognormal_laplace_term_controlled(&ctx, z, mu, sigma);
5707 let explicit = cloglog_survival_term_controlled(&ctx, mu + z.ln(), sigma);
5708 assert_eq!(shifted.1, explicit.1);
5709 assert_relative_eq!(shifted.0, explicit.0, epsilon = 1e-12, max_relative = 1e-12);
5710 }
5711
5712 #[test]
5713 fn test_integrated_dispatch_uses_closed_form_probit() {
5714 let ctx = QuadratureContext::new();
5715 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Probit, 0.7, 1.3)
5716 .expect("probit integrated inverse-link moments should evaluate");
5717 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
5718 let direct = probit_posterior_meanwith_deriv_exact(0.7, 1.3);
5719 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
5720 assert_relative_eq!(out.dmean_dmu, direct.dmean_dmu, epsilon = 1e-12);
5721 }
5722
5723 #[test]
5724 fn test_integrated_probit_jet_matches_closed_form_derivatives() {
5725 let ctx = QuadratureContext::new();
5726 let mu = 0.7;
5727 let sigma = 1.3;
5728 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Probit, mu, sigma)
5729 .expect("probit integrated inverse-link jet should evaluate");
5730 let s = (1.0 + sigma * sigma).sqrt();
5731 let z = mu / s;
5732 let pdf = gam_math::probability::normal_pdf(z);
5733 assert_relative_eq!(
5734 out.mean,
5735 gam_math::probability::normal_cdf(z),
5736 epsilon = 1e-12
5737 );
5738 assert_relative_eq!(out.d1, pdf / s, epsilon = 1e-12);
5739 assert_relative_eq!(out.d2, -z * pdf / (s * s), epsilon = 1e-12);
5740 assert_relative_eq!(out.d3, (z * z - 1.0) * pdf / (s * s * s), epsilon = 1e-12);
5741 }
5742
5743 #[test]
5744 fn test_integrated_logit_jet_matches_central_differences() {
5745 let ctx = QuadratureContext::new();
5758 let mu = 1.1;
5759 let sigma = 0.8;
5760 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
5761 .expect("logit integrated inverse-link jet should evaluate");
5762 assert!(matches!(
5763 out.mode,
5764 IntegratedExpectationMode::ExactSpecialFunction
5765 | IntegratedExpectationMode::QuadratureFallback
5766 ));
5767 let (ref_mean, ref_d1, ref_d2, ref_d3) = logit_reference_jet_highres_simpson(mu, sigma);
5768 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
5769 assert_relative_eq!(out.d1, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
5770 assert_relative_eq!(out.d2, ref_d2, epsilon = 1e-11, max_relative = 1e-10);
5771 assert_relative_eq!(out.d3, ref_d3, epsilon = 1e-11, max_relative = 1e-10);
5772 }
5773
5774 #[test]
5775 fn test_integrated_logit_pirls_jet_matches_general_dispatch() {
5776 let ctx = QuadratureContext::new();
5786 let mu = 1.1;
5787 let sigma = 0.8;
5788
5789 let pirls =
5790 integrated_logit_inverse_link_jet_pirls(&ctx, mu, sigma).expect("PIRLS logit jet");
5791 let general = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
5792 .expect("general logit jet");
5793
5794 assert!(matches!(
5795 pirls.mode,
5796 IntegratedExpectationMode::ExactSpecialFunction
5797 | IntegratedExpectationMode::QuadratureFallback
5798 ));
5799 assert_eq!(pirls.mode, general.mode);
5800 assert_relative_eq!(pirls.mean, general.mean, epsilon = 1e-12);
5801 assert_relative_eq!(pirls.d1, general.d1, epsilon = 1e-12);
5802 assert_relative_eq!(pirls.d2, general.d2, epsilon = 1e-10);
5803 assert_relative_eq!(pirls.d3, general.d3, epsilon = 1e-8);
5804 }
5805
5806 #[test]
5807 fn test_integrated_cloglog_jet_matches_central_differences() {
5808 let ctx = QuadratureContext::new();
5809 let mu = 0.4;
5810 let sigma = 0.6;
5811 let h = 1e-4;
5812 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
5813 .expect("cloglog integrated inverse-link jet should evaluate");
5814 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
5815 .expect("cloglog integrated inverse-link jet should evaluate");
5816 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
5817 .expect("cloglog integrated inverse-link jet should evaluate");
5818 let d1fd = (plus.mean - minus.mean) / (2.0 * h);
5819 let d2fd = (plus.d1 - minus.d1) / (2.0 * h);
5820 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
5821 assert_eq!(out.d1.signum(), d1fd.signum());
5822 assert_eq!(out.d2.signum(), d2fd.signum());
5823 assert_eq!(out.d3.signum(), d3fd.signum());
5824 assert_relative_eq!(out.d1, d1fd, epsilon = 2e-5, max_relative = 3e-4);
5825 assert_relative_eq!(out.d2, d2fd, epsilon = 4e-5, max_relative = 8e-4);
5826 assert_relative_eq!(out.d3, d3fd, epsilon = 8e-5, max_relative = 2e-3);
5827 }
5828
5829 #[test]
5830 fn test_integrated_cloglog_wide_sigma_d3_matches_simpson_and_d2_slope() {
5831 let ctx = QuadratureContext::new();
5832 let cases = [(0.0, 4.0), (-1.0, 4.0), (2.0, 3.0), (3.0, 3.0)];
5833 let h = 1e-4;
5834
5835 for (mu, sigma) in cases {
5836 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
5837 .expect("wide-sigma cloglog integrated jet should evaluate");
5838 let reference = cloglog_reference_jet_highres_simpson(mu, sigma);
5839 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
5840 .expect("wide-sigma cloglog integrated jet should evaluate");
5841 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
5842 .expect("wide-sigma cloglog integrated jet should evaluate");
5843 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
5844
5845 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
5846 assert_relative_eq!(out.mean, reference.0, epsilon = 4e-8, max_relative = 4e-8);
5847 assert_relative_eq!(out.d1, reference.1, epsilon = 4e-8, max_relative = 4e-8);
5848 assert_relative_eq!(out.d2, reference.2, epsilon = 2e-9, max_relative = 2e-7);
5849 assert_relative_eq!(out.d3, reference.3, epsilon = 2e-9, max_relative = 2e-7);
5850 assert_relative_eq!(out.d3, d3fd, epsilon = 2e-7, max_relative = 4e-5);
5851 }
5852 }
5853
5854 #[test]
5855 fn test_latent_cloglog_jet5_matches_higher_order_central_differences() {
5856 let ctx = QuadratureContext::new();
5857 let mu = 0.35;
5858 let sigma = 0.7;
5859 let h = 2e-4;
5860
5861 let out = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu, sigma);
5862 let plus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu + h, sigma);
5863 let minus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu - h, sigma);
5864
5865 let d4fd = (plus.d3 - minus.d3) / (2.0 * h);
5866 let d5fd = (plus.d4 - minus.d4) / (2.0 * h);
5867
5868 assert_eq!(out.d4.signum(), d4fd.signum());
5869 assert_eq!(out.d5.signum(), d5fd.signum());
5870 assert_relative_eq!(out.d4, d4fd, epsilon = 2e-4, max_relative = 5e-3);
5871 assert_relative_eq!(out.d5, d5fd, epsilon = 6e-4, max_relative = 2e-2);
5872 }
5873
5874 #[test]
5875 fn test_logit_exact_derivative_matches_finite_difference() {
5876 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("controlled logit");
5886 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
5887 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
5888 assert!(out.dmean_dmu > 0.0);
5889 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
5890 }
5891
5892 #[test]
5893 fn test_logit_exact_clamped_degenerate_branch_is_locally_flat() {
5894 let out = logit_posterior_meanwith_deriv_exact(-710.0, 0.0).expect("exact logit");
5895 let h = 1e-6;
5896 let plus = logit_posterior_meanwith_deriv_exact(-710.0 + h, 0.0)
5897 .expect("exact logit plus")
5898 .mean;
5899 let minus = logit_posterior_meanwith_deriv_exact(-710.0 - h, 0.0)
5900 .expect("exact logit minus")
5901 .mean;
5902 let fd = (plus - minus) / (2.0 * h);
5903 assert_eq!(fd, 0.0);
5904 assert_eq!(out.dmean_dmu, 0.0);
5905 }
5906
5907 fn simpson_integrate<F>(a: f64, b: f64, n_intervals: usize, f: F) -> f64
5908 where
5909 F: Fn(f64) -> f64,
5910 {
5911 assert_eq!(n_intervals % 2, 0, "Simpson integration requires an even n");
5912 let h = (b - a) / n_intervals as f64;
5913 let mut sum = f(a) + f(b);
5914 for i in 1..n_intervals {
5915 let x = a + i as f64 * h;
5916 let w = if i % 2 == 0 { 2.0 } else { 4.0 };
5917 sum += w * f(x);
5918 }
5919 sum * h / 3.0
5920 }
5921
5922 fn cloglog_reference_mean_and_derivative(mu: f64, sigma: f64) -> (f64, f64) {
5923 if sigma <= CLOGLOG_SIGMA_DEGENERATE {
5924 return (cloglog_mean_exact(mu), cloglog_mean_d1_exact(mu));
5925 }
5926
5927 let z_max = 12.0;
5931 let n_intervals = 4096;
5932 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5933 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5934 let eta = mu + sigma * z;
5935 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_exact(eta)
5936 });
5937 let deriv = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5938 let eta = mu + sigma * z;
5939 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_d1_exact(eta)
5940 });
5941 (mean, deriv)
5942 }
5943
5944 fn logit_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
5957 let z_max = 14.0;
5958 let n_intervals = 16384;
5959 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5960 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
5961 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5962 let eta = mu + sigma * z;
5963 let (p, _, _, _) = component_point_jet(LinkComponent::Logit, eta);
5964 phi(z) * p
5965 });
5966 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5967 let eta = mu + sigma * z;
5968 let (_, p1, _, _) = component_point_jet(LinkComponent::Logit, eta);
5969 phi(z) * p1
5970 });
5971 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5972 let eta = mu + sigma * z;
5973 let (_, _, p2, _) = component_point_jet(LinkComponent::Logit, eta);
5974 phi(z) * p2
5975 });
5976 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5977 let eta = mu + sigma * z;
5978 let (_, _, _, p3) = component_point_jet(LinkComponent::Logit, eta);
5979 phi(z) * p3
5980 });
5981 (mean, d1, d2, d3)
5982 }
5983
5984 fn cloglog_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
5985 let z_max = 14.0;
5986 let n_intervals = 16384;
5987 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
5988 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
5989 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5990 let eta = mu + sigma * z;
5991 let (g, _, _, _, _, _) = cloglog_point_jet5(eta);
5992 phi(z) * g
5993 });
5994 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
5995 let eta = mu + sigma * z;
5996 let (_, g1, _, _, _, _) = cloglog_point_jet5(eta);
5997 phi(z) * g1
5998 });
5999 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6000 let eta = mu + sigma * z;
6001 let (_, _, g2, _, _, _) = cloglog_point_jet5(eta);
6002 phi(z) * g2
6003 });
6004 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6005 let eta = mu + sigma * z;
6006 let (_, _, _, g3, _, _) = cloglog_point_jet5(eta);
6007 phi(z) * g3
6008 });
6009 (mean, d1, d2, d3)
6010 }
6011
6012 #[test]
6013 fn test_cloglog_taylor_negative_tail_matches_mathematical_target() {
6014 let mu = -40.0;
6015 let sigma = 0.1;
6016 let out = cloglog_small_sigma_taylor(mu, sigma);
6017 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
6018
6019 assert!(
6020 out.dmean_dmu > 0.0,
6021 "negative-tail derivative should remain positive"
6022 );
6023 assert_relative_eq!(
6024 out.mean,
6025 expected_mean,
6026 epsilon = 1e-30,
6027 max_relative = 1e-12
6028 );
6029 assert_relative_eq!(
6030 out.dmean_dmu,
6031 expected_deriv,
6032 epsilon = 1e-30,
6033 max_relative = 1e-12
6034 );
6035 }
6036
6037 #[test]
6038 fn test_cloglog_degenerate_negative_tail_matches_pointwise_target() {
6039 let ctx = QuadratureContext::new();
6040 let mu = -40.0;
6041 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, 0.0);
6042
6043 assert!(
6044 out.dmean_dmu > 0.0,
6045 "degenerate negative-tail derivative should remain positive"
6046 );
6047 assert_relative_eq!(
6048 out.mean,
6049 cloglog_mean_exact(mu),
6050 epsilon = 1e-30,
6051 max_relative = 1e-15
6052 );
6053 assert_relative_eq!(
6054 out.dmean_dmu,
6055 cloglog_mean_d1_exact(mu),
6056 epsilon = 1e-30,
6057 max_relative = 1e-15
6058 );
6059 }
6060
6061 #[test]
6062 fn test_degenerate_probit_jet_is_exact_beyond_former_clamp() {
6063 let mu = -30.1;
6064 let probit = integrated_probit_jet(mu, 0.0);
6065 let pdf = gam_math::probability::normal_pdf(mu);
6066 assert!(
6067 pdf > 0.0,
6068 "test point must have a represented Gaussian tail"
6069 );
6070 assert_eq!(probit.mean, gam_math::probability::normal_cdf(mu));
6071 assert_eq!(probit.d1, pdf);
6072 assert_eq!(probit.d2, -mu * pdf);
6073 assert_eq!(probit.d3, (mu * mu - 1.0) * pdf);
6074
6075 let tail = (-710.0_f64).exp();
6091 assert!(
6092 tail > 0.0 && tail < f64::MIN_POSITIVE,
6093 "eta=-710 must sit in the subnormal tail, not underflow"
6094 );
6095 let logit = component_point_jet(LinkComponent::Logit, -710.0);
6096 assert_eq!(logit.0, tail);
6097 assert_eq!(logit.1, tail);
6098 assert_eq!(logit.2, tail);
6099 assert_eq!(logit.3, tail);
6100
6101 assert_eq!(
6103 (-750.0_f64).exp(),
6104 0.0,
6105 "eta=-750 must underflow f64 for this arm to mean anything"
6106 );
6107 let underflowed = component_point_jet(LinkComponent::Logit, -750.0);
6108 assert_eq!(underflowed.1, 0.0);
6109 assert_eq!(underflowed.2, 0.0);
6110 assert_eq!(underflowed.3, 0.0);
6111 }
6112
6113 #[test]
6114 fn test_degenerate_cloglog_component_jet_preserves_smooth_negative_tail() {
6115 let eta: f64 = -40.0;
6116 let t = eta.exp();
6117 let s = (-t).exp();
6118 let cloglog = component_point_jet(LinkComponent::CLogLog, eta);
6119 let expected_mean = -(-t).exp_m1();
6120 let expected_d1 = t * s;
6121 let expected_d2 = (t - t * t) * s;
6122 let expected_d3 = (t - 3.0 * t * t + t * t * t) * s;
6123
6124 assert!(cloglog.1 > 0.0, "negative-tail d1 should remain positive");
6125 assert_relative_eq!(
6126 cloglog.0,
6127 expected_mean,
6128 epsilon = 1e-30,
6129 max_relative = 1e-15
6130 );
6131 assert_relative_eq!(
6132 cloglog.1,
6133 expected_d1,
6134 epsilon = 1e-30,
6135 max_relative = 1e-15
6136 );
6137 assert_relative_eq!(
6138 cloglog.2,
6139 expected_d2,
6140 epsilon = 1e-30,
6141 max_relative = 1e-15
6142 );
6143 assert_relative_eq!(
6144 cloglog.3,
6145 expected_d3,
6146 epsilon = 1e-30,
6147 max_relative = 1e-15
6148 );
6149 }
6150
6151 #[test]
6152 fn test_zero_sigma_logit_and_cloglog_share_component_tail_jets() {
6153 let ctx = QuadratureContext::new();
6154 for (link, component, eta) in [
6155 (LinkFunction::Logit, LinkComponent::Logit, 50.0),
6156 (LinkFunction::CLogLog, LinkComponent::CLogLog, -50.0),
6157 ] {
6158 let integrated = integrated_inverse_link_jet(&ctx, link, eta, 0.0)
6159 .expect("degenerate integrated jet");
6160 let point = component_inverse_link_jet(component, eta);
6161 assert_eq!(integrated.mode, IntegratedExpectationMode::ExactClosedForm);
6162 assert_eq!(integrated.mean, point.mu);
6163 assert_eq!(integrated.d1, point.d1);
6164 assert_eq!(integrated.d2, point.d2);
6165 assert_eq!(integrated.d3, point.d3);
6166 }
6167 }
6168
6169 #[test]
6170 fn test_cloglog_controlled_matches_mathematical_target_on_small_sigma_grid() {
6171 let ctx = QuadratureContext::new();
6172 let cases = [
6176 (-30.0, 1e-10),
6177 (-30.0, 0.1),
6178 (-10.0, 0.24),
6179 (-3.0, 0.2),
6180 (0.0, 0.05),
6181 (0.4, 0.1),
6182 (3.0, 0.24),
6183 (10.0, 0.1),
6184 (30.0, 0.24),
6185 ];
6186
6187 for &(mu, sigma) in &cases {
6188 let approx = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
6189 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
6190 assert_relative_eq!(
6191 approx.mean,
6192 expected_mean,
6193 epsilon = 1e-12,
6194 max_relative = 2e-3
6195 );
6196 assert_relative_eq!(
6197 approx.dmean_dmu,
6198 expected_deriv,
6199 epsilon = 1e-12,
6200 max_relative = 4e-3
6201 );
6202 }
6203 }
6204
6205 #[test]
6206 fn test_cloglog_dispatch_uses_gamma_backend_for_large_sigma_central_regime() {
6207 let ctx = QuadratureContext::new();
6208 let out =
6209 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, -0.2, 0.8)
6210 .expect("cloglog integrated inverse-link moments should evaluate");
6211 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6212 assert!(out.mean.is_finite());
6213 assert!(out.dmean_dmu.is_finite());
6214 assert!(out.dmean_dmu >= 0.0);
6215 }
6216
6217 #[test]
6218 fn test_cloglog_dispatch_uses_large_sigma_asymptotic_without_ghq() {
6219 let ctx = QuadratureContext::new();
6220 let out =
6221 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, 0.0, 20.0)
6222 .expect("cloglog integrated inverse-link moments should evaluate");
6223 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6224 assert!(out.mean.is_finite());
6225 assert!(out.dmean_dmu.is_finite());
6226 assert!(out.dmean_dmu >= 0.0);
6227 }
6228
6229 #[test]
6230 fn test_cloglog_cc_matches_gamma_reference_on_central_case() {
6231 let ctx = QuadratureContext::new();
6232 let mu = -0.2;
6233 let sigma = 0.8;
6234 let cc = cloglog_survival_cc(&ctx, mu, sigma, CLOGLOG_CC_TOL).expect("cc backend");
6235 let gamma = cloglog_survival_gamma_reference(mu, sigma).expect("gamma backend");
6236 assert_relative_eq!(cc, gamma, epsilon = 5e-6, max_relative = 5e-6);
6237 }
6238
6239 #[test]
6240 fn test_cloglog_gamma_reference_matches_seeded_monte_carlo_small_case() {
6241 let mu = -0.2;
6242 let sigma = 0.8;
6243 let gamma =
6244 cloglog_posterior_meanwith_deriv_gamma_reference(mu, sigma).expect("gamma reference");
6245 let mut rng_state = 0x9e3779b97f4a7c15u64;
6246 let mut mean_mc = 0.0f64;
6247 let mut deriv_mc = 0.0f64;
6248 let n_samples = 300_000usize;
6249 for _ in 0..n_samples {
6250 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6251 let u1 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6252 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6253 let u2 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6254 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
6255 let eta = mu + sigma * z;
6256 mean_mc += cloglog_mean_exact(eta);
6257 deriv_mc += cloglog_mean_d1_exact(eta);
6258 }
6259 mean_mc /= n_samples as f64;
6260 deriv_mc /= n_samples as f64;
6261 assert_relative_eq!(gamma.mean, mean_mc, epsilon = 2e-3, max_relative = 2e-3);
6262 assert_relative_eq!(
6263 gamma.dmean_dmu,
6264 deriv_mc,
6265 epsilon = 2e-3,
6266 max_relative = 2e-3
6267 );
6268 }
6269
6270 #[test]
6271 fn test_logit_dispatch_uses_tail_asymptotic_outside_old_guard() {
6272 let ctx = QuadratureContext::new();
6273 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 35.0, 1.0)
6274 .expect("logit integrated inverse-link moments should evaluate");
6275 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6276 assert!(out.mean.is_finite());
6277 assert!(out.dmean_dmu.is_finite());
6278 assert!(out.dmean_dmu >= 0.0);
6279 }
6280
6281 #[test]
6282 fn test_logit_dispatch_prefers_erfcx_in_moderate_regime() {
6283 let ctx = QuadratureContext::new();
6294 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 1.1, 0.8)
6295 .expect("logit integrated inverse-link moments should evaluate");
6296 assert!(matches!(
6297 out.mode,
6298 IntegratedExpectationMode::ExactSpecialFunction
6299 | IntegratedExpectationMode::QuadratureFallback
6300 ));
6301 assert!(out.mean.is_finite());
6302 assert!(out.dmean_dmu.is_finite());
6303 assert!(out.dmean_dmu >= 0.0);
6304 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6305 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6306 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6307 }
6308
6309 #[test]
6310 fn test_logit_dispatch_large_sigma_uses_accurate_quadrature_not_monahan() {
6311 let ctx = QuadratureContext::new();
6320 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 0.5, 20.0)
6321 .expect("logit integrated inverse-link moments should evaluate");
6322 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
6323 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(0.5, 20.0);
6324 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6325 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6326 let kappa = (1.0 + std::f64::consts::PI * 20.0 * 20.0 / 8.0)
6329 .sqrt()
6330 .recip();
6331 let monahan_mean = gam_math::probability::normal_cdf(0.5 * kappa);
6332 assert!(
6333 (out.mean - monahan_mean).abs() > 1e-3,
6334 "dispatcher must not return the inaccurate Monahan mean {monahan_mean}; got {}",
6335 out.mean
6336 );
6337 }
6338
6339 #[test]
6340 fn test_logit_controlled_path_keeps_exact_backend_in_moderate_regime() {
6341 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("logit controlled");
6351 assert!(matches!(
6352 out.mode,
6353 IntegratedExpectationMode::ExactSpecialFunction
6354 | IntegratedExpectationMode::QuadratureFallback
6355 ));
6356 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6357 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6358 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6359 }
6360
6361 #[test]
6362 fn test_logit_dispatch_derivative_correct_at_mu_zero_small_sigma() {
6363 let ctx = QuadratureContext::new();
6372 for &(mu, sigma) in &[(0.0, 0.3), (0.0, 0.4), (0.0, 0.5)] {
6373 let out =
6374 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6375 .expect("logit integrated inverse-link moments should evaluate");
6376 assert_relative_eq!(out.mean, 0.5, epsilon = 1e-10);
6378 assert!(
6380 out.dmean_dmu <= 0.25 + 1e-9,
6381 "E[sigmoid'] must not exceed 0.25 at (μ={mu}, σ={sigma}); got {}",
6382 out.dmean_dmu
6383 );
6384 let (_, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6385 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-6);
6386 }
6387 }
6388
6389 #[test]
6390 fn test_logit_erfcx_exact_branch_is_self_certified() {
6391 for &(mu, sigma) in &[(8.0, 1.0), (10.0, 1.0), (15.0, 2.0)] {
6398 let out = logit_posterior_meanwith_deriv_exact(mu, sigma)
6399 .expect("erfcx branch should certify");
6400 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6401 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6402 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6403 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6404 }
6405 assert!(
6409 logit_posterior_meanwith_deriv_exact(0.0, 0.3).is_err(),
6410 "erfcx branch must not claim ExactSpecialFunction when it cannot certify the derivative"
6411 );
6412 }
6413
6414 #[test]
6415 fn test_logit_integrated_derivative_is_even_in_mu() {
6416 let ctx = QuadratureContext::new();
6422 for &(mu, sigma) in &[(0.3, 0.3), (1.1, 0.8), (10.0, 1.0), (3.0, 3.0), (35.0, 1.0)] {
6423 let pos =
6424 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6425 .expect("logit moments (+μ)");
6426 let neg =
6427 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, -mu, sigma)
6428 .expect("logit moments (-μ)");
6429 assert_relative_eq!(
6430 pos.dmean_dmu,
6431 neg.dmean_dmu,
6432 epsilon = 1e-9,
6433 max_relative = 1e-7
6434 );
6435 assert_relative_eq!(
6437 neg.mean,
6438 1.0 - pos.mean,
6439 epsilon = 1e-9,
6440 max_relative = 1e-7
6441 );
6442 }
6443 }
6444
6445 #[test]
6446 fn test_logit_dmean_dmu_equals_fd_of_mean_across_regimes() {
6447 let ctx = QuadratureContext::new();
6462 let h = 1e-4;
6463 let cases = [
6464 (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), ];
6474 for &(mu, sigma) in &cases {
6475 let at = |m: f64| {
6476 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, m, sigma)
6477 .expect("logit moments")
6478 };
6479 let out = at(mu);
6480 let fd = (at(mu + h).mean - at(mu - h).mean) / (2.0 * h);
6481 assert!(
6482 (out.dmean_dmu - fd).abs() <= 1e-5,
6483 "dmean_dmu must equal d/dμ of mean at (μ={mu}, σ={sigma}): \
6484 returned {}, FD of mean {} (mode {:?})",
6485 out.dmean_dmu,
6486 fd,
6487 out.mode
6488 );
6489 assert!(
6493 out.dmean_dmu <= 0.25 + 1e-9 && out.dmean_dmu >= 0.0,
6494 "dmean_dmu out of [0, 0.25] at (μ={mu}, σ={sigma}): {}",
6495 out.dmean_dmu
6496 );
6497 }
6498 }
6499
6500 #[test]
6501 fn test_logit_scalar_matches_jet_at_large_sigma() {
6502 let ctx = QuadratureContext::new();
6508 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0)] {
6509 let scalar =
6510 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6511 .expect("scalar logit moments");
6512 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6513 .expect("jet logit moments");
6514 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6518 assert_relative_eq!(scalar.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-8);
6519 assert_relative_eq!(
6520 scalar.dmean_dmu,
6521 ref_d1,
6522 epsilon = 1e-9,
6523 max_relative = 1e-8
6524 );
6525 assert_relative_eq!(scalar.mean, jet.mean, epsilon = 1e-12, max_relative = 1e-12);
6533 assert_relative_eq!(
6534 scalar.dmean_dmu,
6535 jet.d1,
6536 epsilon = 1e-12,
6537 max_relative = 1e-12
6538 );
6539 }
6540 }
6541
6542 #[test]
6543 fn test_logit_jet_accurate_at_wide_sigma() {
6544 let ctx = QuadratureContext::new();
6552 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0), (0.5, 20.0)] {
6553 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6554 .expect("wide-σ logit jet");
6555 let (rm, rd1, rd2, rd3) = logit_reference_jet_highres_simpson(mu, sigma);
6556 assert_relative_eq!(jet.mean, rm, epsilon = 1e-8, max_relative = 1e-7);
6557 assert_relative_eq!(jet.d1, rd1, epsilon = 1e-8, max_relative = 1e-6);
6558 assert_relative_eq!(jet.d2, rd2, epsilon = 1e-8, max_relative = 1e-6);
6559 assert_relative_eq!(jet.d3, rd3, epsilon = 1e-8, max_relative = 1e-6);
6560 let scalar =
6562 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6563 .expect("scalar logit moments");
6564 assert_relative_eq!(jet.d1, scalar.dmean_dmu, epsilon = 1e-12);
6565 assert_relative_eq!(jet.mean, scalar.mean, epsilon = 1e-12);
6566 let pirls = integrated_logit_inverse_link_jet_pirls(&ctx, mu, sigma)
6568 .expect("wide-σ PIRLS logit jet");
6569 assert_relative_eq!(pirls.mean, jet.mean, epsilon = 1e-12);
6570 assert_relative_eq!(pirls.d1, jet.d1, epsilon = 1e-12);
6571 assert_relative_eq!(pirls.d2, jet.d2, epsilon = 1e-12);
6572 assert_relative_eq!(pirls.d3, jet.d3, epsilon = 1e-12);
6573 assert_eq!(pirls.mode, jet.mode);
6574 }
6575 }
6576
6577 #[test]
6578 fn test_logit_jet_continuous_across_ghq_simpson_seam() {
6579 let ctx = QuadratureContext::new();
6587 let sigma = LOGIT_JET_GHQ_SIGMA_MAX;
6588 for mu in [-2.0, -0.5, 0.0, 0.7, 1.3, 3.0] {
6589 let ghq = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6591 .expect("jet at seam (GHQ dispatch)");
6592 let simpson = logit_wide_sigma_jet(mu, sigma).expect("jet at seam (Simpson)");
6594 assert_relative_eq!(ghq.mean, simpson.mean, epsilon = 1e-9, max_relative = 1e-8);
6597 assert_relative_eq!(ghq.d1, simpson.d1, epsilon = 1e-9, max_relative = 1e-7);
6598 assert_relative_eq!(ghq.d2, simpson.d2, epsilon = 1e-9, max_relative = 1e-7);
6599 assert_relative_eq!(ghq.d3, simpson.d3, epsilon = 1e-8, max_relative = 1e-6);
6600 }
6601 }
6602
6603 #[test]
6604 fn test_logit_batch_uses_same_dispatchvalues() {
6605 let ctx = QuadratureContext::new();
6606 let eta = ndarray::array![-2.0, 0.0, 1.25, 35.0];
6607 let se = ndarray::array![0.1, 0.5, 1.0, 1.0];
6608 let batch_mean = logit_posterior_mean_batch(&ctx, &eta, &se)
6609 .expect("logit posterior mean batch should evaluate");
6610 let (batchmu, batch_dmu) = logit_posterior_meanwith_deriv_batch(&ctx, &eta, &se)
6611 .expect("logit posterior mean derivative batch should evaluate");
6612 for i in 0..eta.len() {
6613 let direct = integrated_inverse_link_mean_and_derivative(
6614 &ctx,
6615 LinkFunction::Logit,
6616 eta[i],
6617 se[i],
6618 )
6619 .expect("logit integrated inverse-link moments should evaluate");
6620 assert_relative_eq!(batch_mean[i], direct.mean, epsilon = 1e-12);
6621 assert_relative_eq!(batchmu[i], direct.mean, epsilon = 1e-12);
6622 assert_relative_eq!(batch_dmu[i], direct.dmean_dmu, epsilon = 1e-12);
6623 }
6624 }
6625
6626 #[test]
6627 fn exact_logit_small_se_branch_loses_tail_derivative() {
6628 let eta = 50.0_f64;
6629 let stable_z = (-eta).exp();
6630 let stable_dmu = stable_z / (1.0_f64 + stable_z).powi(2);
6631 assert!(stable_dmu > 0.0);
6632 let out = logit_posterior_meanwith_deriv_exact(eta, 0.0).expect("exact branch");
6633 let dmu = out.dmean_dmu;
6634 assert!(
6635 (dmu - stable_dmu).abs() < 1e-30,
6636 "exact logit small-se branch should use the stable derivative z/(1+z)^2 at eta={eta}; got {} vs {}",
6637 dmu,
6638 stable_dmu
6639 );
6640 }
6641
6642 #[test]
6643 fn integrated_family_moments_rejects_latent_cloglog_without_concrete_handler() {
6644 let ctx = QuadratureContext::new();
6650 let latent =
6651 gam_problem::types::LatentCLogLogState::new(0.4).expect("valid latent cloglog state");
6652 let spec =
6653 LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::LatentCLogLog(latent));
6654 let likelihood = GlmLikelihoodSpec::canonical(spec);
6655 let err = integrated_family_moments_jet(
6656 &ctx,
6657 &likelihood,
6658 0.2,
6659 0.5,
6660 )
6661 .expect_err("latent cloglog moments should error in this dispatcher");
6662 assert!(format!("{err}").contains("LatentCLogLog"));
6663 }
6664
6665 #[test]
6666 fn integrated_family_moments_supports_stateful_sas() {
6667 let ctx = QuadratureContext::new();
6668 let sas = crate::mixture_link::state_from_sasspec(gam_problem::types::SasLinkSpec {
6669 initial_epsilon: 0.3,
6670 initial_log_delta: -0.2,
6671 })
6672 .expect("sas state should reconstruct from raw parameters");
6673 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Sas(sas));
6674 let likelihood = GlmLikelihoodSpec::canonical(spec);
6675 let out = integrated_family_moments_jet(
6676 &ctx,
6677 &likelihood,
6678 0.2,
6679 0.5,
6680 )
6681 .expect("stateful SAS integrated moments should evaluate");
6682 assert!(out.mean.is_finite());
6683 assert!(out.d1.is_finite());
6684 assert!(out.d2.is_finite());
6685 assert!(out.d3.is_finite());
6686 assert!(out.mean > 0.0 && out.mean < 1.0);
6687 }
6688
6689 #[test]
6690 fn integrated_family_moments_supports_pure_probit_mixture() {
6691 let ctx = QuadratureContext::new();
6692 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6693 components: vec![gam_problem::types::LinkComponent::Probit],
6694 initial_rho: ndarray::Array1::<f64>::zeros(0),
6695 })
6696 .expect("single-component probit mixture state");
6697 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
6698 let likelihood = GlmLikelihoodSpec::canonical(spec);
6699 let out = integrated_family_moments_jet(
6700 &ctx,
6701 &likelihood,
6702 0.7,
6703 1.3,
6704 )
6705 .expect("pure probit mixture integrated moments should evaluate");
6706 let exact = integrated_probit_jet(0.7, 1.3);
6707 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
6708 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
6709 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
6710 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
6711 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
6712 }
6713
6714 #[test]
6715 fn integrated_family_moments_supports_pure_logit_mixture() {
6716 let ctx = QuadratureContext::new();
6717 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6718 components: vec![gam_problem::types::LinkComponent::Logit],
6719 initial_rho: ndarray::Array1::<f64>::zeros(0),
6720 })
6721 .expect("single-component logit mixture state");
6722 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
6723 let likelihood = GlmLikelihoodSpec::canonical(spec);
6724 let out = integrated_family_moments_jet(
6725 &ctx,
6726 &likelihood,
6727 1.1,
6728 0.8,
6729 )
6730 .expect("pure logit mixture integrated moments should evaluate");
6731 let exact = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, 1.1, 0.8)
6732 .expect("canonical integrated logit jet");
6733 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
6734 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
6735 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
6736 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
6737 assert_eq!(out.mode, exact.mode);
6738 }
6739
6740 #[test]
6741 fn integrated_family_moments_supports_stateful_mixture() {
6742 let ctx = QuadratureContext::new();
6743 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6744 components: vec![
6745 gam_problem::types::LinkComponent::Logit,
6746 gam_problem::types::LinkComponent::Probit,
6747 ],
6748 initial_rho: ndarray::array![0.35],
6749 })
6750 .expect("mixture state should reconstruct from rho");
6751 let spec = LikelihoodSpec::new(
6752 ResponseFamily::Binomial,
6753 InverseLink::Mixture(state.clone()),
6754 );
6755 let likelihood = GlmLikelihoodSpec::canonical(spec);
6756 let out = integrated_family_moments_jet(
6757 &ctx,
6758 &likelihood,
6759 0.2,
6760 0.5,
6761 )
6762 .expect("stateful mixture integrated moments should evaluate");
6763 let direct = integrated_mixture_jet(&ctx, 0.2, 0.5, &state)
6764 .expect("direct integrated mixture jet should evaluate");
6765 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
6766 assert_relative_eq!(out.d1, direct.d1, epsilon = 1e-12);
6767 assert_relative_eq!(out.d2, direct.d2, epsilon = 1e-12);
6768 assert_relative_eq!(out.d3, direct.d3, epsilon = 1e-12);
6769 assert_eq!(out.mode, direct.mode);
6770 }
6771
6772 #[test]
6773 fn integrated_family_moments_use_scale_dispersion_for_tweedie_and_gamma() {
6774 let ctx = QuadratureContext::new();
6778 let e = 0.3_f64;
6780 let se = 0.5_f64;
6781 let m = (e + 0.5 * se * se).exp();
6782
6783 let p = 1.5_f64;
6785 let phi = 2.0_f64;
6786 let tweedie = LikelihoodSpec::tweedie_log(p);
6787 let tweedie_likelihood = GlmLikelihoodSpec {
6788 spec: tweedie.clone(),
6789 scale: LikelihoodScaleMetadata::EstimatedTweediePhi { phi },
6790 };
6791 let out = integrated_family_moments_jet(
6792 &ctx,
6793 &tweedie_likelihood,
6794 e,
6795 se,
6796 )
6797 .expect("tweedie integrated moments should evaluate");
6798 let expected = phi * m.powf(p);
6799 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
6800 assert_relative_eq!(out.variance / m.powf(p), phi, epsilon = 1e-12);
6802
6803 let shape = 4.0_f64;
6805 let gamma = LikelihoodSpec::gamma_log();
6806 let gamma_likelihood = GlmLikelihoodSpec {
6807 spec: gamma.clone(),
6808 scale: LikelihoodScaleMetadata::EstimatedGammaShape { shape },
6809 };
6810 let out = integrated_family_moments_jet(
6811 &ctx,
6812 &gamma_likelihood,
6813 e,
6814 se,
6815 )
6816 .expect("gamma integrated moments should evaluate");
6817 let expected = m * m / shape;
6818 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
6819 assert_relative_eq!(out.variance / (m * m), 1.0 / shape, epsilon = 1e-12);
6821
6822 let poisson = LikelihoodSpec::poisson_log();
6824 let poisson_likelihood = GlmLikelihoodSpec::canonical(poisson);
6825 let out = integrated_family_moments_jet(
6826 &ctx,
6827 &poisson_likelihood,
6828 e,
6829 se,
6830 )
6831 .expect("poisson integrated moments should evaluate");
6832 assert_relative_eq!(out.variance, m, epsilon = 1e-12);
6833
6834 let theta = 3.0_f64;
6836 let nb = LikelihoodSpec::negative_binomial_log(theta);
6837 let nb_likelihood = GlmLikelihoodSpec::canonical(nb);
6838 let out = integrated_family_moments_jet(
6839 &ctx,
6840 &nb_likelihood,
6841 e,
6842 se,
6843 )
6844 .expect("negative-binomial integrated moments should evaluate");
6845 assert_relative_eq!(out.variance, m + m * m / theta, epsilon = 1e-12);
6846
6847 let missing_gamma = GlmLikelihoodSpec {
6849 spec: gamma,
6850 scale: LikelihoodScaleMetadata::Unspecified,
6851 };
6852 let err = integrated_family_moments_jet(
6853 &ctx,
6854 &missing_gamma,
6855 e,
6856 se,
6857 )
6858 .expect_err("gamma without a shape in the scale metadata must error");
6859 assert!(
6860 format!("{err}").contains("GammaShape"),
6861 "unexpected error message: {err}"
6862 );
6863
6864 let missing_tweedie = GlmLikelihoodSpec {
6866 spec: tweedie,
6867 scale: LikelihoodScaleMetadata::Unspecified,
6868 };
6869 let err = integrated_family_moments_jet(
6870 &ctx,
6871 &missing_tweedie,
6872 e,
6873 se,
6874 )
6875 .expect_err("tweedie without a φ in the scale metadata must error");
6876 assert!(
6877 format!("{err}").contains("EstimatedTweediePhi"),
6878 "unexpected error message: {err}"
6879 );
6880 }
6881
6882 #[test]
6885 fn cloglog_g_derivatives_at_zero() {
6886 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(0.0);
6887 let expected_g = 1.0 - (-1.0_f64).exp();
6889 assert_relative_eq!(g, expected_g, epsilon = 1e-14);
6890 let e_neg1 = (-1.0_f64).exp();
6892 assert_relative_eq!(g1, e_neg1, epsilon = 1e-14);
6893 assert_relative_eq!(g2, 0.0, epsilon = 1e-14);
6895 assert_relative_eq!(g3, -e_neg1, epsilon = 1e-14);
6897 assert_relative_eq!(g4, -e_neg1, epsilon = 1e-14);
6899 }
6900
6901 #[test]
6902 fn cloglog_g_derivatives_saturation() {
6903 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(50.0);
6905 assert_relative_eq!(g, 1.0, epsilon = 1e-10);
6906 assert_eq!(g1, 0.0);
6907 assert_eq!(g2, 0.0);
6908 assert_eq!(g3, 0.0);
6909 assert_eq!(g4, 0.0);
6910
6911 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(-50.0);
6913 let expected = (-50.0_f64).exp();
6914 assert_relative_eq!(g, expected, max_relative = 1e-10);
6915 assert_relative_eq!(g1, expected, max_relative = 1e-10);
6916 assert_relative_eq!(g2, expected, max_relative = 1e-10);
6918 assert_relative_eq!(g3, expected, max_relative = 1e-10);
6919 assert_relative_eq!(g4, expected, max_relative = 1e-10);
6920 }
6921
6922 #[test]
6923 fn cloglog_ghq_value_sigma_zero_matches_pointwise() {
6924 let ctx = QuadratureContext::new();
6925 for &mu in &[-2.0, -1.0, 0.0, 0.5, 1.5] {
6927 let val = cloglog_ghq_value(&ctx, mu, 0.0, 21);
6928 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
6929 assert_relative_eq!(val, g, epsilon = 1e-14);
6930 }
6931 }
6932
6933 #[test]
6934 fn cloglog_ghq_value_bounded_zero_one() {
6935 let ctx = QuadratureContext::new();
6936 for &mu in &[-5.0, -2.0, 0.0, 1.0, 3.0, 10.0] {
6938 for &sigma in &[0.1, 0.5, 1.0, 2.0, 5.0] {
6939 let val = cloglog_ghq_value(&ctx, mu, sigma, 31);
6940 assert!((0.0..=1.0).contains(&val), "L({mu},{sigma}) = {val}");
6941 }
6942 }
6943 }
6944
6945 #[test]
6946 fn cloglog_ghq_derivatives_sigma_zero_matches_pointwise() {
6947 let ctx = QuadratureContext::new();
6948 let mu = 0.3;
6949 let d = cloglog_ghq_derivatives(&ctx, mu, 0.0, 21);
6950 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
6951 assert_relative_eq!(d.l, g, epsilon = 1e-14);
6952 assert_relative_eq!(d.l_mu, g1, epsilon = 1e-14);
6953 assert_relative_eq!(d.l_mumu, g2, epsilon = 1e-14);
6954 assert_relative_eq!(d.l_mumumu, g3, epsilon = 1e-14);
6955 assert_relative_eq!(d.l_mumumumu, g4, epsilon = 1e-14);
6956
6957 assert_eq!(d.l_sigma, 0.0);
6959 assert_eq!(d.l_musigma, 0.0);
6960 assert_eq!(d.l_mumusigma, 0.0);
6961 assert_eq!(d.l_mumumusigma, 0.0);
6962 assert_eq!(d.l_sigmasigmasigma, 0.0);
6963 assert_eq!(d.l_musigmasigmasigma, 0.0);
6964
6965 assert_relative_eq!(d.l_sigmasigma, g2, epsilon = 1e-14);
6968 assert_relative_eq!(d.l_musigmasigma, g3, epsilon = 1e-14);
6969 assert_relative_eq!(d.l_mumusigmasigma, g4, epsilon = 1e-14);
6970 assert_relative_eq!(d.l_sigmasigmasigmasigma, 3.0 * g4, epsilon = 1e-14);
6971 }
6972
6973 #[test]
6974 fn cloglog_ghq_derivatives_finite_difference_mu() {
6975 let ctx = QuadratureContext::new();
6977 let mu = 0.5;
6978 let sigma = 0.8;
6979 let h = 1e-6;
6980 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
6981 let l_plus = cloglog_ghq_value(&ctx, mu + h, sigma, 31);
6982 let l_minus = cloglog_ghq_value(&ctx, mu - h, sigma, 31);
6983 let fd_mu = (l_plus - l_minus) / (2.0 * h);
6984 assert_relative_eq!(d.l_mu, fd_mu, epsilon = 1e-5);
6985
6986 let d_plus = cloglog_ghq_derivatives(&ctx, mu + h, sigma, 31);
6988 let d_minus = cloglog_ghq_derivatives(&ctx, mu - h, sigma, 31);
6989 let fd_mumu = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
6990 assert_relative_eq!(d.l_mumu, fd_mumu, epsilon = 1e-4);
6991 }
6992
6993 #[test]
6994 fn cloglog_ghq_derivatives_finite_difference_sigma() {
6995 let ctx = QuadratureContext::new();
6997 let mu = 0.2;
6998 let sigma = 1.0;
6999 let h = 1e-6;
7000 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
7001 let l_plus = cloglog_ghq_value(&ctx, mu, sigma + h, 31);
7002 let l_minus = cloglog_ghq_value(&ctx, mu, sigma - h, 31);
7003 let fd_sigma = (l_plus - l_minus) / (2.0 * h);
7004 assert_relative_eq!(d.l_sigma, fd_sigma, epsilon = 1e-5);
7005 }
7006
7007 #[test]
7008 fn cloglog_ghq_derivatives_finite_difference_cross() {
7009 let ctx = QuadratureContext::new();
7011 let mu = -0.5;
7012 let sigma = 0.6;
7013 let h = 1e-6;
7014 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
7015 let d_plus = cloglog_ghq_derivatives(&ctx, mu, sigma + h, 31);
7016 let d_minus = cloglog_ghq_derivatives(&ctx, mu, sigma - h, 31);
7017 let fd_musigma = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
7018 assert_relative_eq!(d.l_musigma, fd_musigma, epsilon = 1e-4);
7019 }
7020
7021 #[test]
7022 fn cloglog_ghq_l_mu_nonnegative() {
7023 let ctx = QuadratureContext::new();
7025 for &mu in &[-3.0, -1.0, 0.0, 1.0, 3.0] {
7026 for &sigma in &[0.1, 0.5, 1.0, 2.0] {
7027 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 21);
7028 assert!(
7029 d.l_mu >= -1e-14,
7030 "L_mu should be non-negative at mu={mu}, sigma={sigma}: got {}",
7031 d.l_mu
7032 );
7033 }
7034 }
7035 }
7036
7037 #[test]
7038 fn cloglog_ghq_adaptive_matches_explicit() {
7039 let ctx = QuadratureContext::new();
7040 let mu = 0.7;
7041 let sigma = 1.2;
7042 let adaptive = cloglog_ghq_derivatives_adaptive(&ctx, mu, sigma);
7043 let n = adaptive_point_count_from_sd(sigma);
7044 let explicit = cloglog_ghq_derivatives(&ctx, mu, sigma, n);
7045 assert_relative_eq!(adaptive.l, explicit.l, epsilon = 1e-15);
7046 assert_relative_eq!(adaptive.l_mu, explicit.l_mu, epsilon = 1e-15);
7047 assert_relative_eq!(adaptive.l_sigma, explicit.l_sigma, epsilon = 1e-15);
7048 assert_relative_eq!(adaptive.l_mumu, explicit.l_mumu, epsilon = 1e-15);
7049 }
7050
7051 #[test]
7052 fn cloglog_ghq_value_matches_mathematical_target_in_central_regime() {
7053 let ctx = QuadratureContext::new();
7054 for &mu in &[-1.0, 0.0, 0.5, 2.0] {
7055 for &sigma in &[0.1, 0.5, 1.0] {
7056 let ghq = cloglog_ghq_value(&ctx, mu, sigma, 51);
7057 let (expected_mean, _) = cloglog_reference_mean_and_derivative(mu, sigma);
7058 assert_relative_eq!(ghq, expected_mean, epsilon = 1e-12, max_relative = 2e-8);
7059 }
7060 }
7061 }
7062
7063 #[test]
7066 fn cloglog_negative_tail_mean_matches_exact_near_transition() {
7067 let eta: f64 = -30.0;
7071 let exact = {
7072 let ex = eta.exp();
7073 -(-ex).exp_m1()
7074 };
7075 let tail = cloglog_negative_tail_mean(eta);
7076 assert!(
7077 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7078 "tail mean at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7079 );
7080 }
7081
7082 #[inline]
7083 fn cloglog_negative_tail_derivative(eta: f64) -> f64 {
7084 if eta < -745.0 {
7086 0.0
7087 } else {
7088 let ex = safe_exp(eta);
7089 (ex * (-ex).exp()).max(0.0)
7090 }
7091 }
7092
7093 #[test]
7094 fn cloglog_negative_tail_derivative_matches_exact_near_transition() {
7095 let eta: f64 = -30.0;
7097 let ex = eta.exp();
7098 let exact = ex * (-ex).exp();
7099 let tail = cloglog_negative_tail_derivative(eta);
7100 assert!(
7101 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7102 "tail derivative at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7103 );
7104 }
7105
7106 #[test]
7107 fn cloglog_negative_tail_degenerate_branch_matches_target_near_transition() {
7108 let ctx = QuadratureContext::default();
7109 let sigma = 0.0;
7110 for &mu in &[-30.001, -30.0, -29.999] {
7111 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7112 assert_relative_eq!(
7113 out.mean,
7114 cloglog_mean_exact(mu),
7115 epsilon = 1e-28,
7116 max_relative = 1e-15
7117 );
7118 assert_relative_eq!(
7119 out.dmean_dmu,
7120 cloglog_mean_d1_exact(mu),
7121 epsilon = 1e-28,
7122 max_relative = 1e-15
7123 );
7124 }
7125 }
7126
7127 #[test]
7128 fn cloglog_negative_tail_small_sigma_branch_matches_target_near_transition() {
7129 let ctx = QuadratureContext::default();
7130 let sigma = 0.1;
7131 for &mu in &[-30.001, -30.0, -29.999] {
7132 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7133 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
7134 assert_relative_eq!(
7135 out.mean,
7136 expected_mean,
7137 epsilon = 1e-24,
7138 max_relative = 1e-10
7139 );
7140 assert_relative_eq!(
7141 out.dmean_dmu,
7142 expected_deriv,
7143 epsilon = 1e-24,
7144 max_relative = 1e-10
7145 );
7146 }
7147 }
7148
7149 fn ref_cholesky_heap(cov: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
7153 let n = cov.len();
7154 if n == 0 || cov.iter().any(|r| r.len() != n) {
7155 return None;
7156 }
7157 let mut base = cov.to_vec();
7158 for retry in 0..8 {
7159 let jitter = if retry == 0 {
7160 0.0
7161 } else {
7162 1e-12 * 10f64.powi(retry - 1)
7163 };
7164 if jitter > 0.0 {
7165 for i in 0..n {
7166 base[i][i] = cov[i][i] + jitter;
7167 }
7168 }
7169 let mut l = vec![vec![0.0_f64; n]; n];
7170 let mut ok = true;
7171 for i in 0..n {
7172 for j in 0..=i {
7173 let mut sum = base[i][j];
7174 for k in 0..j {
7175 sum -= l[i][k] * l[j][k];
7176 }
7177 if i == j {
7178 if !sum.is_finite() || sum <= 0.0 {
7179 ok = false;
7180 break;
7181 }
7182 l[i][j] = sum.sqrt();
7183 } else {
7184 l[i][j] = sum / l[j][j];
7185 }
7186 }
7187 if !ok {
7188 break;
7189 }
7190 }
7191 if ok {
7192 return Some(l);
7193 }
7194 }
7195 None
7196 }
7197
7198 #[test]
7199 fn cholesky_static_matches_heap_d2() {
7200 let cases: &[[[f64; 2]; 2]] = &[
7203 [[1.0, 0.0], [0.0, 1.0]],
7204 [[2.5, 0.3], [0.3, 0.75]],
7205 [[1.0, 0.9999], [0.9999, 1.0]],
7206 [[1e-10, 0.0], [0.0, 1e-10]],
7207 [[4.0, -1.5], [-1.5, 2.25]],
7208 ];
7209 for cov in cases {
7210 let stack = cholesky_static_with_jitter::<2>(cov).expect("stack cholesky");
7211 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7212 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7213 for i in 0..2 {
7214 for j in 0..2 {
7215 assert_eq!(
7216 stack[i][j].to_bits(),
7217 heap[i][j].to_bits(),
7218 "mismatch at ({i},{j}) for cov={cov:?}"
7219 );
7220 }
7221 }
7222 }
7223 }
7224
7225 #[test]
7226 fn cholesky_static_matches_heap_d3() {
7227 let cases: &[[[f64; 3]; 3]] = &[
7228 [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
7229 [[2.0, 0.5, 0.1], [0.5, 1.5, -0.2], [0.1, -0.2, 0.8]],
7230 [[4.0, 1.0, 0.5], [1.0, 3.0, 0.25], [0.5, 0.25, 2.0]],
7231 ];
7232 for cov in cases {
7233 let stack = cholesky_static_with_jitter::<3>(cov).expect("stack cholesky");
7234 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7235 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7236 for i in 0..3 {
7237 for j in 0..3 {
7238 assert_eq!(
7239 stack[i][j].to_bits(),
7240 heap[i][j].to_bits(),
7241 "mismatch at ({i},{j}) for cov={cov:?}"
7242 );
7243 }
7244 }
7245 }
7246 }
7247
7248 #[test]
7249 fn cholesky_static_d1() {
7250 let l = cholesky_static_with_jitter::<1>(&[[2.25]]).expect("d=1");
7251 assert_eq!(l[0][0], 1.5);
7252 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e-13]]).is_some());
7262 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e3]]).is_none());
7265 }
7266}