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::quadrature::{GaussHermiteRule, gauss_hermite_rule};
182use gam_math::special::stable_polynomial_times_exp_neg as cloglog_stable_poly_times_exp_neg;
183use gam_problem::types::{
184 GlmLikelihoodSpec, InverseLink, LinkComponent, LinkFunction, MixtureLinkState, ResponseFamily,
185 SasLinkState, StandardLink,
186};
187const N_POINTS: usize = 7;
189const SQRT_2: f64 = std::f64::consts::SQRT_2;
190const QUADRATURE_EXP_LOG_MAX: f64 = 700.0;
191
192#[inline]
195fn safe_exp(x: f64) -> f64 {
196 if x.is_nan() {
197 f64::NAN
198 } else {
199 x.min(QUADRATURE_EXP_LOG_MAX).exp()
200 }
201}
202
203#[inline]
204fn safe_expwith_saturation(x: f64) -> (f64, bool) {
205 (safe_exp(x), x > QUADRATURE_EXP_LOG_MAX)
206}
207
208#[derive(Clone, Copy, Debug, Default)]
209struct Complex {
210 re: f64,
211 im: f64,
212}
213
214pub struct QuadratureContext {
216 gh_cache: OnceLock<GaussHermiteRule>,
217 gh15_cache: OnceLock<GaussHermiteRule>,
218 gh21_cache: OnceLock<GaussHermiteRule>,
219 gh31_cache: OnceLock<GaussHermiteRule>,
220 gh51_cache: OnceLock<GaussHermiteRule>,
221 cc_cache: Mutex<HashMap<usize, Arc<ClenshawCurtisRule>>>,
225}
226
227#[derive(Clone, Copy, Debug, Eq, PartialEq)]
228pub enum IntegratedExpectationMode {
229 ExactClosedForm,
230 ExactSpecialFunction,
231 ControlledAsymptotic,
232 QuadratureFallback,
233}
234
235impl IntegratedExpectationMode {
236 #[inline]
240 pub const fn rank(self) -> u8 {
241 match self {
242 Self::ExactClosedForm => 0,
243 Self::ExactSpecialFunction => 1,
244 Self::ControlledAsymptotic => 2,
245 Self::QuadratureFallback => 3,
246 }
247 }
248}
249
250#[derive(Clone, Copy, Debug)]
251pub struct IntegratedMeanDerivative {
252 pub mean: f64,
253 pub dmean_dmu: f64,
254 pub mode: IntegratedExpectationMode,
255}
256
257#[derive(Clone, Copy, Debug)]
258pub struct IntegratedInverseLinkJet {
259 pub mean: f64,
260 pub d1: f64,
261 pub d2: f64,
262 pub d3: f64,
263 pub mode: IntegratedExpectationMode,
264}
265
266#[derive(Clone, Copy, Debug)]
267pub(crate) struct IntegratedInverseLinkJet5 {
268 pub mean: f64,
269 pub d1: f64,
270 pub d2: f64,
271 pub d3: f64,
272 pub d4: f64,
273 pub d5: f64,
274 pub mode: IntegratedExpectationMode,
275}
276
277#[inline]
278pub(crate) fn validate_latent_cloglog_inputs(eta: f64, sigma: f64) -> Result<(), EstimationError> {
279 if !eta.is_finite() || !sigma.is_finite() || sigma < 0.0 {
280 crate::bail_invalid_estim!(
281 "latent cloglog jet requires finite eta and sigma >= 0, got eta={eta}, sigma={sigma}"
282 );
283 }
284 Ok::<(), _>(())
285}
286
287#[derive(Clone, Copy, Debug)]
292pub struct IntegratedMomentsJet {
293 pub mean: f64,
294 pub variance: f64,
295 pub d1: f64,
296 pub d2: f64,
297 pub d3: f64,
298 pub mode: IntegratedExpectationMode,
299}
300
301const LOGIT_SIGMA_DEGENERATE: f64 = 1e-10;
302const LOGIT_ERFCX_SIGMA_MIN: f64 = 2.5e-1;
303const LOGIT_TAIL_LOG_MAX: f64 = -18.0;
304const LOGIT_ERFCX_MU_MAX: f64 = 40.0;
305const LOGIT_ERFCX_SIGMA_MAX: f64 = 6.0;
306const LOGIT_JET_GHQ_SIGMA_MAX: f64 = 1.0;
330const CLOGLOG_SIGMA_DEGENERATE: f64 = 1e-10;
331const CLOGLOG_SIGMA_TAYLOR_MAX: f64 = 0.25;
332const CLOGLOG_JET_MOMENT_SIGMA_MAX: f64 = 1.0;
339const CLOGLOG_RARE_EVENT_LOG_MAX: f64 = -18.0;
340const CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN: f64 = 8.0;
341const CLOGLOG_POSITIVE_SATURATION_EDGE: f64 = 5.0;
342const CLOGLOG_POSITIVE_SATURATION_SIGMAS: f64 = 8.0;
343const LOG_SURVIVAL_PANEL_LOG_DROP: f64 = 60.0;
393const LOG_SURVIVAL_PANEL_ORDER_LOG_DROP: f64 = 2.0;
400const LOG_SURVIVAL_PANEL_ARCLENGTH_NODE_DENSITY: f64 = 4.6;
409const LOG_SURVIVAL_PANEL_SIGMA_NODE_SCALE: f64 = 70.0;
418const LOG_SURVIVAL_PANEL_ORDER_NODES: f64 = 8.0;
421const LOG_SURVIVAL_PANEL_MIN_NODES: usize = 65;
424const LOG_SURVIVAL_PANEL_MAX_NODES: usize = 4097;
428pub const LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER: usize = 8;
434const LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION: f64 = 6.1;
456const SERIES_CONSECUTIVE_SMALL_TERMS: usize = 6;
457const LOGIT_MAX_TERMS: usize = 160;
458const LOGIT_ERFCX_ACCURACY_TARGET: f64 = 1.0e-11;
474const CLOGLOG_MILES_ALPHA: f64 = 60.0;
475const CLOGLOG_MILES_MAX_TERMS: usize = 256;
476const CLOGLOG_MILES_PEAK_LOG_MAX: f64 = 0.0;
493const CLOGLOG_GAMMA_K_REF: f64 = 0.5;
494const CLOGLOG_GAMMA_T_MAX_REF: f64 = 24.0;
495const CLOGLOG_GAMMA_H_REF: f64 = 0.01;
496const CLOGLOG_CC_TOL: f64 = 1e-12;
500const CLOGLOG_CC_NODE_CAP: usize = 1025;
504const CLOGLOG_GAMMA_SAMPLE_COUNT: usize =
508 (CLOGLOG_GAMMA_T_MAX_REF / CLOGLOG_GAMMA_H_REF) as usize + 1;
509const CLOGLOG_CC_PREFER_THRESHOLD: usize = CLOGLOG_GAMMA_SAMPLE_COUNT / 3;
514const CLOGLOG_CC_MIN_N: usize = 17;
518
519impl QuadratureContext {
520 pub fn new() -> Self {
521 Self {
522 gh_cache: OnceLock::new(),
523 gh15_cache: OnceLock::new(),
524 gh21_cache: OnceLock::new(),
525 gh31_cache: OnceLock::new(),
526 gh51_cache: OnceLock::new(),
527 cc_cache: Mutex::new(HashMap::new()),
528 }
529 }
530
531 fn gauss_hermite(&self) -> &GaussHermiteRule {
532 self.gh_cache.get_or_init(compute_gauss_hermite)
533 }
534
535 fn gauss_hermite_n(&self, n: usize) -> &GaussHermiteRule {
536 match n {
537 7 => self.gh15_cache.get_or_init(|| compute_gauss_hermite_n(15)),
540 15 => self.gh15_cache.get_or_init(|| compute_gauss_hermite_n(15)),
541 21 => self.gh21_cache.get_or_init(|| compute_gauss_hermite_n(21)),
542 31 => self.gh31_cache.get_or_init(|| compute_gauss_hermite_n(31)),
543 51 => self.gh51_cache.get_or_init(|| compute_gauss_hermite_n(51)),
544 _ => self.gh21_cache.get_or_init(|| compute_gauss_hermite_n(21)),
545 }
546 }
547
548 fn clenshaw_curtis_n(&self, n: usize) -> Arc<ClenshawCurtisRule> {
549 let mut cache = match self.cc_cache.lock() {
550 Ok(guard) => guard,
551 Err(poisoned) => poisoned.into_inner(),
552 };
553 cache
554 .entry(n)
555 .or_insert_with(|| Arc::new(compute_clenshaw_curtis_n(n)))
556 .clone()
557 }
558}
559
560impl Default for QuadratureContext {
561 fn default() -> Self {
562 Self::new()
563 }
564}
565
566#[derive(Clone)]
567struct ClenshawCurtisRule {
568 nodes: Vec<f64>,
569 weights: Vec<f64>,
570}
571
572fn compute_clenshaw_curtis_n(n: usize) -> ClenshawCurtisRule {
573 assert!(
574 n >= 2,
575 "Clenshaw-Curtis rule requires at least two nodes: n={n}"
576 );
577 let m = n - 1;
592 let theta: Vec<f64> = (0..=m)
593 .map(|j| std::f64::consts::PI * (j as f64) / (m as f64))
594 .collect();
595 let nodes: Vec<f64> = theta.iter().map(|&th| th.cos()).collect();
596
597 if n == 2 {
598 return ClenshawCurtisRule {
599 nodes,
600 weights: vec![1.0, 1.0],
601 };
602 }
603
604 let mut weights = vec![0.0_f64; n];
605 let mut v = vec![1.0_f64; m - 1];
606
607 if m.is_multiple_of(2) {
608 let w0 = 1.0 / ((m * m - 1) as f64);
609 weights[0] = w0;
610 weights[m] = w0;
611 for k in 1..(m / 2) {
612 let denom = (4 * k * k - 1) as f64;
613 for j in 1..m {
614 v[j - 1] -= 2.0 * (2.0 * (k as f64) * theta[j]).cos() / denom;
615 }
616 }
617 for j in 1..m {
618 v[j - 1] -= ((m as f64) * theta[j]).cos() / ((m * m - 1) as f64);
619 }
620 } else {
621 let w0 = 1.0 / ((m * m) as f64);
622 weights[0] = w0;
623 weights[m] = w0;
624 for k in 1..=((m - 1) / 2) {
625 let denom = (4 * k * k - 1) as f64;
626 for j in 1..m {
627 v[j - 1] -= 2.0 * (2.0 * (k as f64) * theta[j]).cos() / denom;
628 }
629 }
630 }
631
632 for j in 1..m {
633 weights[j] = 2.0 * v[j - 1] / (m as f64);
634 }
635
636 for j in 0..=(m / 2) {
640 let jj = m - j;
641 let avg = 0.5 * (weights[j] + weights[jj]);
642 weights[j] = avg;
643 weights[jj] = avg;
644 }
645 let weight_sum: f64 = weights.iter().sum();
646 if weight_sum.is_finite() && weight_sum != 0.0 {
647 let scale = 2.0 / weight_sum;
648 for w in &mut weights {
649 *w *= scale;
650 }
651 }
652
653 ClenshawCurtisRule { nodes, weights }
654}
655
656fn cloglog_cc_required_nodes(mu: f64, sigma: f64, tol: f64) -> Result<usize, EstimationError> {
657 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0 && tol.is_finite() && tol > 0.0) {
658 crate::bail_invalid_estim!(
659 "CC cloglog backend requires finite mu, positive sigma, and positive tolerance"
660 .to_string(),
661 );
662 }
663
664 let p_tail = (tol / 8.0).clamp(1e-300, 0.25);
669 let a = gam_math::probability::standard_normal_quantile(p_tail)
670 .map(|z| -z)
671 .unwrap_or(8.0)
672 .max(1.0);
673
674 let ay = a * sigma;
675 let y = if ay > 0.0 {
676 1.0_f64.min(std::f64::consts::PI / (4.0 * ay))
677 } else {
678 1.0
679 };
680 let rho = y + (1.0 + y * y).sqrt();
681 let m_s = (0.5 * (a * y) * (a * y)).exp() / (2.0 * std::f64::consts::PI).sqrt();
682 let eps_quad = (tol / 4.0).max(1e-300);
683 let numer = ((8.0 * a * m_s) / ((rho - 1.0).max(1e-12) * eps_quad)).max(1.0);
684 let denom = rho.ln();
685 if !denom.is_finite() || denom <= 0.0 {
686 crate::bail_invalid_estim!("CC cloglog backend ellipse bound became degenerate");
687 }
688
689 let mut n = (1.0 + numer.ln() / denom).ceil() as usize;
690 n = n.max(CLOGLOG_CC_MIN_N);
691 if n.is_multiple_of(2) {
692 n += 1;
693 }
694 Ok(n)
695}
696
697#[inline]
698fn cloglog_should_prefer_cc(mu: f64, sigma: f64, tol: f64) -> bool {
699 match cloglog_cc_required_nodes(mu, sigma, tol) {
705 Ok(n) => n <= CLOGLOG_CC_PREFER_THRESHOLD,
706 Err(_) => false,
707 }
708}
709
710fn compute_gauss_hermite() -> GaussHermiteRule {
713 compute_gauss_hermite_n(N_POINTS)
714}
715
716pub(crate) fn compute_gauss_hermite_n(n: usize) -> GaussHermiteRule {
717 gauss_hermite_rule(n)
722 .unwrap_or_else(|error| panic!("shared Gauss-Hermite construction failed: {error}"))
723}
724
725#[inline]
736pub fn logit_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
737 match logit_posterior_meanwith_deriv_controlled(eta, se_eta) {
738 Ok(out) => out.mean,
739 Err(_) => integrate_normal_ghq_adaptive(ctx, eta, se_eta, sigmoid),
740 }
741}
742
743#[inline]
751pub fn logit_posterior_meanwith_deriv(
752 eta: f64,
753 se_eta: f64,
754) -> Result<(f64, f64), EstimationError> {
755 let out = logit_posterior_meanwith_deriv_controlled(eta, se_eta)?;
767 Ok((out.mean, out.dmean_dmu))
768}
769
770#[inline]
771pub fn probit_posterior_meanwith_deriv_exact(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
772 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 1e-12 {
800 let mean = gam_math::probability::normal_cdf(mu);
801 let dmean_dmu = gam_math::probability::normal_pdf(mu);
802 return IntegratedMeanDerivative {
803 mean,
804 dmean_dmu,
805 mode: IntegratedExpectationMode::ExactClosedForm,
806 };
807 }
808 let denom = (1.0 + sigma * sigma).sqrt();
809 let z = mu / denom;
810 IntegratedMeanDerivative {
811 mean: gam_math::probability::normal_cdf(z),
812 dmean_dmu: gam_math::probability::normal_pdf(z) / denom,
813 mode: IntegratedExpectationMode::ExactClosedForm,
814 }
815}
816
817#[inline]
818fn logistic_normal_exact_eligible(mu: f64, sigma: f64) -> bool {
819 mu.is_finite()
820 && sigma.is_finite()
821 && mu.abs() <= LOGIT_ERFCX_MU_MAX
822 && (LOGIT_ERFCX_SIGMA_MIN..=LOGIT_ERFCX_SIGMA_MAX).contains(&sigma)
823}
824
825#[inline]
869fn logistic_normal_series_cutoff(mu: f64, sigma: f64, target_accuracy: f64) -> Option<usize> {
870 assert!(sigma > 0.0);
871 assert!(target_accuracy > 0.0);
872 let m = mu.abs();
873 let s = sigma;
874 let gauss = (-(m * m) / (2.0 * s * s)).exp();
875 let coeff_mean = m * (2.0_f64 / std::f64::consts::PI).sqrt() * gauss / (s * s * s);
876 let coeff_deriv =
877 2.0 * gauss * (m * m - s * s).abs() / ((2.0 * std::f64::consts::PI).sqrt() * s.powi(5));
878 let asymptotic_index = |coeff: f64| -> f64 {
882 if !coeff.is_finite() || coeff <= target_accuracy {
883 0.0
884 } else {
885 (coeff / target_accuracy).sqrt() - 1.0
886 }
887 };
888 let peak_floor = m / (s * s) + 1.0;
892 let required = asymptotic_index(coeff_mean)
893 .max(asymptotic_index(coeff_deriv))
894 .max(peak_floor);
895 if !required.is_finite() || required > LOGIT_MAX_TERMS as f64 {
896 return None;
897 }
898 Some((required.ceil() as usize).max(4))
901}
902
903#[inline]
904fn stable_sigmoidwith_derivative(x: f64) -> (f64, f64) {
905 let x_clamped = x.clamp(-QUADRATURE_EXP_LOG_MAX, QUADRATURE_EXP_LOG_MAX);
906 if x_clamped != x {
907 return (sigmoid(x), 0.0);
908 }
909 if x_clamped >= 0.0 {
910 let z = (-x_clamped).exp();
911 let denom = 1.0 + z;
912 (1.0 / denom, z / (denom * denom))
913 } else {
914 let z = x_clamped.exp();
915 let denom = 1.0 + z;
916 (z / denom, z / (denom * denom))
917 }
918}
919
920#[inline]
921fn logit_tail_asymptotic(mu: f64, sigma: f64) -> Option<IntegratedMeanDerivative> {
922 if mu <= 0.0 {
927 let log_mean = mu + 0.5 * sigma * sigma;
928 if log_mean <= LOGIT_TAIL_LOG_MAX {
929 let mean = safe_exp(log_mean);
930 return Some(IntegratedMeanDerivative {
931 mean,
932 dmean_dmu: mean,
933 mode: IntegratedExpectationMode::ControlledAsymptotic,
934 });
935 }
936 } else {
937 let log_tail = -mu + 0.5 * sigma * sigma;
938 if log_tail <= LOGIT_TAIL_LOG_MAX {
939 let tail = safe_exp(log_tail);
940 return Some(IntegratedMeanDerivative {
941 mean: 1.0 - tail,
942 dmean_dmu: tail,
943 mode: IntegratedExpectationMode::ControlledAsymptotic,
944 });
945 }
946 }
947 None
948}
949
950#[inline]
951fn scaled_erfcx_termwith_derivative(m: f64, s: f64, x: f64, dxdm: f64) -> (f64, f64) {
952 let pref = 0.5 * (-(m * m) / (2.0 * s * s)).exp();
953 if x >= 0.0 {
954 let ex = erfcx_nonnegative(x);
955 let term = pref * ex;
956 let ex_prime = 2.0 * x * ex - std::f64::consts::FRAC_2_SQRT_PI;
957 let dterm = pref * ((-m / (s * s)) * ex + ex_prime * dxdm);
958 (term, dterm)
959 } else {
960 let lead = (x * x - (m * m) / (2.0 * s * s)).exp();
961 let dlead = lead * (2.0 * x * dxdm - m / (s * s));
962 let (rest, drest) = scaled_erfcx_termwith_derivative(m, s, -x, -dxdm);
963 (lead - rest, dlead - drest)
964 }
965}
966
967pub(crate) fn logit_posterior_meanwith_deriv_exact(
968 mu: f64,
969 sigma: f64,
970) -> Result<IntegratedMeanDerivative, EstimationError> {
971 if !(mu.is_finite() && sigma.is_finite()) {
990 crate::bail_invalid_estim!("logit exact expectation requires finite mu and sigma");
991 }
992 if sigma <= LOGIT_SIGMA_DEGENERATE {
993 let (mean, dmean_dmu) = stable_sigmoidwith_derivative(mu);
994 return Ok(IntegratedMeanDerivative {
995 mean,
996 dmean_dmu,
997 mode: IntegratedExpectationMode::ExactClosedForm,
998 });
999 }
1000 if let Some(out) = logit_tail_asymptotic(mu, sigma) {
1001 return Ok(out);
1002 }
1003 if logistic_normal_exact_eligible(mu, sigma)
1004 && let Ok(out) = logit_posterior_meanwith_deriv_exact_erfcx(mu, sigma)
1005 {
1006 return Ok(out);
1007 }
1008 Err(EstimationError::InvalidInput(
1017 "logit analytic expectation has no certified representation in this regime".to_string(),
1018 ))
1019}
1020
1021fn logit_posterior_meanwith_deriv_exact_erfcx(
1022 mu: f64,
1023 sigma: f64,
1024) -> Result<IntegratedMeanDerivative, EstimationError> {
1025 let m = mu.abs();
1051 let s = sigma;
1052 let z = SQRT_2 * s;
1053 let phi_term = gam_math::probability::normal_cdf(m / s);
1054 let phi_prime = gam_math::probability::normal_pdf(m / s) / s;
1055 let Some(max_k) = logistic_normal_series_cutoff(mu, sigma, LOGIT_ERFCX_ACCURACY_TARGET) else {
1056 crate::bail_invalid_estim!(
1057 "logit erfcx series truncation bound exceeds LOGIT_MAX_TERMS at the required accuracy"
1058 .to_string(),
1059 );
1060 };
1061
1062 let mut sum = 0.0_f64;
1063 let mut dsum = 0.0_f64;
1064 let mut k = 1usize;
1072 while k <= max_k {
1073 for kk in [k, k + 1].into_iter().filter(|kk| *kk <= max_k) {
1074 let kf = kk as f64;
1075 let a = (kf * s * s + m) / z;
1076 let b = (kf * s * s - m) / z;
1077 let sign = if kk % 2 == 1 { 1.0 } else { -1.0 };
1078 let (va, dva) = scaled_erfcx_termwith_derivative(m, s, a, 1.0 / z);
1079 let (vb, dvb) = scaled_erfcx_termwith_derivative(m, s, b, -1.0 / z);
1080 sum += sign * (va - vb);
1081 dsum += sign * (dva - dvb);
1082 }
1083 k += 2;
1084 }
1085
1086 let mut mean = phi_term + sum;
1087 let dmean = (phi_prime + dsum).max(0.0);
1088 if mu < 0.0 {
1089 mean = 1.0 - mean;
1090 }
1091 if !(mean.is_finite() && dmean.is_finite() && dmean >= 0.0) {
1092 crate::bail_invalid_estim!("logit erfcx expectation produced non-finite values");
1093 }
1094 Ok(IntegratedMeanDerivative {
1095 mean,
1096 dmean_dmu: dmean,
1097 mode: IntegratedExpectationMode::ExactSpecialFunction,
1098 })
1099}
1100
1101#[inline]
1106fn logit_posterior_meanwith_deriv_quadrature(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1107 let mean = integrate_normal_adaptive(mu, sigma, |x| stable_sigmoidwith_derivative(x).0);
1108 let dmean_dmu =
1109 integrate_normal_adaptive(mu, sigma, |x| stable_sigmoidwith_derivative(x).1).max(0.0);
1110 IntegratedMeanDerivative {
1111 mean,
1112 dmean_dmu,
1113 mode: IntegratedExpectationMode::QuadratureFallback,
1114 }
1115}
1116
1117#[inline]
1118fn logit_posterior_meanwith_deriv_controlled(
1119 mu: f64,
1120 sigma: f64,
1121) -> Result<IntegratedMeanDerivative, EstimationError> {
1122 if !(mu.is_finite() && sigma.is_finite()) {
1123 crate::bail_invalid_estim!("logit integrated moments require finite mu and sigma");
1124 }
1125 let candidate = match logit_posterior_meanwith_deriv_exact(mu, sigma) {
1126 Ok(out) => out,
1127 Err(_) => return Ok(logit_posterior_meanwith_deriv_quadrature(mu, sigma)),
1128 };
1129 match candidate.mode {
1141 IntegratedExpectationMode::ExactSpecialFunction
1142 | IntegratedExpectationMode::ControlledAsymptotic => {
1143 let reference = logit_posterior_meanwith_deriv_quadrature(mu, sigma);
1144 if integrated_mean_derivative_drift_exceeds(
1145 &candidate, &reference, 1e-6, 1e-4, 1e-7, 1e-3,
1146 ) {
1147 Ok(reference)
1148 } else {
1149 Ok(candidate)
1150 }
1151 }
1152 _ => Ok(candidate),
1153 }
1154}
1155
1156#[inline]
1157fn cloglog_extreme_asymptotic(mu: f64, sigma: f64) -> Option<IntegratedMeanDerivative> {
1158 let rare_log = mu + 0.5 * sigma * sigma;
1168 if rare_log <= CLOGLOG_RARE_EVENT_LOG_MAX {
1169 let mean = safe_exp(rare_log);
1170 return Some(IntegratedMeanDerivative {
1171 mean,
1172 dmean_dmu: mean,
1173 mode: IntegratedExpectationMode::ControlledAsymptotic,
1174 });
1175 }
1176 if mu - CLOGLOG_POSITIVE_SATURATION_SIGMAS * sigma >= CLOGLOG_POSITIVE_SATURATION_EDGE {
1177 return Some(IntegratedMeanDerivative {
1178 mean: 1.0,
1179 dmean_dmu: 0.0,
1180 mode: IntegratedExpectationMode::ControlledAsymptotic,
1181 });
1182 }
1183 None
1189}
1190
1191#[inline]
1192fn cloglog_survival_extreme_asymptotic(
1193 mu: f64,
1194 sigma: f64,
1195) -> Option<(f64, IntegratedExpectationMode)> {
1196 let rare_log = mu + 0.5 * sigma * sigma;
1197 if rare_log <= CLOGLOG_RARE_EVENT_LOG_MAX {
1198 let mean = safe_exp(rare_log);
1199 return Some((
1200 (1.0 - mean).clamp(0.0, 1.0),
1201 IntegratedExpectationMode::ControlledAsymptotic,
1202 ));
1203 }
1204 if mu - CLOGLOG_POSITIVE_SATURATION_SIGMAS * sigma >= CLOGLOG_POSITIVE_SATURATION_EDGE {
1205 return Some((0.0, IntegratedExpectationMode::ControlledAsymptotic));
1210 }
1211 None
1215}
1216
1217#[derive(Clone, Copy, Debug)]
1226struct LogSurvivalPanel {
1227 z_lo: f64,
1228 z_hi: f64,
1229 nodes: usize,
1230}
1231
1232#[derive(Clone, Copy, Debug, PartialEq, Eq)]
1240enum LogSurvivalBranch {
1241 Survival,
1243 Complement,
1245}
1246
1247impl LogSurvivalBranch {
1248 #[inline]
1251 fn log_tilt(self, mu: f64, sigma: f64, z: f64) -> f64 {
1252 let u = safe_exp(mu + sigma * z);
1253 match self {
1254 Self::Survival => -u,
1255 Self::Complement => {
1258 let m = (-u).exp_m1();
1259 if m == 0.0 {
1260 mu + sigma * z
1262 } else {
1263 (-m).ln()
1264 }
1265 }
1266 }
1267 }
1268
1269 #[inline]
1272 fn log_integrand(self, mu: f64, sigma: f64, z: f64) -> f64 {
1273 -0.5 * z * z + self.log_tilt(mu, sigma, z)
1274 }
1275
1276 #[inline]
1282 fn log_integrand_slope(self, mu: f64, sigma: f64, z: f64) -> f64 {
1283 let u = safe_exp(mu + sigma * z);
1284 match self {
1285 Self::Survival => -z - sigma * u,
1286 Self::Complement => {
1287 let em1 = u.exp_m1();
1288 let tilt_slope = if em1.is_finite() && em1 > 0.0 {
1289 sigma * u / em1
1290 } else if em1 == 0.0 {
1291 sigma
1293 } else {
1294 0.0
1295 };
1296 -z + tilt_slope
1297 }
1298 }
1299 }
1300}
1301
1302fn log_survival_peak_z(mu: f64, sigma: f64) -> f64 {
1312 let log_sigma = sigma.ln();
1313 let residual = |z: f64| log_sigma + mu + sigma * z - (-z).ln();
1314 let mut hi = -f64::MIN_POSITIVE;
1315 let mut lo = -1.0;
1316 let mut widen = 0;
1317 while residual(lo) > 0.0 && widen < 4096 {
1318 lo *= 2.0;
1319 widen += 1;
1320 if !lo.is_finite() {
1321 return f64::MIN;
1322 }
1323 }
1324 let mut z = if lo > -1.0e6 { 0.5 * (lo + hi) } else { lo };
1325 for _ in 0..200 {
1326 let r = residual(z);
1327 if r > 0.0 {
1328 hi = z;
1329 } else {
1330 lo = z;
1331 }
1332 let slope = sigma + 1.0 / (-z);
1334 let mut next = z - r / slope;
1335 if !(next > lo && next < hi) {
1336 next = 0.5 * (lo + hi);
1337 }
1338 if (next - z).abs() <= f64::EPSILON * (1.0 + z.abs()) {
1339 return next;
1340 }
1341 z = next;
1342 }
1343 z
1344}
1345
1346fn log_complement_peak_z(mu: f64, sigma: f64) -> f64 {
1354 let branch = LogSurvivalBranch::Complement;
1355 let (mut lo, mut hi) = (0.0_f64, sigma.max(f64::MIN_POSITIVE));
1356 if branch.log_integrand_slope(mu, sigma, hi) > 0.0 {
1357 return hi;
1358 }
1359 for _ in 0..200 {
1360 let mid = 0.5 * (lo + hi);
1361 if branch.log_integrand_slope(mu, sigma, mid) > 0.0 {
1362 lo = mid;
1363 } else {
1364 hi = mid;
1365 }
1366 if hi - lo <= f64::EPSILON * (1.0 + hi.abs()) {
1367 break;
1368 }
1369 }
1370 0.5 * (lo + hi)
1371}
1372
1373fn log_survival_panel_edge(
1380 branch: LogSurvivalBranch,
1381 mu: f64,
1382 sigma: f64,
1383 z_peak: f64,
1384 drop: f64,
1385 direction: f64,
1386) -> f64 {
1387 let peak_log = branch.log_integrand(mu, sigma, z_peak);
1388 let fallen = |z: f64| peak_log - branch.log_integrand(mu, sigma, z) >= drop;
1389 let mut step = 1.0_f64;
1390 let mut inner = z_peak;
1391 let mut outer = z_peak + direction * step;
1392 let mut widen = 0;
1393 while !fallen(outer) && widen < 4096 {
1394 inner = outer;
1395 step *= 2.0;
1396 outer = z_peak + direction * step;
1397 widen += 1;
1398 if !outer.is_finite() {
1399 return inner;
1400 }
1401 }
1402 let (mut lo, mut hi) = if direction > 0.0 {
1403 (inner, outer)
1404 } else {
1405 (outer, inner)
1406 };
1407 for _ in 0..200 {
1408 let mid = 0.5 * (lo + hi);
1409 if fallen(mid) == (direction > 0.0) {
1410 hi = mid;
1411 } else {
1412 lo = mid;
1413 }
1414 if hi - lo <= f64::EPSILON * (1.0 + mid.abs()) {
1415 break;
1416 }
1417 }
1418 if direction > 0.0 { hi } else { lo }
1419}
1420
1421fn log_survival_panel_arclength(mu: f64, sigma: f64, z_lo: f64, z_hi: f64) -> f64 {
1429 let scale_root = |z: f64| (1.0 + safe_exp(2.0 * sigma.ln() + mu + sigma * z)).sqrt();
1430 let (r_lo, r_hi) = (scale_root(z_lo), scale_root(z_hi));
1431 if !(r_lo.is_finite() && r_hi.is_finite()) {
1432 return f64::INFINITY;
1433 }
1434 (z_hi - z_lo) + 2.0 / sigma * ((r_hi - r_lo) - ((1.0 + r_hi) / (1.0 + r_lo)).ln())
1435}
1436
1437const LOG_SURVIVAL_TOWER_PANEL_ORDER: usize = LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER;
1460
1461#[inline]
1470fn log_survival_panel_placement_order(requested_order: usize) -> usize {
1471 if requested_order == 0 {
1472 0
1473 } else {
1474 LOG_SURVIVAL_TOWER_PANEL_ORDER
1475 }
1476}
1477
1478fn log_survival_panel(
1485 branch: LogSurvivalBranch,
1486 mu: f64,
1487 sigma: f64,
1488 requested_order: usize,
1489) -> LogSurvivalPanel {
1490 let order = log_survival_panel_placement_order(requested_order);
1491 let z_peak = match branch {
1492 LogSurvivalBranch::Survival => log_survival_peak_z(mu, sigma),
1493 LogSurvivalBranch::Complement => log_complement_peak_z(mu, sigma),
1494 };
1495 let drop = LOG_SURVIVAL_PANEL_LOG_DROP
1496 + LOG_SURVIVAL_PANEL_ORDER_LOG_DROP * (order as f64) * (2.0 + z_peak.abs()).ln();
1497 let z_lo = log_survival_panel_edge(branch, mu, sigma, z_peak, drop, -1.0);
1498 let z_hi = log_survival_panel_edge(branch, mu, sigma, z_peak, drop, 1.0);
1499 let arclength = log_survival_panel_arclength(mu, sigma, z_lo, z_hi);
1500 let requested = LOG_SURVIVAL_PANEL_ARCLENGTH_NODE_DENSITY * arclength
1501 + LOG_SURVIVAL_PANEL_SIGMA_NODE_SCALE * sigma.sqrt()
1502 + LOG_SURVIVAL_PANEL_ORDER_NODES * (order as f64);
1503 let mut nodes = if requested.is_finite() {
1504 (requested.ceil() as usize).clamp(
1505 LOG_SURVIVAL_PANEL_MIN_NODES,
1506 LOG_SURVIVAL_PANEL_MAX_NODES,
1507 )
1508 } else {
1509 LOG_SURVIVAL_PANEL_MAX_NODES
1510 };
1511 if nodes.is_multiple_of(2) {
1512 nodes += 1;
1513 }
1514 LogSurvivalPanel { z_lo, z_hi, nodes }
1515}
1516
1517#[derive(Clone, Copy, Debug)]
1525pub struct LogSurvivalSignedValue {
1526 pub log_abs: f64,
1527 pub sign: f64,
1528 pub log_cancellation: f64,
1529}
1530
1531impl LogSurvivalSignedValue {
1532 const ZERO: Self = Self {
1533 log_abs: f64::NEG_INFINITY,
1534 sign: 0.0,
1535 log_cancellation: f64::INFINITY,
1536 };
1537}
1538
1539#[derive(Clone, Debug)]
1541pub struct LogSurvivalJet {
1542 pub log_survival: f64,
1544 pub scaled_mu_derivatives: [LogSurvivalSignedValue; LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER + 1],
1546 pub order: usize,
1548 pub mode: IntegratedExpectationMode,
1549}
1550
1551impl LogSurvivalJet {
1552 pub fn certified_scaled_mu_derivatives(
1568 &self,
1569 order: usize,
1570 ) -> Option<&[LogSurvivalSignedValue]> {
1571 let prefix = self.certified_prefix_order(order)?;
1572 (prefix == order).then(|| &self.scaled_mu_derivatives[..=order])
1573 }
1574
1575 pub fn certified_prefix_order(&self, order: usize) -> Option<usize> {
1598 let order = order.min(self.order);
1599 let mut certified = None;
1600 for (index, entry) in self.scaled_mu_derivatives[..=order].iter().enumerate() {
1601 let usable = entry.log_abs.is_finite()
1602 && (index == 0 || entry.sign != 0.0)
1603 && entry.log_cancellation <= LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION;
1604 if !usable {
1605 break;
1606 }
1607 certified = Some(index);
1608 }
1609 certified
1610 }
1611}
1612
1613struct SignedLogAccumulator {
1615 running_max: f64,
1616 signed_sum: f64,
1617 abs_sum: f64,
1618}
1619
1620impl SignedLogAccumulator {
1621 #[inline]
1622 fn new() -> Self {
1623 Self {
1624 running_max: f64::NEG_INFINITY,
1625 signed_sum: 0.0,
1626 abs_sum: 0.0,
1627 }
1628 }
1629
1630 #[inline]
1631 fn push(&mut self, log_abs: f64, sign: f64) {
1632 if !log_abs.is_finite() {
1633 return;
1634 }
1635 if log_abs > self.running_max {
1636 let rescale = (self.running_max - log_abs).exp();
1637 self.signed_sum = self.signed_sum * rescale + sign;
1638 self.abs_sum = self.abs_sum * rescale + 1.0;
1639 self.running_max = log_abs;
1640 } else {
1641 let weight = (log_abs - self.running_max).exp();
1642 self.signed_sum += sign * weight;
1643 self.abs_sum += weight;
1644 }
1645 }
1646
1647 #[inline]
1648 fn finish(self) -> LogSurvivalSignedValue {
1649 if self.running_max == f64::NEG_INFINITY || self.signed_sum == 0.0 {
1650 return LogSurvivalSignedValue::ZERO;
1651 }
1652 LogSurvivalSignedValue {
1653 log_abs: self.running_max + self.signed_sum.abs().ln(),
1654 sign: self.signed_sum.signum(),
1655 log_cancellation: (self.abs_sum / self.signed_sum.abs()).ln().max(0.0),
1656 }
1657 }
1658}
1659
1660fn log_survival_panel_moments(
1666 ctx: &QuadratureContext,
1667 branch: LogSurvivalBranch,
1668 mu: f64,
1669 sigma: f64,
1670 order: usize,
1671) -> ([LogSurvivalSignedValue; LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER + 1], usize) {
1672 let panel = log_survival_panel(branch, mu, sigma, order);
1673 let rule = ctx.clenshaw_curtis_n(panel.nodes);
1674 let half = 0.5 * (panel.z_hi - panel.z_lo);
1675 let mid = 0.5 * (panel.z_hi + panel.z_lo);
1676 let log_gaussian_norm = -0.5 * (2.0 * std::f64::consts::PI).ln();
1677 let mut accumulators: Vec<SignedLogAccumulator> =
1678 (0..=order).map(|_| SignedLogAccumulator::new()).collect();
1679 for (&node, &weight) in rule.nodes.iter().zip(rule.weights.iter()) {
1680 let z = half * node + mid;
1681 let base = (weight * half).ln() + log_gaussian_norm - 0.5 * z * z
1682 + branch.log_tilt(mu, sigma, z);
1683 if !base.is_finite() {
1684 continue;
1685 }
1686 for (j, accumulator) in accumulators.iter_mut().enumerate() {
1687 let hermite = hermite_he(j, z);
1688 if hermite == 0.0 {
1689 continue;
1690 }
1691 accumulator.push(base + hermite.abs().ln(), hermite.signum());
1692 }
1693 }
1694 let mut out = [LogSurvivalSignedValue::ZERO; LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER + 1];
1695 for (slot, accumulator) in out.iter_mut().zip(accumulators.into_iter()) {
1696 *slot = accumulator.finish();
1697 }
1698 (out, panel.nodes)
1699}
1700
1701pub fn log_survival_jet(
1718 ctx: &QuadratureContext,
1719 mu: f64,
1720 sigma: f64,
1721 order: usize,
1722) -> LogSurvivalJet {
1723 let order = order.min(LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER);
1724 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
1725 let mut scaled = [LogSurvivalSignedValue::ZERO;
1728 LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER + 1];
1729 let log_survival = -safe_exp(mu);
1730 scaled[0] = LogSurvivalSignedValue {
1731 log_abs: log_survival,
1732 sign: 1.0,
1733 log_cancellation: 0.0,
1734 };
1735 return LogSurvivalJet {
1736 log_survival,
1737 scaled_mu_derivatives: scaled,
1738 order,
1739 mode: IntegratedExpectationMode::ExactClosedForm,
1740 };
1741 }
1742
1743 let (survival, _) =
1744 log_survival_panel_moments(ctx, LogSurvivalBranch::Survival, mu, sigma, order);
1745 let use_complement = survival[0].sign > 0.0 && survival[0].log_abs > -0.5;
1752 if !use_complement {
1753 let log_survival = if survival[0].sign > 0.0 {
1754 survival[0].log_abs
1755 } else {
1756 f64::NEG_INFINITY
1757 };
1758 return LogSurvivalJet {
1759 log_survival,
1760 scaled_mu_derivatives: survival,
1761 order,
1762 mode: IntegratedExpectationMode::ControlledAsymptotic,
1763 };
1764 }
1765
1766 let (complement, _) =
1767 log_survival_panel_moments(ctx, LogSurvivalBranch::Complement, mu, sigma, order);
1768 let mut scaled = survival;
1769 let log_survival = if complement[0].sign > 0.0 && complement[0].log_abs < 0.0 {
1773 (-safe_exp(complement[0].log_abs)).ln_1p()
1774 } else {
1775 survival[0].log_abs
1776 };
1777 scaled[0] = LogSurvivalSignedValue {
1778 log_abs: log_survival,
1779 sign: 1.0,
1780 log_cancellation: complement[0].log_cancellation,
1781 };
1782 for j in 1..=order {
1787 let flipped = LogSurvivalSignedValue {
1788 log_abs: complement[j].log_abs,
1789 sign: -complement[j].sign,
1790 log_cancellation: complement[j].log_cancellation,
1791 };
1792 if flipped.sign != 0.0 && flipped.log_cancellation < scaled[j].log_cancellation {
1793 scaled[j] = flipped;
1794 }
1795 }
1796 LogSurvivalJet {
1797 log_survival,
1798 scaled_mu_derivatives: scaled,
1799 order,
1800 mode: IntegratedExpectationMode::ControlledAsymptotic,
1801 }
1802}
1803
1804#[inline]
1805fn hermite_he(n: usize, x: f64) -> f64 {
1806 let mut previous = 1.0_f64;
1807 if n == 0 {
1808 return previous;
1809 }
1810 let mut current = x;
1811 for order in 1..n {
1812 let next = x * current - (order as f64) * previous;
1813 previous = current;
1814 current = next;
1815 }
1816 current
1817}
1818
1819pub fn log_survival_scaled_mu_derivative(
1830 ctx: &QuadratureContext,
1831 mu: f64,
1832 sigma: f64,
1833 order: usize,
1834) -> (f64, f64) {
1835 assert!(
1836 order >= 1,
1837 "the log-space survival mu-derivative tower starts at order 1; \
1838 order 0 is ln S, which log_survival_jet returns"
1839 );
1840 let jet = log_survival_jet(ctx, mu, sigma, order);
1841 let entry = jet.scaled_mu_derivatives[order.min(jet.order)];
1842 (entry.log_abs, entry.sign)
1843}
1844
1845pub(crate) fn cloglog_log_survival_term_controlled(
1863 ctx: &QuadratureContext,
1864 mu: f64,
1865 sigma: f64,
1866) -> (f64, IntegratedExpectationMode) {
1867 let jet = log_survival_jet(ctx, mu, sigma, 0);
1868 (jet.log_survival, jet.mode)
1869}
1870
1871#[inline]
1890fn gumbel_survival(x: f64) -> f64 {
1891 (-safe_exp(x)).exp()
1892}
1893
1894#[inline]
1899fn cloglog_mean_d1_exact(x: f64) -> f64 {
1900 let ex = safe_exp(x);
1901 if ex.is_infinite() {
1902 0.0
1903 } else {
1904 ex * (-ex).exp()
1905 }
1906}
1907
1908#[inline]
1918fn cloglog_mean_exact(x: f64) -> f64 {
1919 cloglog_negative_tail_mean(x)
1920}
1921
1922#[inline]
1938fn cloglog_negative_tail_mean(eta: f64) -> f64 {
1939 if eta < -745.0 {
1943 0.0
1945 } else {
1946 let ex = safe_exp(eta);
1949 -(-ex).exp_m1()
1950 }
1951}
1952
1953#[inline]
1958fn cloglog_small_sigma_taylor(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
1959 if sigma <= CLOGLOG_SIGMA_DEGENERATE {
1983 return IntegratedMeanDerivative {
1984 mean: cloglog_mean_exact(mu),
1985 dmean_dmu: cloglog_mean_d1_exact(mu),
1986 mode: IntegratedExpectationMode::ExactClosedForm,
1987 };
1988 }
1989
1990 let ex = safe_exp(mu);
1991 if !ex.is_finite() {
1992 return IntegratedMeanDerivative {
1994 mean: 1.0,
1995 dmean_dmu: 0.0,
1996 mode: IntegratedExpectationMode::ControlledAsymptotic,
1997 };
1998 }
1999 let surv = (-ex).exp();
2000 if surv == 0.0 {
2001 return IntegratedMeanDerivative {
2003 mean: 1.0,
2004 dmean_dmu: 0.0,
2005 mode: IntegratedExpectationMode::ControlledAsymptotic,
2006 };
2007 }
2008
2009 let s2 = sigma * sigma;
2010 let s4 = s2 * s2;
2011 let s6 = s4 * s2;
2012 let s8 = s4 * s4;
2013 let e2x = ex * ex;
2014 let e3x = e2x * ex;
2015 let e4x = e3x * ex;
2016 let e5x = e4x * ex;
2017 let e6x = e5x * ex;
2018 let e7x = e6x * ex;
2019 let e8x = e7x * ex;
2020 let e9x = e8x * ex;
2021 let f0 = -(-ex).exp_m1();
2035 let f1 = ex * surv;
2036 let f2 = surv * (ex - e2x);
2037 let f3 = surv * (ex - 3.0 * e2x + e3x);
2038 let f4 = surv * (ex - 7.0 * e2x + 6.0 * e3x - e4x);
2039 let f5 = surv * (ex - 15.0 * e2x + 25.0 * e3x - 10.0 * e4x + e5x);
2040 let f6 = surv * (ex - 31.0 * e2x + 90.0 * e3x - 65.0 * e4x + 15.0 * e5x - e6x);
2041 let f7 = surv * (ex - 63.0 * e2x + 301.0 * e3x - 350.0 * e4x + 140.0 * e5x - 21.0 * e6x + e7x);
2042 let f8 = surv
2043 * (ex - 127.0 * e2x + 966.0 * e3x - 1701.0 * e4x + 1050.0 * e5x - 266.0 * e6x + 28.0 * e7x
2044 - e8x);
2045 let f9 = surv
2046 * (ex - 255.0 * e2x + 3025.0 * e3x - 7770.0 * e4x + 6951.0 * e5x - 2646.0 * e6x
2047 + 462.0 * e7x
2048 - 36.0 * e8x
2049 + e9x);
2050 IntegratedMeanDerivative {
2054 mean: f0 + 0.5 * s2 * f2 + (s4 / 8.0) * f4 + (s6 / 48.0) * f6 + (s8 / 384.0) * f8,
2055 dmean_dmu: (f1 + 0.5 * s2 * f3 + (s4 / 8.0) * f5 + (s6 / 48.0) * f7 + (s8 / 384.0) * f9)
2056 .max(0.0),
2057 mode: IntegratedExpectationMode::ControlledAsymptotic,
2058 }
2059}
2060
2061#[inline]
2062fn adaptive_simpson_refine(
2067 g: &impl Fn(f64) -> f64,
2068 a: f64,
2069 b: f64,
2070 fa: f64,
2071 fb: f64,
2072 fm: f64,
2073 whole: f64,
2074 tol: f64,
2075 depth: i32,
2076) -> f64 {
2077 let m = 0.5 * (a + b);
2078 let lm = 0.5 * (a + m);
2079 let rm = 0.5 * (m + b);
2080 let flm = g(lm);
2081 let frm = g(rm);
2082 let left = (m - a) / 6.0 * (fa + 4.0 * flm + fm);
2083 let right = (b - m) / 6.0 * (fm + 4.0 * frm + fb);
2084 let est = left + right;
2085 if depth <= 0 || (est - whole).abs() <= 15.0 * tol {
2086 return est + (est - whole) / 15.0;
2087 }
2088 adaptive_simpson_refine(g, a, m, fa, fm, flm, left, 0.5 * tol, depth - 1)
2089 + adaptive_simpson_refine(g, m, b, fm, fb, frm, right, 0.5 * tol, depth - 1)
2090}
2091
2092fn integrate_normal_adaptive(mu: f64, sigma: f64, f: impl Fn(f64) -> f64) -> f64 {
2107 if !(sigma.is_finite()) || sigma < 1e-10 {
2108 return f(mu);
2109 }
2110 const K: f64 = 15.0;
2111 const INITIAL_PANELS: usize = 24;
2112 const TOL: f64 = 1e-12;
2113 const MAX_DEPTH: i32 = 40;
2114 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
2115 let g = |u: f64| f(mu + sigma * u) * inv_sqrt_2pi * (-0.5 * u * u).exp();
2119 let panel = 2.0 * K / INITIAL_PANELS as f64;
2120 let mut total = 0.0;
2121 for p in 0..INITIAL_PANELS {
2122 let a = -K + p as f64 * panel;
2123 let b = a + panel;
2124 let fa = g(a);
2125 let fb = g(b);
2126 let fm = g(0.5 * (a + b));
2127 let whole = (b - a) / 6.0 * (fa + 4.0 * fm + fb);
2128 total += adaptive_simpson_refine(&g, a, b, fa, fb, fm, whole, TOL, MAX_DEPTH);
2129 }
2130 total
2131}
2132
2133fn cloglog_posterior_meanwith_deriv_quadrature(mu: f64, sigma: f64) -> IntegratedMeanDerivative {
2134 if sigma < 1e-10 {
2135 return IntegratedMeanDerivative {
2136 mean: cloglog_mean_exact(mu),
2137 dmean_dmu: cloglog_mean_d1_exact(mu),
2138 mode: IntegratedExpectationMode::ExactClosedForm,
2139 };
2140 }
2141 let mean = cloglog_mean_from_survival(survival_posterior_mean_quadrature(mu, sigma));
2142 let dmean_dmu = integrate_normal_adaptive(mu, sigma, cloglog_mean_d1_exact).max(0.0);
2143 IntegratedMeanDerivative {
2144 mean,
2145 dmean_dmu,
2146 mode: IntegratedExpectationMode::QuadratureFallback,
2147 }
2148}
2149
2150#[inline]
2151fn survival_posterior_mean_quadrature(eta: f64, se_eta: f64) -> f64 {
2152 integrate_normal_adaptive(eta, se_eta, gumbel_survival).clamp(0.0, 1.0)
2153}
2154
2155fn cloglog_survival_term_controlled(
2156 ctx: &QuadratureContext,
2157 mu: f64,
2158 sigma: f64,
2159) -> (f64, IntegratedExpectationMode) {
2160 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
2197 return (
2198 gumbel_survival(mu).clamp(0.0, 1.0),
2199 IntegratedExpectationMode::ExactClosedForm,
2200 );
2201 }
2202 if sigma < CLOGLOG_SIGMA_TAYLOR_MAX {
2203 let mean = cloglog_small_sigma_taylor(mu, sigma).mean;
2204 return (
2205 (1.0 - mean).clamp(0.0, 1.0),
2206 IntegratedExpectationMode::ControlledAsymptotic,
2207 );
2208 }
2209 if let Some(out) = cloglog_survival_extreme_asymptotic(mu, sigma) {
2210 return out;
2211 }
2212 if sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN {
2213 let log_s = log_survival_jet(ctx, mu, sigma, 0).log_survival;
2220 return (
2221 safe_exp(log_s).clamp(0.0, 1.0),
2222 IntegratedExpectationMode::ControlledAsymptotic,
2223 );
2224 }
2225 if cloglog_survival_miles_is_reliable(mu, sigma)
2226 && let Ok(out) = cloglog_survival_miles(mu, sigma)
2227 {
2228 return (
2229 out.clamp(0.0, 1.0),
2230 IntegratedExpectationMode::ExactSpecialFunction,
2231 );
2232 }
2233 if cloglog_should_prefer_cc(mu, sigma, CLOGLOG_CC_TOL)
2234 && let Ok(out) = cloglog_survival_cc(ctx, mu, sigma, CLOGLOG_CC_TOL)
2235 {
2236 return (
2237 out.clamp(0.0, 1.0),
2238 IntegratedExpectationMode::ExactSpecialFunction,
2239 );
2240 }
2241 if let Ok(out) = cloglog_survival_gamma_reference(mu, sigma) {
2242 return (
2243 out.clamp(0.0, 1.0),
2244 IntegratedExpectationMode::ExactSpecialFunction,
2245 );
2246 }
2247 (
2248 survival_posterior_mean_quadrature(mu, sigma),
2249 IntegratedExpectationMode::QuadratureFallback,
2250 )
2251}
2252
2253#[inline]
2254fn lognormal_laplace_term_controlled(
2255 ctx: &QuadratureContext,
2256 z: f64,
2257 mu: f64,
2258 sigma: f64,
2259) -> (f64, IntegratedExpectationMode) {
2260 if !(z.is_finite() && z > 0.0) {
2286 return (f64::NAN, IntegratedExpectationMode::QuadratureFallback);
2287 }
2288 lognormal_laplace_unit_term_shared(ctx, mu + z.ln(), sigma)
2289}
2290
2291#[inline]
2292pub(crate) fn lognormal_laplace_unit_term_shared(
2293 ctx: &QuadratureContext,
2294 shifted_mu: f64,
2295 sigma: f64,
2296) -> (f64, IntegratedExpectationMode) {
2297 cloglog_survival_term_controlled(ctx, shifted_mu, sigma)
2298}
2299
2300#[inline]
2304pub fn lognormal_laplace_unit_log_term_shared(
2305 ctx: &QuadratureContext,
2306 shifted_mu: f64,
2307 sigma: f64,
2308) -> (f64, IntegratedExpectationMode) {
2309 cloglog_log_survival_term_controlled(ctx, shifted_mu, sigma)
2310}
2311
2312#[inline]
2313fn cloglog_survivalsecond_moment_controlled(
2314 ctx: &QuadratureContext,
2315 mu: f64,
2316 sigma: f64,
2317) -> (f64, IntegratedExpectationMode) {
2318 lognormal_laplace_term_controlled(ctx, 2.0, mu, sigma)
2333}
2334
2335#[inline]
2336fn cloglog_survival_pair_controlled(
2337 ctx: &QuadratureContext,
2338 mu: f64,
2339 sigma: f64,
2340) -> (
2341 (f64, IntegratedExpectationMode),
2342 (f64, IntegratedExpectationMode),
2343) {
2344 let shiftedmu = mu + sigma * sigma;
2345
2346 if cloglog_survival_miles_is_reliable(mu, sigma)
2357 && cloglog_survival_miles_is_reliable(shiftedmu, sigma)
2358 && let (Ok(base), Ok(shifted)) = (
2359 cloglog_survival_miles(mu, sigma),
2360 cloglog_survival_miles(shiftedmu, sigma),
2361 )
2362 {
2363 return (
2364 (
2365 base.clamp(0.0, 1.0),
2366 IntegratedExpectationMode::ExactSpecialFunction,
2367 ),
2368 (
2369 shifted.clamp(0.0, 1.0),
2370 IntegratedExpectationMode::ExactSpecialFunction,
2371 ),
2372 );
2373 }
2374
2375 if cloglog_should_prefer_cc(mu, sigma, CLOGLOG_CC_TOL)
2376 && cloglog_should_prefer_cc(shiftedmu, sigma, CLOGLOG_CC_TOL)
2377 && let (Ok(base), Ok(shifted)) = (
2378 cloglog_survival_cc(ctx, mu, sigma, CLOGLOG_CC_TOL),
2379 cloglog_survival_cc(ctx, shiftedmu, sigma, CLOGLOG_CC_TOL),
2380 )
2381 {
2382 return (
2383 (
2384 base.clamp(0.0, 1.0),
2385 IntegratedExpectationMode::ExactSpecialFunction,
2386 ),
2387 (
2388 shifted.clamp(0.0, 1.0),
2389 IntegratedExpectationMode::ExactSpecialFunction,
2390 ),
2391 );
2392 }
2393
2394 if let (Ok(base), Ok(shifted)) = (
2395 cloglog_survival_gamma_reference(mu, sigma),
2396 cloglog_survival_gamma_reference(shiftedmu, sigma),
2397 ) {
2398 return (
2399 (
2400 base.clamp(0.0, 1.0),
2401 IntegratedExpectationMode::ExactSpecialFunction,
2402 ),
2403 (
2404 shifted.clamp(0.0, 1.0),
2405 IntegratedExpectationMode::ExactSpecialFunction,
2406 ),
2407 );
2408 }
2409
2410 (
2411 cloglog_survival_term_controlled(ctx, mu, sigma),
2412 cloglog_survival_term_controlled(ctx, shiftedmu, sigma),
2413 )
2414}
2415
2416#[inline]
2417fn cloglog_mean_from_survival(survival: f64) -> f64 {
2418 let survival = survival.clamp(0.0, 1.0);
2419 if survival > 0.5 {
2420 -survival.ln().exp_m1()
2430 } else {
2431 1.0 - survival
2432 }
2433}
2434
2435#[inline]
2436fn cloglog_shift_identity_derivative(mu: f64, sigma: f64, shifted_survival: f64) -> f64 {
2437 if !(mu.is_finite() && sigma.is_finite()) || shifted_survival <= 0.0 {
2451 return 0.0;
2452 }
2453 cloglog_shift_identity_derivative_log(mu, sigma, shifted_survival.ln())
2454}
2455
2456#[inline]
2466fn cloglog_shift_identity_derivative_log(mu: f64, sigma: f64, log_shifted_survival: f64) -> f64 {
2467 if !(mu.is_finite() && sigma.is_finite()) || log_shifted_survival == f64::NEG_INFINITY {
2468 return 0.0;
2469 }
2470 let log_derivative = mu + 0.5 * sigma * sigma + log_shifted_survival;
2471 let upper = 1.0 / std::f64::consts::E;
2472 if !log_derivative.is_finite() {
2473 return upper;
2476 }
2477 safe_exp(log_derivative).clamp(0.0, upper)
2478}
2479
2480#[inline]
2481fn log_half_erfc_stable(u: f64) -> f64 {
2482 if u > 0.0 {
2494 -u * u + (0.5 * erfcx_nonnegative(u)).ln()
2495 } else {
2496 normal_logcdf(-u * SQRT_2)
2497 }
2498}
2499
2500#[inline]
2540fn cloglog_survival_miles_is_reliable(mu: f64, sigma: f64) -> bool {
2541 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0) {
2542 return false;
2543 }
2544 let alpha_ln = CLOGLOG_MILES_ALPHA.ln();
2545 let shifted = mu - alpha_ln;
2546 let peak_log = CLOGLOG_MILES_ALPHA - 0.5 * shifted * shifted / (sigma * sigma);
2547 peak_log.is_finite() && peak_log <= CLOGLOG_MILES_PEAK_LOG_MAX
2548}
2549
2550fn cloglog_survival_miles(mu: f64, sigma: f64) -> Result<f64, EstimationError> {
2551 let alpha_ln = CLOGLOG_MILES_ALPHA.ln();
2581 let mut s_sum = 0.0_f64;
2582 let mut stable_pairs = 0usize;
2583
2584 for pair_start in (0..CLOGLOG_MILES_MAX_TERMS).step_by(2) {
2585 let mut pair_s = 0.0_f64;
2586 for n in pair_start..(pair_start + 2).min(CLOGLOG_MILES_MAX_TERMS) {
2587 let nf = n as f64;
2588 let sign = if n % 2 == 0 { 1.0 } else { -1.0 };
2589 let base_log = nf * mu + 0.5 * sigma * sigma * nf * nf
2590 - statrs::function::gamma::ln_gamma(nf + 1.0);
2591 let u = (mu - alpha_ln + sigma * sigma * nf) / (SQRT_2 * sigma);
2592 let log_half_erfc = log_half_erfc_stable(u);
2593 let term_log = base_log + log_half_erfc;
2594 if term_log > QUADRATURE_EXP_LOG_MAX {
2595 crate::bail_invalid_estim!("Miles cloglog series term exceeded finite exp range");
2596 }
2597 let term = sign * safe_exp(term_log);
2598 pair_s += term;
2599 }
2600 s_sum += pair_s;
2601
2602 let s_scale = s_sum.abs().max(1.0);
2603 if pair_s.abs() <= 2e-15 * s_scale {
2604 stable_pairs += 1;
2605 if stable_pairs >= SERIES_CONSECUTIVE_SMALL_TERMS {
2606 if s_sum.is_finite() && (-1e-10..=1.0 + 1e-10).contains(&s_sum) {
2607 return Ok(s_sum.clamp(0.0, 1.0));
2608 }
2609 break;
2610 }
2611 } else {
2612 stable_pairs = 0;
2613 }
2614 }
2615
2616 Err(EstimationError::InvalidInput(
2617 "Miles cloglog series did not converge safely".to_string(),
2618 ))
2619}
2620
2621fn cloglog_survival_cc(
2622 ctx: &QuadratureContext,
2623 mu: f64,
2624 sigma: f64,
2625 tol: f64,
2626) -> Result<f64, EstimationError> {
2627 if !(mu.is_finite() && sigma.is_finite() && sigma > 0.0 && tol.is_finite() && tol > 0.0) {
2628 crate::bail_invalid_estim!(
2629 "CC cloglog backend requires finite mu, positive sigma, and positive tolerance"
2630 .to_string(),
2631 );
2632 }
2633
2634 let p_tail = (tol / 8.0).clamp(1e-300, 0.25);
2672 let a = gam_math::probability::standard_normal_quantile(p_tail)
2673 .map(|z| -z)
2674 .unwrap_or(8.0)
2675 .max(1.0);
2676 let n = cloglog_cc_required_nodes(mu, sigma, tol)?;
2677 if n > CLOGLOG_CC_NODE_CAP {
2678 crate::bail_invalid_estim!("CC cloglog backend requires too many nodes");
2679 }
2680
2681 let rule = ctx.clenshaw_curtis_n(n);
2682 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
2683 let mut sum = 0.0_f64;
2684 let mut c = 0.0_f64;
2685 for (&x, &w) in rule.nodes.iter().zip(rule.weights.iter()) {
2686 let t = a * x;
2687 let u = mu + sigma * t;
2688 let e = safe_exp(u);
2689 let w0 = (-0.5 * t * t).exp() * inv_sqrt_2pi;
2690 let yk = w * w0 * (-e).exp() - c;
2691 let tk = sum + yk;
2692 c = (tk - sum) - yk;
2693 sum = tk;
2694 }
2695
2696 let survival = (a * sum).clamp(0.0, 1.0);
2697 if !survival.is_finite() {
2698 crate::bail_invalid_estim!("CC cloglog backend produced non-finite values");
2699 }
2700 Ok(survival)
2701}
2702
2703#[inline]
2704fn complex_add(a: Complex, b: Complex) -> Complex {
2705 Complex {
2706 re: a.re + b.re,
2707 im: a.im + b.im,
2708 }
2709}
2710
2711#[inline]
2712fn complex_sub(a: Complex, b: Complex) -> Complex {
2713 Complex {
2714 re: a.re - b.re,
2715 im: a.im - b.im,
2716 }
2717}
2718
2719#[inline]
2720fn complexmul(a: Complex, b: Complex) -> Complex {
2721 Complex {
2722 re: a.re * b.re - a.im * b.im,
2723 im: a.re * b.im + a.im * b.re,
2724 }
2725}
2726
2727#[inline]
2728fn complex_div(a: Complex, b: Complex) -> Complex {
2729 let den = (b.re * b.re + b.im * b.im).max(1e-300);
2730 Complex {
2731 re: (a.re * b.re + a.im * b.im) / den,
2732 im: (a.im * b.re - a.re * b.im) / den,
2733 }
2734}
2735
2736#[inline]
2737fn complex_abs(z: Complex) -> f64 {
2738 z.re.hypot(z.im)
2739}
2740
2741#[inline]
2742fn complex_ln(z: Complex) -> Complex {
2743 Complex {
2744 re: complex_abs(z).ln(),
2745 im: z.im.atan2(z.re),
2746 }
2747}
2748
2749#[inline]
2750fn complex_exp(z: Complex) -> Complex {
2751 let e = z.re.exp();
2752 Complex {
2753 re: e * z.im.cos(),
2754 im: e * z.im.sin(),
2755 }
2756}
2757
2758#[inline]
2759fn complex_sin(z: Complex) -> Complex {
2760 Complex {
2761 re: z.re.sin() * z.im.cosh(),
2762 im: z.re.cos() * z.im.sinh(),
2763 }
2764}
2765
2766fn complex_log_gamma_lanczos(z: Complex) -> Complex {
2767 const G: f64 = 7.0;
2771 const COEFFS: [f64; 9] = [
2772 0.999_999_999_999_809_9,
2773 676.520_368_121_885_1,
2774 -1_259.139_216_722_402_8,
2775 771.323_428_777_653_1,
2776 -176.615_029_162_140_6,
2777 12.507_343_278_686_905,
2778 -0.138_571_095_265_720_12,
2779 9.984_369_578_019_572e-6,
2780 1.505_632_735_149_311_6e-7,
2781 ];
2782
2783 if z.re < 0.5 {
2784 let piz = Complex {
2785 re: std::f64::consts::PI * z.re,
2786 im: std::f64::consts::PI * z.im,
2787 };
2788 let one_minusz = Complex {
2789 re: 1.0 - z.re,
2790 im: -z.im,
2791 };
2792 return complex_sub(
2793 complex_sub(
2794 Complex {
2795 re: std::f64::consts::PI.ln(),
2796 im: 0.0,
2797 },
2798 complex_ln(complex_sin(piz)),
2799 ),
2800 complex_log_gamma_lanczos(one_minusz),
2801 );
2802 }
2803
2804 let z1 = Complex {
2805 re: z.re - 1.0,
2806 im: z.im,
2807 };
2808 let mut x = Complex {
2809 re: COEFFS[0],
2810 im: 0.0,
2811 };
2812 for (i, c) in COEFFS.iter().enumerate().skip(1) {
2813 x = complex_add(
2814 x,
2815 complex_div(
2816 Complex { re: *c, im: 0.0 },
2817 Complex {
2818 re: z1.re + i as f64,
2819 im: z1.im,
2820 },
2821 ),
2822 );
2823 }
2824 let t = Complex {
2825 re: z1.re + G + 0.5,
2826 im: z1.im,
2827 };
2828 complex_add(
2829 complex_add(
2830 Complex {
2831 re: 0.5 * (2.0 * std::f64::consts::PI).ln(),
2832 im: 0.0,
2833 },
2834 complexmul(
2835 Complex {
2836 re: z1.re + 0.5,
2837 im: z1.im,
2838 },
2839 complex_ln(t),
2840 ),
2841 ),
2842 complex_sub(complex_ln(x), t),
2843 )
2844}
2845
2846fn cloglog_survival_gamma_reference(mu: f64, sigma: f64) -> Result<f64, EstimationError> {
2850 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 0.0 {
2851 crate::bail_invalid_estim!(
2852 "Gamma cloglog reference backend requires finite mu and positive sigma"
2853 );
2854 }
2855
2856 let n = (CLOGLOG_GAMMA_T_MAX_REF / CLOGLOG_GAMMA_H_REF).round() as usize;
2890 let n = if n.is_multiple_of(2) { n } else { n + 1 };
2891 let h = CLOGLOG_GAMMA_T_MAX_REF / n as f64;
2892
2893 let eval = |t: f64| -> f64 {
2894 let z = Complex {
2895 re: CLOGLOG_GAMMA_K_REF,
2896 im: t,
2897 };
2898 let log_gamma = complex_log_gamma_lanczos(z);
2899 let z_sq = complexmul(z, z);
2900 let exponent = complex_sub(
2901 complex_add(
2902 log_gamma,
2903 Complex {
2904 re: 0.5 * sigma * sigma * z_sq.re,
2905 im: 0.5 * sigma * sigma * z_sq.im,
2906 },
2907 ),
2908 Complex {
2909 re: mu * z.re,
2910 im: mu * z.im,
2911 },
2912 );
2913 complex_exp(exponent).re
2914 };
2915
2916 let f0 = eval(0.0);
2917 let fn_ = eval(CLOGLOG_GAMMA_T_MAX_REF);
2918 let mut sum_s = f0 + fn_;
2919 for i in 1..n {
2920 let t = i as f64 * h;
2921 let fi = eval(t);
2922 let w = if i % 2 == 0 { 2.0 } else { 4.0 };
2923 sum_s += w * fi;
2924 }
2925 let sval = ((h / 3.0) * sum_s / std::f64::consts::PI).clamp(0.0, 1.0);
2926 if !sval.is_finite() {
2927 crate::bail_invalid_estim!("Gamma cloglog reference backend produced non-finite values");
2928 }
2929 Ok(sval)
2930}
2931
2932pub(crate) fn cloglog_posterior_meanwith_deriv_controlled(
2933 ctx: &QuadratureContext,
2934 mu: f64,
2935 sigma: f64,
2936) -> IntegratedMeanDerivative {
2937 if !(mu.is_finite() && sigma.is_finite()) || sigma <= CLOGLOG_SIGMA_DEGENERATE {
2982 return IntegratedMeanDerivative {
2983 mean: cloglog_mean_exact(mu),
2986 dmean_dmu: cloglog_mean_d1_exact(mu),
2989 mode: IntegratedExpectationMode::ExactClosedForm,
2990 };
2991 }
2992 if sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN {
2993 let (log_base, base_mode) = cloglog_log_survival_term_controlled(ctx, mu, sigma);
2999 let (log_shift, shift_mode) =
3000 cloglog_log_survival_term_controlled(ctx, mu + sigma * sigma, sigma);
3001 let mean = (-log_base.exp_m1()).clamp(0.0, 1.0);
3003 let dmean = cloglog_shift_identity_derivative_log(mu, sigma, log_shift);
3004 return IntegratedMeanDerivative {
3005 mean,
3006 dmean_dmu: dmean.max(0.0),
3007 mode: worse_integrated_expectation_mode(base_mode, shift_mode),
3008 };
3009 }
3010 let candidate = if sigma < CLOGLOG_SIGMA_TAYLOR_MAX {
3011 cloglog_small_sigma_taylor(mu, sigma)
3012 } else if let Some(out) = cloglog_extreme_asymptotic(mu, sigma) {
3013 out
3014 } else {
3015 let ((survival, mode), (shifted_survival, shifted_mode)) =
3016 cloglog_survival_pair_controlled(ctx, mu, sigma);
3017 if matches!(mode, IntegratedExpectationMode::QuadratureFallback)
3018 || matches!(shifted_mode, IntegratedExpectationMode::QuadratureFallback)
3019 {
3020 return cloglog_posterior_meanwith_deriv_quadrature(mu, sigma);
3021 }
3022 let mean = cloglog_mean_from_survival(survival);
3023 let dmean = cloglog_shift_identity_derivative(mu, sigma, shifted_survival);
3024 let mode = if matches!(mode, IntegratedExpectationMode::ControlledAsymptotic)
3025 || matches!(
3026 shifted_mode,
3027 IntegratedExpectationMode::ControlledAsymptotic
3028 ) {
3029 IntegratedExpectationMode::ControlledAsymptotic
3030 } else {
3031 mode
3032 };
3033 IntegratedMeanDerivative {
3034 mean,
3035 dmean_dmu: dmean.max(0.0),
3036 mode,
3037 }
3038 };
3039 if matches!(
3045 candidate.mode,
3046 IntegratedExpectationMode::ControlledAsymptotic
3047 ) && sigma >= CLOGLOG_LARGE_SIGMA_ASYMPTOTIC_MIN
3048 {
3049 return candidate;
3050 }
3051 let ghq = cloglog_posterior_meanwith_deriv_quadrature(mu, sigma);
3052 if integrated_mean_derivative_drift_exceeds(&candidate, &ghq, 1e-6, 1e-4, 1e-7, 1e-3) {
3059 ghq
3060 } else {
3061 candidate
3062 }
3063}
3064
3065pub fn integrated_inverse_link_mean_and_derivative(
3066 quadctx: &QuadratureContext,
3067 link: LinkFunction,
3068 mu: f64,
3069 sigma: f64,
3070) -> Result<IntegratedMeanDerivative, EstimationError> {
3071 match link {
3100 LinkFunction::Log => {
3101 let (mean, saturated) = safe_expwith_saturation(mu + 0.5 * sigma * sigma);
3102 Ok(IntegratedMeanDerivative {
3103 mean,
3104 dmean_dmu: mean,
3105 mode: if saturated {
3106 IntegratedExpectationMode::ControlledAsymptotic
3107 } else {
3108 IntegratedExpectationMode::ExactClosedForm
3109 },
3110 })
3111 }
3112 LinkFunction::Probit => Ok(probit_posterior_meanwith_deriv_exact(mu, sigma)),
3113 LinkFunction::Logit => logit_posterior_meanwith_deriv_controlled(mu, sigma),
3114 LinkFunction::CLogLog => Ok(cloglog_posterior_meanwith_deriv_controlled(quadctx, mu, sigma)),
3115 LinkFunction::LogLog | LinkFunction::Cauchit => {
3116 let component = if matches!(link, LinkFunction::LogLog) {
3118 LinkComponent::LogLog
3119 } else {
3120 LinkComponent::Cauchit
3121 };
3122 let (mean, dmean_dmu, _, _) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
3123 component_point_jet(component, x)
3124 });
3125 Ok(IntegratedMeanDerivative {
3126 mean,
3127 dmean_dmu,
3128 mode: if sigma <= 1e-10 {
3129 IntegratedExpectationMode::ExactClosedForm
3130 } else {
3131 IntegratedExpectationMode::QuadratureFallback
3132 },
3133 })
3134 }
3135 LinkFunction::Sas => Err(EstimationError::InvalidInput(
3136 "state-less integrated SAS moments are unsupported; use SAS-aware prediction APIs with explicit (epsilon, log_delta)".to_string(),
3137 )),
3138 LinkFunction::BetaLogistic => Err(EstimationError::InvalidInput(
3139 "state-less integrated Beta-Logistic moments are unsupported; use link-aware prediction APIs with explicit (delta, epsilon)".to_string(),
3140 )),
3141 LinkFunction::Identity => Ok(IntegratedMeanDerivative {
3142 mean: mu,
3143 dmean_dmu: 1.0,
3144 mode: IntegratedExpectationMode::ExactClosedForm,
3145 }),
3146 }
3147}
3148
3149#[inline]
3150pub fn integrated_inverse_link_jet(
3151 quadctx: &QuadratureContext,
3152 link: LinkFunction,
3153 mu: f64,
3154 sigma: f64,
3155) -> Result<IntegratedInverseLinkJet, EstimationError> {
3156 match link {
3157 LinkFunction::Log => {
3158 let (mean, saturated) = safe_expwith_saturation(mu + 0.5 * sigma * sigma);
3159 Ok(IntegratedInverseLinkJet {
3160 mean,
3161 d1: mean,
3162 d2: mean,
3163 d3: mean,
3164 mode: if saturated {
3165 IntegratedExpectationMode::ControlledAsymptotic
3166 } else {
3167 IntegratedExpectationMode::ExactClosedForm
3168 },
3169 })
3170 }
3171 LinkFunction::Probit => Ok(integrated_probit_jet(mu, sigma)),
3172 LinkFunction::Logit => {
3173 if sigma > LOGIT_JET_GHQ_SIGMA_MAX {
3174 return logit_wide_sigma_jet(mu, sigma);
3178 }
3179 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
3184 component_point_jet(LinkComponent::Logit, x)
3185 });
3186 let mode = if sigma <= 1e-10 {
3187 IntegratedExpectationMode::ExactClosedForm
3188 } else {
3189 match logit_posterior_meanwith_deriv_controlled(mu, sigma) {
3193 Ok(scalar) => scalar.mode,
3194 Err(_) => IntegratedExpectationMode::QuadratureFallback,
3195 }
3196 };
3197 Ok(IntegratedInverseLinkJet {
3198 mean,
3199 d1: d1.max(0.0),
3200 d2,
3201 d3,
3202 mode,
3203 })
3204 }
3205 LinkFunction::CLogLog => {
3206 validate_latent_cloglog_inputs(mu, sigma)?;
3207 Ok(integrated_cloglog_inverse_link_jet_controlled(
3208 quadctx, mu, sigma,
3209 ))
3210 }
3211 LinkFunction::LogLog | LinkFunction::Cauchit => {
3212 let component = if matches!(link, LinkFunction::LogLog) {
3214 LinkComponent::LogLog
3215 } else {
3216 LinkComponent::Cauchit
3217 };
3218 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(quadctx, mu, sigma, |x| {
3219 component_point_jet(component, x)
3220 });
3221 Ok(IntegratedInverseLinkJet {
3222 mean,
3223 d1,
3224 d2,
3225 d3,
3226 mode: if sigma <= 1e-10 {
3227 IntegratedExpectationMode::ExactClosedForm
3228 } else {
3229 IntegratedExpectationMode::QuadratureFallback
3230 },
3231 })
3232 }
3233 LinkFunction::Sas => Err(EstimationError::InvalidInput(
3234 "state-less integrated SAS jet is unsupported; use SAS-aware prediction APIs with explicit (epsilon, log_delta)".to_string(),
3235 )),
3236 LinkFunction::BetaLogistic => Err(EstimationError::InvalidInput(
3237 "state-less integrated Beta-Logistic jet is unsupported; use link-aware prediction APIs with explicit (delta, epsilon)".to_string(),
3238 )),
3239 LinkFunction::Identity => Ok(IntegratedInverseLinkJet {
3240 mean: mu,
3241 d1: 1.0,
3242 d2: 0.0,
3243 d3: 0.0,
3244 mode: IntegratedExpectationMode::ExactClosedForm,
3245 }),
3246 }
3247}
3248
3249#[inline]
3260fn logit_wide_sigma_jet(mu: f64, sigma: f64) -> Result<IntegratedInverseLinkJet, EstimationError> {
3261 let scalar = logit_posterior_meanwith_deriv_controlled(mu, sigma)?;
3262 let d2 = integrate_normal_adaptive(mu, sigma, |x| {
3263 component_point_jet(LinkComponent::Logit, x).2
3264 });
3265 let d3 = integrate_normal_adaptive(mu, sigma, |x| {
3266 component_point_jet(LinkComponent::Logit, x).3
3267 });
3268 Ok(IntegratedInverseLinkJet {
3269 mean: scalar.mean,
3270 d1: scalar.dmean_dmu.max(0.0),
3271 d2,
3272 d3,
3273 mode: scalar.mode,
3274 })
3275}
3276
3277#[inline]
3278fn sas_point_jet(x: f64, epsilon: f64, log_delta: f64) -> (f64, f64, f64, f64) {
3279 let jet = sas_inverse_link_jet(x, epsilon, log_delta)
3280 .expect("normal quadrature nodes must be finite");
3281 (jet.mu, jet.d1, jet.d2, jet.d3)
3282}
3283
3284#[inline]
3285fn beta_logistic_point_jet(x: f64, log_shape_center: f64, epsilon: f64) -> (f64, f64, f64, f64) {
3286 let jet = beta_logistic_inverse_link_jet(x, log_shape_center, epsilon);
3287 (jet.mu, jet.d1, jet.d2, jet.d3)
3288}
3289
3290#[inline]
3291fn worse_integrated_expectation_mode(
3292 lhs: IntegratedExpectationMode,
3293 rhs: IntegratedExpectationMode,
3294) -> IntegratedExpectationMode {
3295 if lhs.rank() >= rhs.rank() { lhs } else { rhs }
3296}
3297
3298#[inline]
3299fn integrated_scalar_drift_exceeds(
3300 candidate: f64,
3301 reference: f64,
3302 abs_tol: f64,
3303 rel_tol: f64,
3304) -> bool {
3305 if !(candidate.is_finite() && reference.is_finite()) {
3306 return true;
3307 }
3308 (candidate - reference).abs() > abs_tol.max(rel_tol * reference.abs().max(candidate.abs()))
3309}
3310
3311#[inline]
3312fn integrated_mean_derivative_drift_exceeds(
3313 candidate: &IntegratedMeanDerivative,
3314 reference: &IntegratedMeanDerivative,
3315 mean_abs_tol: f64,
3316 mean_rel_tol: f64,
3317 deriv_abs_tol: f64,
3318 deriv_rel_tol: f64,
3319) -> bool {
3320 integrated_scalar_drift_exceeds(candidate.mean, reference.mean, mean_abs_tol, mean_rel_tol)
3321 || integrated_scalar_drift_exceeds(
3322 candidate.dmean_dmu,
3323 reference.dmean_dmu,
3324 deriv_abs_tol,
3325 deriv_rel_tol,
3326 )
3327}
3328
3329#[inline]
3330fn component_point_jet(component: LinkComponent, x: f64) -> (f64, f64, f64, f64) {
3331 let jet = component_inverse_link_jet(component, x);
3334 (jet.mu, jet.d1, jet.d2, jet.d3)
3335}
3336
3337#[inline]
3338fn integrated_mixture_component_jet(
3339 ctx: &QuadratureContext,
3340 component: LinkComponent,
3341 mu: f64,
3342 sigma: f64,
3343) -> IntegratedInverseLinkJet {
3344 match component {
3349 LinkComponent::Logit => integrated_inverse_link_jet(ctx, LinkFunction::Logit, mu, sigma)
3350 .unwrap_or_else(|error| {
3351 log::debug!(
3352 "integrated logit jet at (mu={mu}, sigma={sigma}) fell back to GHQ: {error}"
3353 );
3354 integrated_logit_jet_ghq(ctx, mu, sigma)
3355 }),
3356 LinkComponent::Probit => integrated_probit_jet(mu, sigma),
3357 LinkComponent::CLogLog => integrated_cloglog_inverse_link_jet_controlled(ctx, mu, sigma),
3358 LinkComponent::LogLog | LinkComponent::Cauchit => {
3359 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3360 component_point_jet(component, x)
3361 });
3362 IntegratedInverseLinkJet {
3363 mean,
3364 d1: d1.max(0.0),
3365 d2,
3366 d3,
3367 mode: if sigma <= 1e-10 {
3368 IntegratedExpectationMode::ExactClosedForm
3369 } else {
3370 IntegratedExpectationMode::QuadratureFallback
3371 },
3372 }
3373 }
3374 }
3375}
3376
3377#[inline]
3378fn integrated_mixture_jet(
3379 ctx: &QuadratureContext,
3380 mu: f64,
3381 sigma: f64,
3382 mixture_state: &MixtureLinkState,
3383) -> Result<IntegratedInverseLinkJet, EstimationError> {
3384 if mixture_state.components.is_empty() {
3389 crate::bail_invalid_estim!(
3390 "integrated mixture-link jet requires at least one blended component"
3391 );
3392 }
3393 if mixture_state.components.len() != mixture_state.pi.len() {
3394 crate::bail_invalid_estim!(
3395 "integrated mixture-link jet requires matching component and weight counts"
3396 );
3397 }
3398
3399 let mut mean = 0.0_f64;
3404 let mut d1 = 0.0_f64;
3405 let mut d2 = 0.0_f64;
3406 let mut d3 = 0.0_f64;
3407 let mut mode = IntegratedExpectationMode::ExactClosedForm;
3408 let mut saw_positive_weight = false;
3409
3410 for (&component, &weight) in mixture_state.components.iter().zip(mixture_state.pi.iter()) {
3411 if weight <= 0.0 {
3412 continue;
3413 }
3414 let jet = integrated_mixture_component_jet(ctx, component, mu, sigma);
3415 mean += weight * jet.mean;
3416 d1 += weight * jet.d1;
3417 d2 += weight * jet.d2;
3418 d3 += weight * jet.d3;
3419 if jet.mode.rank() > mode.rank() {
3420 mode = jet.mode;
3421 }
3422 saw_positive_weight = true;
3423 }
3424
3425 if !saw_positive_weight {
3426 crate::bail_invalid_estim!(
3427 "integrated mixture-link jet requires at least one positive component weight"
3428 .to_string(),
3429 );
3430 }
3431
3432 Ok(IntegratedInverseLinkJet {
3433 mean,
3434 d1: d1.max(0.0),
3435 d2,
3436 d3,
3437 mode,
3438 })
3439}
3440
3441#[inline]
3442fn integrated_sas_jet_ghq(
3443 ctx: &QuadratureContext,
3444 mu: f64,
3445 sigma: f64,
3446 sas_state: &SasLinkState,
3447) -> IntegratedInverseLinkJet {
3448 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3449 sas_point_jet(x, sas_state.epsilon, sas_state.log_delta)
3450 });
3451 IntegratedInverseLinkJet {
3452 mean,
3453 d1: d1.max(0.0),
3454 d2,
3455 d3,
3456 mode: if sigma <= 1e-10 {
3457 IntegratedExpectationMode::ExactClosedForm
3458 } else {
3459 IntegratedExpectationMode::QuadratureFallback
3460 },
3461 }
3462}
3463
3464#[inline]
3465fn integrated_beta_logistic_jet_ghq(
3466 ctx: &QuadratureContext,
3467 mu: f64,
3468 sigma: f64,
3469 beta_state: &SasLinkState,
3470) -> IntegratedInverseLinkJet {
3471 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3472 beta_logistic_point_jet(x, beta_state.log_delta, beta_state.epsilon)
3473 });
3474 IntegratedInverseLinkJet {
3475 mean,
3476 d1: d1.max(0.0),
3477 d2,
3478 d3,
3479 mode: if sigma <= 1e-10 {
3480 IntegratedExpectationMode::ExactClosedForm
3481 } else {
3482 IntegratedExpectationMode::QuadratureFallback
3483 },
3484 }
3485}
3486
3487#[inline]
3489pub fn integrated_inverse_link_jetwith_state(
3490 quadctx: &QuadratureContext,
3491 link: LinkFunction,
3492 mu: f64,
3493 sigma: f64,
3494 mixture_link_state: Option<&MixtureLinkState>,
3495 sas_link_state: Option<&SasLinkState>,
3496) -> Result<IntegratedInverseLinkJet, EstimationError> {
3497 if let Some(state) = mixture_link_state {
3498 return integrated_mixture_jet(quadctx, mu, sigma, state);
3499 }
3500 if matches!(link, LinkFunction::Sas) {
3501 let sas = sas_link_state.ok_or_else(|| {
3502 EstimationError::InvalidInput(
3503 "state-less integrated SAS jet is unsupported; explicit SasLinkState is required"
3504 .to_string(),
3505 )
3506 })?;
3507 return Ok(integrated_sas_jet_ghq(quadctx, mu, sigma, sas));
3508 }
3509 if matches!(link, LinkFunction::BetaLogistic) {
3510 let state = sas_link_state.ok_or_else(|| {
3511 EstimationError::InvalidInput(
3512 "state-less integrated Beta-Logistic jet is unsupported; explicit link state is required"
3513 .to_string(),
3514 )
3515 })?;
3516 return Ok(integrated_beta_logistic_jet_ghq(quadctx, mu, sigma, state));
3517 }
3518 integrated_inverse_link_jet(quadctx, link, mu, sigma)
3519}
3520
3521#[inline]
3531pub fn integrated_family_moments_jet(
3532 quadctx: &QuadratureContext,
3533 likelihood: &GlmLikelihoodSpec,
3534 eta: f64,
3535 se_eta: f64,
3536) -> Result<IntegratedMomentsJet, EstimationError> {
3537 const PROB_EPS: f64 = 1e-12;
3538 if !(eta.is_finite() && (-700.0..=700.0).contains(&eta)) {
3539 crate::bail_invalid_estim!(
3540 "integrated moments eta must be finite and within [-700, 700]; got {eta}"
3541 );
3542 }
3543 let e = eta;
3544 let se = se_eta.max(0.0);
3545 let resolved_scale = likelihood
3549 .resolved_scale()
3550 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3551 let spec = &likelihood.spec;
3552 let mixture_link_state: Option<&MixtureLinkState> = spec.link.mixture_state();
3553 let sas_link_state: Option<&SasLinkState> = spec.link.sas_state();
3554 match &spec.response {
3555 ResponseFamily::Binomial => match &spec.link {
3556 InverseLink::Standard(StandardLink::Logit) => {
3557 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3558 let mean = jet.mean;
3559 Ok(IntegratedMomentsJet {
3560 mean,
3561 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3562 d1: jet.d1,
3563 d2: jet.d2,
3564 d3: jet.d3,
3565 mode: jet.mode,
3566 })
3567 }
3568 InverseLink::Standard(StandardLink::Probit) => {
3569 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Probit, e, se)?;
3570 let mean = jet.mean;
3571 Ok(IntegratedMomentsJet {
3572 mean,
3573 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3574 d1: jet.d1,
3575 d2: jet.d2,
3576 d3: jet.d3,
3577 mode: jet.mode,
3578 })
3579 }
3580 InverseLink::Standard(StandardLink::CLogLog) => {
3581 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::CLogLog, e, se)?;
3582 let mean = jet.mean;
3583 Ok(IntegratedMomentsJet {
3584 mean,
3585 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3586 d1: jet.d1,
3587 d2: jet.d2,
3588 d3: jet.d3,
3589 mode: jet.mode,
3590 })
3591 }
3592 InverseLink::LatentCLogLog(_) => Err(EstimationError::InvalidInput(
3593 "Binomial+LatentCLogLog integrated moments require an explicit latent cloglog inverse-link state"
3594 .to_string(),
3595 )),
3596 InverseLink::Sas(_) => {
3597 let jet = integrated_inverse_link_jetwith_state(
3598 quadctx,
3599 LinkFunction::Sas,
3600 e,
3601 se,
3602 mixture_link_state,
3603 sas_link_state,
3604 )?;
3605 let mean = jet.mean;
3606 Ok(IntegratedMomentsJet {
3607 mean,
3608 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3609 d1: jet.d1,
3610 d2: jet.d2,
3611 d3: jet.d3,
3612 mode: jet.mode,
3613 })
3614 }
3615 InverseLink::BetaLogistic(_) => {
3616 let jet = integrated_inverse_link_jetwith_state(
3617 quadctx,
3618 LinkFunction::BetaLogistic,
3619 e,
3620 se,
3621 mixture_link_state,
3622 sas_link_state,
3623 )?;
3624 let mean = jet.mean;
3625 Ok(IntegratedMomentsJet {
3626 mean,
3627 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3628 d1: jet.d1,
3629 d2: jet.d2,
3630 d3: jet.d3,
3631 mode: jet.mode,
3632 })
3633 }
3634 InverseLink::Mixture(state) => {
3635 let jet = integrated_mixture_jet(quadctx, e, se, &state)?;
3636 let mean = jet.mean;
3637 Ok(IntegratedMomentsJet {
3638 mean,
3639 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3640 d1: jet.d1,
3641 d2: jet.d2,
3642 d3: jet.d3,
3643 mode: jet.mode,
3644 })
3645 }
3646 InverseLink::Standard(other) => Err(EstimationError::InvalidInput(format!(
3647 "Binomial response paired with unsupported standard link {other:?} for integrated moments"
3648 ))),
3649 },
3650 ResponseFamily::Gaussian => {
3651 let variance = resolved_scale
3652 .gaussian_phi()
3653 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3654 Ok(IntegratedMomentsJet {
3655 mean: e,
3656 variance,
3657 d1: 1.0,
3658 d2: 0.0,
3659 d3: 0.0,
3660 mode: IntegratedExpectationMode::ExactClosedForm,
3661 })
3662 }
3663 ResponseFamily::RoystonParmar => {
3664 let jet = integrated_inverse_link_jetwith_state(
3665 quadctx,
3666 LinkFunction::CLogLog,
3667 e,
3668 se,
3669 mixture_link_state,
3670 sas_link_state,
3671 )?;
3672 let mean = (1.0 - jet.mean).clamp(0.0, 1.0);
3673 Ok(IntegratedMomentsJet {
3674 mean,
3675 variance: (mean * (1.0 - mean)).max(PROB_EPS),
3676 d1: -jet.d1,
3677 d2: -jet.d2,
3678 d3: -jet.d3,
3679 mode: jet.mode,
3680 })
3681 }
3682 ResponseFamily::Beta { .. } => {
3683 let precision = resolved_scale
3684 .beta_precision()
3685 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3686 let jet = integrated_inverse_link_jet(quadctx, LinkFunction::Logit, e, se)?;
3687 let mean = jet.mean.clamp(PROB_EPS, 1.0 - PROB_EPS);
3688 Ok(IntegratedMomentsJet {
3689 mean,
3690 variance: (mean * (1.0 - mean) / (1.0 + precision)).max(PROB_EPS),
3691 d1: jet.d1,
3692 d2: jet.d2,
3693 d3: jet.d3,
3694 mode: jet.mode,
3695 })
3696 }
3697 ResponseFamily::Poisson
3698 | ResponseFamily::Tweedie { .. }
3699 | ResponseFamily::NegativeBinomial { .. }
3700 | ResponseFamily::Gamma => {
3701 let s2 = se * se;
3706 let (mean, saturated) = safe_expwith_saturation(e + 0.5 * s2);
3707 let variance = match &spec.response {
3718 ResponseFamily::Poisson => mean,
3719 ResponseFamily::Tweedie { p } => {
3720 let phi = resolved_scale
3721 .tweedie_phi()
3722 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3723 phi * mean.powf(*p)
3724 }
3725 ResponseFamily::NegativeBinomial { .. } => {
3726 let theta = resolved_scale.negative_binomial_theta().map_err(|error| {
3727 EstimationError::InvalidInput(error.to_string())
3728 })?;
3729 mean + mean * mean / theta
3730 }
3731 ResponseFamily::Gamma => {
3732 let phi = resolved_scale
3733 .gamma_phi()
3734 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
3735 phi * mean * mean
3736 }
3737 other => {
3741 return Err(EstimationError::InvalidInput(format!(
3742 "integrated log-normal moments reached unexpected family {other:?}"
3743 )));
3744 }
3745 };
3746 if !(variance.is_finite() && variance >= 0.0) {
3747 return Err(EstimationError::InvalidInput(format!(
3748 "integrated {} variance is not representable: {variance:?}",
3749 spec.response.name()
3750 )));
3751 }
3752 Ok(IntegratedMomentsJet {
3753 mean,
3754 variance,
3755 d1: mean,
3756 d2: mean,
3757 d3: mean,
3758 mode: if saturated {
3759 IntegratedExpectationMode::ControlledAsymptotic
3760 } else {
3761 IntegratedExpectationMode::ExactClosedForm
3762 },
3763 })
3764 }
3765 }
3766}
3767
3768pub fn logit_posterior_meanwith_deriv_batch(
3771 ctx: &QuadratureContext,
3772 eta: &ndarray::Array1<f64>,
3773 se_eta: &ndarray::Array1<f64>,
3774) -> Result<(ndarray::Array1<f64>, ndarray::Array1<f64>), EstimationError> {
3775 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3776 let n = eta.len();
3777 let pairs: Result<Vec<(f64, f64)>, _> = (0..n)
3779 .into_par_iter()
3780 .map(|i| {
3781 let integrated = integrated_inverse_link_mean_and_derivative(
3782 ctx,
3783 LinkFunction::Logit,
3784 eta[i],
3785 se_eta[i],
3786 )?;
3787 Ok::<_, EstimationError>((integrated.mean, integrated.dmean_dmu))
3788 })
3789 .collect();
3790 let pairs = pairs?;
3791 let mut mu = ndarray::Array1::<f64>::zeros(n);
3792 let mut dmu = ndarray::Array1::<f64>::zeros(n);
3793 for (i, (m, d)) in pairs.into_iter().enumerate() {
3794 mu[i] = m;
3795 dmu[i] = d;
3796 }
3797
3798 Ok((mu, dmu))
3799}
3800
3801pub fn logit_posterior_mean_batch(
3805 ctx: &QuadratureContext,
3806 eta: &ndarray::Array1<f64>,
3807 se_eta: &ndarray::Array1<f64>,
3808) -> Result<ndarray::Array1<f64>, EstimationError> {
3809 use rayon::iter::{IntoParallelIterator, ParallelIterator};
3810 let n = eta.len();
3811 let values: Result<Vec<f64>, EstimationError> = (0..n)
3812 .into_par_iter()
3813 .map(|i| {
3814 integrated_inverse_link_mean_and_derivative(ctx, LinkFunction::Logit, eta[i], se_eta[i])
3815 .map(|integrated| integrated.mean)
3816 })
3817 .collect();
3818 Ok(ndarray::Array1::from_vec(values?))
3819}
3820
3821pub trait GhqValue: Sized {
3822 fn zero() -> Self;
3823 fn addweighted(&mut self, weight: f64, value: Self);
3824 fn scale(self, factor: f64) -> Self;
3825}
3826
3827impl GhqValue for f64 {
3828 #[inline]
3829 fn zero() -> Self {
3830 0.0
3831 }
3832
3833 #[inline]
3834 fn addweighted(&mut self, weight: f64, value: Self) {
3835 *self += weight * value;
3836 }
3837
3838 #[inline]
3839 fn scale(self, factor: f64) -> Self {
3840 self * factor
3841 }
3842}
3843
3844impl GhqValue for (f64, f64) {
3845 #[inline]
3846 fn zero() -> Self {
3847 (0.0, 0.0)
3848 }
3849
3850 #[inline]
3851 fn addweighted(&mut self, weight: f64, value: Self) {
3852 self.0 += weight * value.0;
3853 self.1 += weight * value.1;
3854 }
3855
3856 #[inline]
3857 fn scale(self, factor: f64) -> Self {
3858 (self.0 * factor, self.1 * factor)
3859 }
3860}
3861
3862impl GhqValue for (f64, f64, f64, f64) {
3863 #[inline]
3864 fn zero() -> Self {
3865 (0.0, 0.0, 0.0, 0.0)
3866 }
3867
3868 #[inline]
3869 fn addweighted(&mut self, weight: f64, value: Self) {
3870 self.0 += weight * value.0;
3871 self.1 += weight * value.1;
3872 self.2 += weight * value.2;
3873 self.3 += weight * value.3;
3874 }
3875
3876 #[inline]
3877 fn scale(self, factor: f64) -> Self {
3878 (
3879 self.0 * factor,
3880 self.1 * factor,
3881 self.2 * factor,
3882 self.3 * factor,
3883 )
3884 }
3885}
3886
3887impl GhqValue for (f64, f64, f64, f64, f64, f64) {
3888 #[inline]
3889 fn zero() -> Self {
3890 (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
3891 }
3892
3893 #[inline]
3894 fn addweighted(&mut self, weight: f64, value: Self) {
3895 self.0 += weight * value.0;
3896 self.1 += weight * value.1;
3897 self.2 += weight * value.2;
3898 self.3 += weight * value.3;
3899 self.4 += weight * value.4;
3900 self.5 += weight * value.5;
3901 }
3902
3903 #[inline]
3904 fn scale(self, factor: f64) -> Self {
3905 (
3906 self.0 * factor,
3907 self.1 * factor,
3908 self.2 * factor,
3909 self.3 * factor,
3910 self.4 * factor,
3911 self.5 * factor,
3912 )
3913 }
3914}
3915
3916#[inline]
3917fn integrate_normal_ghq_adaptive<F, R>(ctx: &QuadratureContext, eta: f64, se_eta: f64, f: F) -> R
3918where
3919 F: Fn(f64) -> R,
3920 R: GhqValue,
3921{
3922 if se_eta < 1e-10 {
3923 return f(eta);
3924 }
3925 let n = adaptive_point_count_from_sd(se_eta.abs());
3926 with_gh_nodesweights(ctx, n, |nodes, weights| {
3927 let scale = SQRT_2 * se_eta;
3928 let mut sum = R::zero();
3929 for i in 0..n {
3930 sum.addweighted(weights[i], f(eta + scale * nodes[i]));
3931 }
3932 sum.scale(1.0 / std::f64::consts::PI.sqrt())
3933 })
3934}
3935
3936#[inline]
3937fn integrated_probit_jet(mu: f64, sigma: f64) -> IntegratedInverseLinkJet {
3938 let s = sigma.hypot(1.0);
3944 let z = mu / s;
3945 let mean = gam_math::probability::normal_cdf(z);
3946 let pdf = gam_math::probability::normal_pdf(z);
3947 if pdf == 0.0 {
3948 return IntegratedInverseLinkJet {
3949 mean,
3950 d1: 0.0,
3951 d2: 0.0,
3952 d3: 0.0,
3953 mode: IntegratedExpectationMode::ExactClosedForm,
3954 };
3955 }
3956 IntegratedInverseLinkJet {
3957 mean,
3958 d1: pdf / s,
3959 d2: -z * pdf / (s * s),
3960 d3: (z * z - 1.0) * pdf / (s * s * s),
3961 mode: IntegratedExpectationMode::ExactClosedForm,
3962 }
3963}
3964
3965#[inline]
3966fn integrated_logit_jet_ghq(
3967 ctx: &QuadratureContext,
3968 mu: f64,
3969 sigma: f64,
3970) -> IntegratedInverseLinkJet {
3971 let (mean, d1, d2, d3) = integrate_normal_ghq_adaptive(ctx, mu, sigma, |x| {
3972 component_point_jet(LinkComponent::Logit, x)
3973 });
3974 IntegratedInverseLinkJet {
3975 mean,
3976 d1: d1.max(0.0),
3977 d2,
3978 d3,
3979 mode: if sigma <= 1e-10 {
3980 IntegratedExpectationMode::ExactClosedForm
3981 } else {
3982 IntegratedExpectationMode::QuadratureFallback
3983 },
3984 }
3985}
3986
3987#[inline]
3988fn cloglog_inverse_link_controlled_values(
3989 ctx: &QuadratureContext,
3990 mu: f64,
3991 sigma: f64,
3992 max_order: usize,
3993) -> ([f64; 6], IntegratedExpectationMode) {
3994 assert!(max_order <= 5);
3995 if sigma <= 1e-10 {
3996 let (mean, d1, d2, d3, d4, d5) = cloglog_point_jet5(mu);
3997 return (
3998 [mean, d1, d2, d3, d4, d5],
3999 IntegratedExpectationMode::ExactClosedForm,
4000 );
4001 }
4002
4003 let (k, log_k0, mode) = latent_cloglog_kernel_terms(ctx, mu, sigma, max_order);
4004 let mut values = [0.0; 6];
4005 values[0] = if log_k0.is_finite() {
4006 -log_k0.exp_m1()
4007 } else {
4008 1.0
4009 };
4010 values[1] = k[1].max(0.0);
4011 if sigma > CLOGLOG_JET_MOMENT_SIGMA_MAX {
4012 if max_order >= 2 {
4013 values[2] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).2);
4014 }
4015 if max_order >= 3 {
4016 values[3] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).3);
4017 }
4018 if max_order >= 4 {
4019 values[4] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).4);
4020 }
4021 if max_order >= 5 {
4022 values[5] = integrate_normal_adaptive(mu, sigma, |x| cloglog_point_jet5(x).5);
4023 }
4024 return (
4025 values,
4026 worse_integrated_expectation_mode(mode, IntegratedExpectationMode::QuadratureFallback),
4027 );
4028 }
4029 if max_order >= 2 {
4030 values[2] = k[1] - k[2];
4031 }
4032 if max_order >= 3 {
4033 values[3] = k[1] - 3.0 * k[2] + k[3];
4034 }
4035 if max_order >= 4 {
4036 values[4] = k[1] - 7.0 * k[2] + 6.0 * k[3] - k[4];
4037 }
4038 if max_order >= 5 {
4039 values[5] = k[1] - 15.0 * k[2] + 25.0 * k[3] - 10.0 * k[4] + k[5];
4040 }
4041 (values, mode)
4042}
4043
4044#[inline]
4045pub(crate) fn latent_cloglog_inverse_link_jet5_controlled(
4046 ctx: &QuadratureContext,
4047 mu: f64,
4048 sigma: f64,
4049) -> IntegratedInverseLinkJet5 {
4050 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 5);
4051 IntegratedInverseLinkJet5 {
4052 mean: values[0],
4053 d1: values[1],
4054 d2: values[2],
4055 d3: values[3],
4056 d4: values[4],
4057 d5: values[5],
4058 mode,
4059 }
4060}
4061
4062#[derive(Clone, Copy, Debug)]
4072pub struct LatentCLogLogJet5 {
4073 pub mean: f64,
4074 pub d1: f64,
4075 pub d2: f64,
4076 pub d3: f64,
4077 pub d4: f64,
4078 pub d5: f64,
4079 pub mode: IntegratedExpectationMode,
4080}
4081
4082pub fn latent_cloglog_jet5(
4083 quadctx: &QuadratureContext,
4084 eta: f64,
4085 sigma: f64,
4086) -> Result<LatentCLogLogJet5, EstimationError> {
4087 validate_latent_cloglog_inputs(eta, sigma)?;
4088 let jet = latent_cloglog_inverse_link_jet5_controlled(quadctx, eta, sigma);
4094 Ok(LatentCLogLogJet5 {
4095 mean: jet.mean,
4096 d1: jet.d1,
4097 d2: jet.d2,
4098 d3: jet.d3,
4099 d4: jet.d4,
4100 d5: jet.d5,
4101 mode: jet.mode,
4102 })
4103}
4104
4105#[inline]
4106pub fn latent_cloglog_inverse_link_jet(
4107 quadctx: &QuadratureContext,
4108 eta: f64,
4109 sigma: f64,
4110) -> Result<IntegratedInverseLinkJet, EstimationError> {
4111 let jet = latent_cloglog_jet5(quadctx, eta, sigma)?;
4112 Ok(IntegratedInverseLinkJet {
4113 mean: jet.mean,
4114 d1: jet.d1,
4115 d2: jet.d2,
4116 d3: jet.d3,
4117 mode: jet.mode,
4118 })
4119}
4120
4121#[inline]
4122fn integrated_cloglog_inverse_link_jet_controlled(
4123 ctx: &QuadratureContext,
4124 mu: f64,
4125 sigma: f64,
4126) -> IntegratedInverseLinkJet {
4127 let (values, mode) = cloglog_inverse_link_controlled_values(ctx, mu, sigma, 3);
4128 IntegratedInverseLinkJet {
4129 mean: values[0],
4130 d1: values[1],
4131 d2: values[2],
4132 d3: values[3],
4133 mode,
4134 }
4135}
4136
4137#[inline]
4138fn latent_cloglog_kernel_terms(
4139 ctx: &QuadratureContext,
4140 mu: f64,
4141 sigma: f64,
4142 max_order: usize,
4143) -> ([f64; 6], f64, IntegratedExpectationMode) {
4144 let sigma2 = sigma * sigma;
4145 let mut k = [0.0; 6];
4146 let mut log_k0 = f64::NEG_INFINITY;
4147 let mut mode = IntegratedExpectationMode::ExactClosedForm;
4148
4149 for (order, out) in k.iter_mut().enumerate().take(max_order + 1) {
4150 let kf = order as f64;
4151 let shifted_mu = mu + kf * sigma2;
4152 let (log_survival, term_mode) =
4160 cloglog_log_survival_term_controlled(ctx, shifted_mu, sigma);
4161 mode = worse_integrated_expectation_mode(mode, term_mode);
4162
4163 let log_value = kf * mu + 0.5 * kf * kf * sigma2 + log_survival;
4164 if order == 0 {
4165 log_k0 = log_value;
4166 }
4167 if !log_value.is_finite() {
4168 *out = 0.0;
4169 continue;
4170 }
4171 let upper = if order == 0 {
4172 1.0
4173 } else {
4174 let k_over_e = kf / std::f64::consts::E;
4175 k_over_e.powf(kf)
4176 };
4177 *out = safe_exp(log_value).clamp(0.0, upper);
4178 }
4179
4180 (k, log_k0, mode)
4181}
4182
4183#[inline]
4184pub fn normal_expectation_1d_adaptive<F>(
4185 ctx: &QuadratureContext,
4186 eta: f64,
4187 se_eta: f64,
4188 f: F,
4189) -> f64
4190where
4191 F: Fn(f64) -> f64,
4192{
4193 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
4194}
4195
4196#[inline]
4197pub fn normal_expectation_1d_adaptive_pair<F>(
4198 ctx: &QuadratureContext,
4199 eta: f64,
4200 se_eta: f64,
4201 f: F,
4202) -> (f64, f64)
4203where
4204 F: Fn(f64) -> (f64, f64),
4205{
4206 integrate_normal_ghq_adaptive(ctx, eta, se_eta, f)
4207}
4208
4209fn adaptive_point_count_from_sd(max_sd: f64) -> usize {
4210 if max_sd.is_finite() && max_sd > 2.5 {
4220 51
4221 } else if max_sd.is_finite() && max_sd > 0.5 {
4222 31
4223 } else if max_sd.is_finite() && max_sd > 0.35 {
4224 21
4225 } else if max_sd.is_finite() && max_sd > 0.1 {
4226 15
4227 } else {
4228 7
4229 }
4230}
4231
4232#[inline]
4233fn with_gh_nodesweights<R>(
4234 ctx: &QuadratureContext,
4235 n: usize,
4236 f: impl FnOnce(&[f64], &[f64]) -> R,
4237) -> R {
4238 if n == 7 {
4239 let gh = ctx.gauss_hermite();
4240 f(&gh.nodes, &gh.weights)
4241 } else {
4242 let gh = ctx.gauss_hermite_n(n);
4243 f(&gh.nodes, &gh.weights)
4244 }
4245}
4246
4247#[inline]
4257fn cholesky_static<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
4258 let mut l = [[0.0_f64; D]; D];
4259 for i in 0..D {
4260 for j in 0..=i {
4261 let mut sum = cov[i][j];
4262 for k in 0..j {
4263 sum -= l[i][k] * l[j][k];
4264 }
4265 if i == j {
4266 if !sum.is_finite() || sum <= 0.0 {
4267 return None;
4268 }
4269 l[i][j] = sum.sqrt();
4270 } else {
4271 l[i][j] = sum / l[j][j];
4272 }
4273 }
4274 }
4275 Some(l)
4276}
4277
4278#[inline]
4281fn cholesky_static_with_jitter<const D: usize>(cov: &[[f64; D]; D]) -> Option<[[f64; D]; D]> {
4282 if D == 0 {
4283 return None;
4284 }
4285 for retry in 0..8 {
4286 let jitter = if retry == 0 {
4287 0.0
4288 } else {
4289 1e-12 * 10f64.powi(retry - 1)
4290 };
4291 if jitter == 0.0 {
4292 if let Some(l) = cholesky_static::<D>(cov) {
4293 return Some(l);
4294 }
4295 } else {
4296 let mut base = *cov;
4297 for i in 0..D {
4298 base[i][i] = cov[i][i] + jitter;
4299 }
4300 if let Some(l) = cholesky_static::<D>(&base) {
4301 return Some(l);
4302 }
4303 }
4304 }
4305 None
4306}
4307
4308#[inline]
4309fn adaptive_point_countwith_cap(max_sd: f64, max_n: usize) -> usize {
4310 adaptive_point_count_from_sd(max_sd).min(max_n)
4311}
4312
4313#[inline]
4314fn ghq_nd_integrate_try<const D: usize, F, R, E>(
4315 ctx: &QuadratureContext,
4316 mu: [f64; D],
4317 cov: [[f64; D]; D],
4318 max_n: usize,
4319 f: F,
4320) -> Result<Option<R>, E>
4321where
4322 F: Fn([f64; D]) -> Result<R, E>,
4323 R: GhqValue,
4324{
4325 let mut maxvar = 0.0_f64;
4326 for (i, row) in cov.iter().enumerate() {
4327 maxvar = maxvar.max(row[i]).max(0.0);
4328 }
4329 let n = adaptive_point_countwith_cap(maxvar.sqrt(), max_n);
4330
4331 let mut cov_arr = cov;
4336 for i in 0..D {
4337 cov_arr[i][i] = cov_arr[i][i].max(0.0);
4338 }
4339 let Some(l) = cholesky_static_with_jitter::<D>(&cov_arr) else {
4340 return Ok(None);
4341 };
4342 let norm = 1.0 / std::f64::consts::PI.powf(0.5 * D as f64);
4343
4344 with_gh_nodesweights(ctx, n, |nodes, weights| {
4345 let mut acc = R::zero();
4346 let mut idx = [0usize; D];
4347 loop {
4348 let mut z = [0.0_f64; D];
4349 let mut weight = 1.0_f64;
4350 for d in 0..D {
4351 z[d] = SQRT_2 * nodes[idx[d]];
4352 weight *= weights[idx[d]];
4353 }
4354
4355 let mut x = mu;
4356 for row in 0..D {
4357 let mut dot = 0.0_f64;
4358 for (col, zc) in z.iter().enumerate().take(row + 1) {
4359 dot += l[row][col] * *zc;
4360 }
4361 x[row] += dot;
4362 }
4363 acc.addweighted(weight, f(x)?);
4364
4365 let mut carry = true;
4366 for d in (0..D).rev() {
4367 idx[d] += 1;
4368 if idx[d] < n {
4369 carry = false;
4370 break;
4371 }
4372 idx[d] = 0;
4373 }
4374 if carry {
4375 break;
4376 }
4377 }
4378 Ok(Some(acc.scale(norm)))
4379 })
4380}
4381
4382#[inline]
4383fn ghq_nd_integrate<const D: usize, F, R>(
4384 ctx: &QuadratureContext,
4385 mu: [f64; D],
4386 cov: [[f64; D]; D],
4387 max_n: usize,
4388 f: F,
4389) -> Option<R>
4390where
4391 F: Fn([f64; D]) -> R,
4392 R: GhqValue,
4393{
4394 match ghq_nd_integrate_try::<D, _, R, Infallible>(ctx, mu, cov, max_n, |x| Ok(f(x))) {
4395 Ok(v) => v,
4396 Err(e) => match e {},
4397 }
4398}
4399
4400#[inline]
4401fn ghq_nd_integrate_result<const D: usize, F, R, E>(
4402 ctx: &QuadratureContext,
4403 mu: [f64; D],
4404 cov: [[f64; D]; D],
4405 max_n: usize,
4406 f: F,
4407) -> Result<Option<R>, E>
4408where
4409 F: Fn([f64; D]) -> Result<R, E>,
4410 R: GhqValue,
4411{
4412 ghq_nd_integrate_try::<D, _, R, E>(ctx, mu, cov, max_n, f)
4413}
4414
4415pub fn normal_expectation_nd_adaptive<const D: usize, F>(
4417 ctx: &QuadratureContext,
4418 mu: [f64; D],
4419 cov: [[f64; D]; D],
4420 max_n: usize,
4421 f: F,
4422) -> f64
4423where
4424 F: Fn([f64; D]) -> f64,
4425{
4426 match ghq_nd_integrate::<D, _, f64>(ctx, mu, cov, max_n, &f) {
4427 Some(v) => v,
4428 None => f(mu),
4429 }
4430}
4431
4432pub fn normal_expectation_nd_adaptive_result<const D: usize, F, R, E>(
4434 ctx: &QuadratureContext,
4435 mu: [f64; D],
4436 cov: [[f64; D]; D],
4437 max_n: usize,
4438 f: F,
4439) -> Result<R, E>
4440where
4441 F: Fn([f64; D]) -> Result<R, E>,
4442 R: GhqValue,
4443{
4444 match ghq_nd_integrate_result::<D, _, R, E>(ctx, mu, cov, max_n, &f)? {
4445 Some(v) => Ok(v),
4446 None => f(mu),
4447 }
4448}
4449
4450pub fn normal_expectation_2d_adaptive_result<F, E>(
4452 ctx: &QuadratureContext,
4453 mu: [f64; 2],
4454 cov: [[f64; 2]; 2],
4455 f: F,
4456) -> Result<f64, E>
4457where
4458 F: Fn(f64, f64) -> Result<f64, E>,
4459{
4460 normal_expectation_nd_adaptive_result::<2, _, _, E>(ctx, mu, cov, 21, |x| f(x[0], x[1]))
4461}
4462
4463pub fn normal_expectation_3d_adaptive<F>(
4465 ctx: &QuadratureContext,
4466 mu: [f64; 3],
4467 cov: [[f64; 3]; 3],
4468 f: F,
4469) -> f64
4470where
4471 F: Fn(f64, f64, f64) -> f64,
4472{
4473 normal_expectation_nd_adaptive::<3, _>(ctx, mu, cov, 15, |x| f(x[0], x[1], x[2]))
4475}
4476
4477#[inline]
4496pub fn probit_posterior_mean(eta: f64, se_eta: f64) -> f64 {
4497 if se_eta < 1e-10 {
4498 return gam_math::probability::normal_cdf(eta);
4499 }
4500 let denom = (1.0 + se_eta * se_eta).sqrt();
4501 gam_math::probability::normal_cdf(eta / denom)
4502}
4503
4504#[inline]
4505pub fn logit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4506 let (m1, m2) = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4507 let p = sigmoid(x);
4508 (p, p * p)
4509 });
4510 let m1 = m1.clamp(0.0, 1.0);
4511 let m2 = m2.clamp(0.0, 1.0);
4512 (m1, (m2 - m1 * m1).max(0.0))
4513}
4514
4515#[inline]
4516pub fn probit_posterior_meanvariance(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> (f64, f64) {
4517 let m1 = probit_posterior_mean(eta, se_eta);
4518 let m2 = integrate_normal_ghq_adaptive(ctx, eta, se_eta, |x| {
4519 let p = gam_math::probability::normal_cdf(x);
4520 p * p
4521 })
4522 .clamp(0.0, 1.0);
4523 (m1, (m2 - m1 * m1).max(0.0))
4524}
4525
4526#[inline]
4527pub fn cloglog_posterior_meanvariance(
4528 ctx: &QuadratureContext,
4529 eta: f64,
4530 se_eta: f64,
4531) -> (f64, f64) {
4532 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4552 return (cloglog_mean_exact(eta), 0.0);
4553 }
4554 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4555 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4556 let mean = cloglog_mean_from_survival(survival);
4557 let variance = (survival_sq - survival * survival).max(0.0);
4558 (mean, variance)
4559}
4560
4561#[inline]
4595pub fn cloglog_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4596 if !(eta.is_finite() && se_eta.is_finite()) || se_eta <= CLOGLOG_SIGMA_DEGENERATE {
4600 return cloglog_mean_exact(eta);
4601 }
4602 let (survival, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4603 cloglog_mean_from_survival(survival)
4604}
4605
4606#[inline]
4620pub fn survival_posterior_mean(ctx: &QuadratureContext, eta: f64, se_eta: f64) -> f64 {
4621 cloglog_survival_term_controlled(ctx, eta, se_eta)
4622 .0
4623 .clamp(0.0, 1.0)
4624}
4625
4626#[inline]
4627pub fn survival_posterior_meanvariance(
4628 ctx: &QuadratureContext,
4629 eta: f64,
4630 se_eta: f64,
4631) -> (f64, f64) {
4632 let (m1, _) = cloglog_survival_term_controlled(ctx, eta, se_eta);
4633 let (m2, _) = cloglog_survivalsecond_moment_controlled(ctx, eta, se_eta);
4634 (m1.clamp(0.0, 1.0), (m2 - m1 * m1).max(0.0))
4635}
4636
4637pub fn logit_posterior_mean_exact(mu: f64, sigma: f64) -> f64 {
4713 if !(mu.is_finite() && sigma.is_finite()) || sigma <= 0.0 {
4714 return sigmoid(mu);
4715 }
4716 if sigma < LOGIT_SIGMA_DEGENERATE {
4717 return sigmoid(mu);
4720 }
4721
4722 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI; let sqrt2_sigma = SQRT_2 * sigma;
4724 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;
4728
4729 let mut corr = 0.0_f64;
4735 let mut n = 1usize;
4736 let tail_start = loop {
4737 let b = (2.0 * (n as f64) - 1.0) * beta;
4738 let abs_xi2 = c * c + b * b;
4739 if abs_xi2 > r2 && n >= FADDEEVA_TAIL_MIN_INDEX {
4740 break n;
4741 }
4742 let xi = Complex { re: c, im: b };
4743 let d = if abs_xi2 > r2 {
4744 inv_sqrt_pi * faddeeva_asymptotic_a(xi).re
4746 } else {
4747 faddeeva_upper_halfplane(xi).im - inv_sqrt_pi * c / abs_xi2
4748 };
4749 corr += d;
4750 n += 1;
4751 };
4752
4753 corr += faddeeva_pole_series_em_tail(c, beta, tail_start, inv_sqrt_pi);
4754
4755 sigmoid(mu) - coeff * corr
4756}
4757
4758const FADDEEVA_TAIL_MIN_INDEX: usize = 48;
4762const FADDEEVA_ASYMPTOTIC_RADIUS: f64 = 7.0;
4765const FADDEEVA_ASYMPTOTIC_TERMS: usize = 14;
4768
4769fn faddeeva_asymptotic_a(xi: Complex) -> Complex {
4773 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4774 let inv2 = complexmul(inv, inv);
4775 let mut xp = complexmul(inv2, inv); let mut cm = 0.5_f64; let mut s = Complex::default();
4778 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4779 s = complex_add(
4780 s,
4781 Complex {
4782 re: cm * xp.re,
4783 im: cm * xp.im,
4784 },
4785 );
4786 cm *= (2.0 * (m as f64) + 1.0) / 2.0; xp = complexmul(xp, inv2);
4788 }
4789 s
4790}
4791
4792fn faddeeva_pole_series_em_tail(c: f64, beta: f64, tail_start: usize, inv_sqrt_pi: f64) -> f64 {
4801 let b_a = (2.0 * (tail_start as f64) - 1.0) * beta;
4802 let xi = Complex { re: c, im: b_a };
4803 let inv = complex_div(Complex { re: 1.0, im: 0.0 }, xi);
4804 let inv2 = complexmul(inv, inv);
4805 let two_i_beta = Complex {
4807 re: 0.0,
4808 im: 2.0 * beta,
4809 };
4810
4811 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;
4819 for m in 1..=FADDEEVA_ASYMPTOTIC_TERMS {
4820 let mf = m as f64;
4821 let inv_4ibm = Complex {
4823 re: 0.0,
4824 im: -1.0 / (4.0 * beta * mf),
4825 };
4826 s = complex_add(
4827 s,
4828 complexmul(
4829 Complex {
4830 re: cm * x2m.re,
4831 im: cm * x2m.im,
4832 },
4833 inv_4ibm,
4834 ),
4835 );
4836 a_acc = complex_add(
4837 a_acc,
4838 Complex {
4839 re: cm * x2m1.re,
4840 im: cm * x2m1.im,
4841 },
4842 );
4843 let fc = cm * (-(2.0 * mf + 1.0));
4844 fp_inner = complex_add(
4845 fp_inner,
4846 Complex {
4847 re: fc * x2m2.re,
4848 im: fc * x2m2.im,
4849 },
4850 );
4851 cm *= (2.0 * mf + 1.0) / 2.0;
4852 x2m = complexmul(x2m, inv2);
4853 x2m1 = complexmul(x2m1, inv2);
4854 x2m2 = complexmul(x2m2, inv2);
4855 }
4856
4857 s = complex_add(
4859 s,
4860 Complex {
4861 re: 0.5 * a_acc.re,
4862 im: 0.5 * a_acc.im,
4863 },
4864 );
4865 let fprime = complexmul(two_i_beta, fp_inner);
4867 s = complex_add(
4868 s,
4869 Complex {
4870 re: -fprime.re / 12.0,
4871 im: -fprime.im / 12.0,
4872 },
4873 );
4874
4875 inv_sqrt_pi * s.re
4878}
4879
4880fn faddeeva_upper_halfplane(z: Complex) -> Complex {
4893 let (l, coeffs) = faddeeva_weideman_coeffs();
4894 let iz = Complex {
4895 re: -z.im,
4896 im: z.re,
4897 }; let l_minus = Complex {
4899 re: l - iz.re,
4900 im: -iz.im,
4901 }; let l_plus = Complex {
4903 re: l + iz.re,
4904 im: iz.im,
4905 }; let zz = complex_div(l_plus, l_minus); let mut p = Complex {
4909 re: coeffs[0],
4910 im: 0.0,
4911 };
4912 for &c in &coeffs[1..] {
4913 p = complex_add(complexmul(p, zz), Complex { re: c, im: 0.0 });
4914 }
4915 let l_minus_sq = complexmul(l_minus, l_minus);
4916 let term1 = complex_div(
4917 Complex {
4918 re: 2.0 * p.re,
4919 im: 2.0 * p.im,
4920 },
4921 l_minus_sq,
4922 );
4923 let inv_sqrt_pi = 0.5 * std::f64::consts::FRAC_2_SQRT_PI;
4924 let term2 = complex_div(
4925 Complex {
4926 re: inv_sqrt_pi,
4927 im: 0.0,
4928 },
4929 l_minus,
4930 );
4931 complex_add(term1, term2)
4932}
4933
4934const FADDEEVA_WEIDEMAN_N: usize = 44;
4937
4938fn faddeeva_weideman_coeffs() -> &'static (f64, [f64; FADDEEVA_WEIDEMAN_N]) {
4944 static CACHE: OnceLock<(f64, [f64; FADDEEVA_WEIDEMAN_N])> = OnceLock::new();
4945 CACHE.get_or_init(|| {
4946 let n = FADDEEVA_WEIDEMAN_N;
4947 let l = (n as f64 / SQRT_2).sqrt();
4948 let m = 2 * n;
4949 let m2 = 2 * m; let mut f = vec![0.0_f64; m2];
4953 for (idx, fi) in f.iter_mut().enumerate().skip(1) {
4954 let k = (idx as isize - 1) - (m as isize - 1);
4955 let theta = (k as f64) * std::f64::consts::PI / (m as f64);
4956 let t = l * (0.5 * theta).tan();
4957 *fi = (-t * t).exp() * (l * l + t * t);
4958 }
4959 let half = m2 / 2;
4962 let mut coeffs = [0.0_f64; FADDEEVA_WEIDEMAN_N];
4963 for j in 1..=n {
4964 let mut acc = 0.0_f64;
4965 for (p, _) in f.iter().enumerate() {
4966 let fp = f[(p + half) % m2];
4967 if fp != 0.0 {
4968 acc += fp
4969 * (-2.0 * std::f64::consts::PI * (j as f64) * (p as f64) / (m2 as f64))
4970 .cos();
4971 }
4972 }
4973 coeffs[n - j] = acc / (m2 as f64);
4975 }
4976 (l, coeffs)
4977 })
4978}
4979
4980#[inline]
4982fn sigmoid(x: f64) -> f64 {
4983 let x_clamped = x.clamp(-QUADRATURE_EXP_LOG_MAX, QUADRATURE_EXP_LOG_MAX);
4984 1.0 / (1.0 + f64::exp(-x_clamped))
4985}
4986
4987#[derive(Clone, Copy, Debug)]
5003pub struct CLogLogConvolutionDerivatives {
5004 pub l: f64,
5006
5007 pub l_mu: f64,
5009 pub l_sigma: f64,
5010
5011 pub l_mumu: f64,
5013 pub l_musigma: f64,
5014 pub l_sigmasigma: f64,
5015
5016 pub l_mumumu: f64,
5018 pub l_mumusigma: f64,
5019 pub l_musigmasigma: f64,
5020 pub l_sigmasigmasigma: f64,
5021
5022 pub l_mumumumu: f64,
5024 pub l_mumumusigma: f64,
5025 pub l_mumusigmasigma: f64,
5026 pub l_musigmasigmasigma: f64,
5027 pub l_sigmasigmasigmasigma: f64,
5028}
5029
5030#[inline]
5031pub(crate) fn cloglog_point_jet5(t: f64) -> (f64, f64, f64, f64, f64, f64) {
5032 if t.is_nan() {
5033 return (f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN);
5034 }
5035 let et = safe_exp(t);
5036
5037 (
5038 -(-et).exp_m1(),
5039 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0]),
5040 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -1.0]),
5041 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -3.0, 1.0]),
5042 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -7.0, 6.0, -1.0]),
5043 cloglog_stable_poly_times_exp_neg(et, &[0.0, 1.0, -15.0, 25.0, -10.0, 1.0]),
5044 )
5045}
5046
5047#[inline]
5059fn cloglog_g_derivatives(t: f64) -> (f64, f64, f64, f64, f64) {
5060 let (g, g1, g2, g3, g4, _) = cloglog_point_jet5(t);
5061 (g, g1, g2, g3, g4)
5062}
5063
5064pub fn cloglog_ghq_value(ctx: &QuadratureContext, mu: f64, sigma: f64, n_nodes: usize) -> f64 {
5072 if sigma.abs() < 1e-14 {
5073 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
5074 return g.clamp(0.0, 1.0);
5075 }
5076 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
5077
5078 let inv_sig2 = 1.0 / (sigma * sigma);
5110 let mut eta_hat = mu;
5111 let mut converged = false;
5112 for _ in 0..100 {
5113 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
5114 if !(g > 0.0) || !g1.is_finite() || !g2.is_finite() {
5115 break;
5116 }
5117 let r = g1 / g;
5118 let lp = r - (eta_hat - mu) * inv_sig2;
5119 let lpp = g2 / g - r * r - inv_sig2;
5120 if !lpp.is_finite() || lpp >= 0.0 {
5121 break;
5122 }
5123 let step = lp / lpp;
5124 eta_hat -= step;
5125 if step.abs() <= 1e-13 * (1.0 + eta_hat.abs()) {
5126 converged = true;
5127 break;
5128 }
5129 }
5130
5131 let tau = if converged {
5135 let (g, g1, g2, _, _, _) = cloglog_point_jet5(eta_hat);
5136 if g > 0.0 {
5137 let r = g1 / g;
5138 let lpp = g2 / g - r * r - inv_sig2;
5139 let tau2 = -1.0 / lpp;
5140 if tau2.is_finite() && tau2 > 0.0 {
5141 Some(tau2.sqrt())
5142 } else {
5143 None
5144 }
5145 } else {
5146 None
5147 }
5148 } else {
5149 None
5150 };
5151
5152 let eval_at = |n: usize| -> f64 {
5154 match tau {
5155 Some(tau) => {
5156 let pref = tau * inv_sqrt_pi / sigma;
5157 with_gh_nodesweights(ctx, n, |nodes, weights| {
5158 let mut sum = 0.0_f64;
5159 for i in 0..nodes.len() {
5160 let t = nodes[i];
5161 let eta_i = eta_hat + SQRT_2 * tau * t;
5162 let (g, _, _, _, _, _) = cloglog_point_jet5(eta_i);
5163 let dev = eta_i - mu;
5164 sum += weights[i] * (t * t - 0.5 * dev * dev * inv_sig2).exp() * g;
5165 }
5166 (pref * sum).clamp(0.0, 1.0)
5167 })
5168 }
5169 None => {
5170 let scale = SQRT_2 * sigma;
5171 with_gh_nodesweights(ctx, n, |nodes, weights| {
5172 let mut sum = 0.0_f64;
5173 for i in 0..nodes.len() {
5174 let t = mu + scale * nodes[i];
5175 let (g, _, _, _, _) = cloglog_g_derivatives(t);
5176 sum += weights[i] * g;
5177 }
5178 (sum * inv_sqrt_pi).clamp(0.0, 1.0)
5179 })
5180 }
5181 }
5182 };
5183
5184 const CLOGLOG_GHQ_ORDER_LADDER: [usize; 5] = [7, 15, 21, 31, 51];
5189 const CLOGLOG_GHQ_CONV_TOL: f64 = 1e-10;
5190 let floor = n_nodes.min(CLOGLOG_GHQ_ORDER_LADDER[CLOGLOG_GHQ_ORDER_LADDER.len() - 1]);
5191 let mut prev: Option<f64> = None;
5192 let mut result = 0.0_f64;
5193 for &n in CLOGLOG_GHQ_ORDER_LADDER.iter().filter(|&&n| n >= floor) {
5194 let cur = eval_at(n);
5195 result = cur;
5196 if let Some(p) = prev
5197 && (cur - p).abs() < CLOGLOG_GHQ_CONV_TOL
5198 {
5199 break;
5200 }
5201 prev = Some(cur);
5202 }
5203 result
5204}
5205
5206pub fn cloglog_ghq_derivatives(
5217 ctx: &QuadratureContext,
5218 mu: f64,
5219 sigma: f64,
5220 n_nodes: usize,
5221) -> CLogLogConvolutionDerivatives {
5222 let inv_sqrt_pi = 1.0 / std::f64::consts::PI.sqrt();
5223
5224 if sigma.abs() < 1e-14 {
5231 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
5232 return CLogLogConvolutionDerivatives {
5233 l: g,
5234 l_mu: g1,
5235 l_sigma: 0.0,
5236 l_mumu: g2,
5237 l_musigma: 0.0,
5238 l_sigmasigma: g2,
5239 l_mumumu: g3,
5240 l_mumusigma: 0.0,
5241 l_musigmasigma: g3,
5242 l_sigmasigmasigma: 0.0,
5243 l_mumumumu: g4,
5244 l_mumumusigma: 0.0,
5245 l_mumusigmasigma: g4,
5246 l_musigmasigmasigma: 0.0,
5247 l_sigmasigmasigmasigma: 3.0 * g4,
5248 };
5249 }
5250
5251 let scale = SQRT_2 * sigma;
5252 let sqrt2 = SQRT_2;
5253
5254 with_gh_nodesweights(ctx, n_nodes, |nodes, weights| {
5255 let mut s = [[0.0_f64; 5]; 5];
5267
5268 for i in 0..nodes.len() {
5269 let x = nodes[i];
5270 let t = mu + scale * x;
5271 let (g0, g1, g2, g3, g4) = cloglog_g_derivatives(t);
5272 let w = weights[i];
5273
5274 let x2 = x * x;
5276 let x3 = x2 * x;
5277 let x4 = x3 * x;
5278
5279 s[0][0] += w * g0;
5281
5282 s[1][0] += w * g1;
5284 s[1][1] += w * x * g1;
5285
5286 s[2][0] += w * g2;
5288 s[2][1] += w * x * g2;
5289 s[2][2] += w * x2 * g2;
5290
5291 s[3][0] += w * g3;
5293 s[3][1] += w * x * g3;
5294 s[3][2] += w * x2 * g3;
5295 s[3][3] += w * x3 * g3;
5296
5297 s[4][0] += w * g4;
5299 s[4][1] += w * x * g4;
5300 s[4][2] += w * x2 * g4;
5301 s[4][3] += w * x3 * g4;
5302 s[4][4] += w * x4 * g4;
5303 }
5304
5305 let sqrt2_1 = sqrt2;
5308 let sqrt2_2 = 2.0; let sqrt2_3 = 2.0 * sqrt2; let sqrt2_4 = 4.0; CLogLogConvolutionDerivatives {
5313 l: inv_sqrt_pi * s[0][0],
5315
5316 l_mu: inv_sqrt_pi * s[1][0],
5318 l_sigma: inv_sqrt_pi * sqrt2_1 * s[1][1],
5319
5320 l_mumu: inv_sqrt_pi * s[2][0],
5322 l_musigma: inv_sqrt_pi * sqrt2_1 * s[2][1],
5323 l_sigmasigma: inv_sqrt_pi * sqrt2_2 * s[2][2],
5324
5325 l_mumumu: inv_sqrt_pi * s[3][0],
5327 l_mumusigma: inv_sqrt_pi * sqrt2_1 * s[3][1],
5328 l_musigmasigma: inv_sqrt_pi * sqrt2_2 * s[3][2],
5329 l_sigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[3][3],
5330
5331 l_mumumumu: inv_sqrt_pi * s[4][0],
5333 l_mumumusigma: inv_sqrt_pi * sqrt2_1 * s[4][1],
5334 l_mumusigmasigma: inv_sqrt_pi * sqrt2_2 * s[4][2],
5335 l_musigmasigmasigma: inv_sqrt_pi * sqrt2_3 * s[4][3],
5336 l_sigmasigmasigmasigma: inv_sqrt_pi * sqrt2_4 * s[4][4],
5337 }
5338 })
5339}
5340
5341pub fn cloglog_ghq_derivatives_adaptive(
5347 ctx: &QuadratureContext,
5348 mu: f64,
5349 sigma: f64,
5350) -> CLogLogConvolutionDerivatives {
5351 let n = adaptive_point_count_from_sd(sigma.abs());
5352 cloglog_ghq_derivatives(ctx, mu, sigma, n)
5353}
5354
5355#[cfg(test)]
5356mod tests {
5357 use super::*;
5358 use approx::assert_relative_eq;
5359 use gam_problem::LikelihoodScaleMetadata;
5360 use gam_spec::LikelihoodSpec;
5361
5362 #[test]
5369 fn log_half_erfc_stable_matches_high_precision_reference() {
5370 let refs: &[(f64, f64)] = &[
5371 (-3.0, -1.1045309498499094e-5),
5372 (-1.5, -0.017092677825984745),
5373 (-0.5, -0.27410803278438573),
5374 (0.0, -0.69314718055994531),
5375 (0.7, -1.8257336940742865),
5376 (2.0, -6.0580884451765829),
5377 (5.0, -27.89403672609738),
5378 (12.0, -147.75386135854695),
5379 ];
5380 for &(u, reference) in refs {
5381 let got = log_half_erfc_stable(u);
5382 let rel = (got - reference).abs() / reference.abs().max(1.0e-6);
5383 assert!(
5384 rel < 1.0e-12,
5385 "log_half_erfc_stable({u}) = {got:.17e}, reference {reference:.17e}, \
5386 rel {rel:.3e} >= 1e-12"
5387 );
5388 }
5389 }
5390
5391 pub(crate) fn cloglog_posterior_meanwith_deriv_gamma_reference(
5392 mu: f64,
5393 sigma: f64,
5394 ) -> Result<IntegratedMeanDerivative, EstimationError> {
5395 let survival = cloglog_survival_gamma_reference(mu, sigma)?;
5398 let shifted_survival = cloglog_survival_gamma_reference(mu + sigma * sigma, sigma)?;
5399 let mean = cloglog_mean_from_survival(survival);
5400 let dmean = cloglog_shift_identity_derivative(mu, sigma, shifted_survival);
5401 if !(mean.is_finite() && dmean.is_finite()) {
5402 crate::bail_invalid_estim!(
5403 "Gamma cloglog reference backend produced non-finite values"
5404 );
5405 }
5406 Ok(IntegratedMeanDerivative {
5407 mean,
5408 dmean_dmu: dmean.max(0.0),
5409 mode: IntegratedExpectationMode::ExactSpecialFunction,
5410 })
5411 }
5412
5413 fn even_moment_exp_neg_x2(power: usize) -> f64 {
5414 assert!(power.is_multiple_of(2));
5415 let m = power / 2;
5416 let mut odd_double_factorial = 1.0_f64;
5417 for k in 0..m {
5418 odd_double_factorial *= (2 * k + 1) as f64;
5419 }
5420 odd_double_factorial * std::f64::consts::PI.sqrt() / 2.0_f64.powi(m as i32)
5421 }
5422
5423 fn normal_pdf(z: f64) -> f64 {
5424 (-(z * z) * 0.5).exp() / (2.0 * std::f64::consts::PI).sqrt()
5425 }
5426
5427 fn high_res_sigmoid_integral(eta: f64, se: f64) -> f64 {
5428 let a = -12.0_f64;
5430 let b = 12.0_f64;
5431 let n = 20_000usize; let h = (b - a) / n as f64;
5433
5434 let integrand = |z: f64| -> f64 { sigmoid(eta + se * z) * normal_pdf(z) };
5435
5436 let mut sum = integrand(a) + integrand(b);
5437 for i in 1..n {
5438 let x = a + (i as f64) * h;
5439 if i % 2 == 0 {
5440 sum += 2.0 * integrand(x);
5441 } else {
5442 sum += 4.0 * integrand(x);
5443 }
5444 }
5445 sum * h / 3.0
5446 }
5447
5448 #[test]
5449 fn test_computed_nodes_symmetric() {
5450 let ctx = QuadratureContext::new();
5452 let gh = ctx.gauss_hermite();
5453 for i in 0..N_POINTS / 2 {
5454 let j = N_POINTS - 1 - i;
5455 assert_relative_eq!(gh.nodes[i], -gh.nodes[j], epsilon = 1e-12);
5456 }
5457 assert_relative_eq!(gh.nodes[N_POINTS / 2], 0.0, epsilon = 1e-12);
5459 }
5460
5461 #[test]
5462 fn test_computedweights_symmetric() {
5463 let ctx = QuadratureContext::new();
5465 let gh = ctx.gauss_hermite();
5466 for i in 0..N_POINTS / 2 {
5467 let j = N_POINTS - 1 - i;
5468 assert_relative_eq!(gh.weights[i], gh.weights[j], epsilon = 1e-12);
5469 }
5470 }
5471
5472 #[test]
5473 fn testweights_sum_to_sqrt_pi() {
5474 let ctx = QuadratureContext::new();
5476 let gh = ctx.gauss_hermite();
5477 let sum: f64 = gh.weights.iter().sum();
5478 assert_relative_eq!(sum, std::f64::consts::PI.sqrt(), epsilon = 1e-10);
5479 }
5480
5481 #[test]
5482 fn test_clenshaw_curtisweights_are_symmetric_and_integrate_constants() {
5483 let rule = compute_clenshaw_curtis_n(33);
5484 let m = rule.weights.len() - 1;
5485 for j in 0..=m / 2 {
5486 assert_relative_eq!(rule.nodes[j], -rule.nodes[m - j], epsilon = 1e-14);
5487 assert_relative_eq!(rule.weights[j], rule.weights[m - j], epsilon = 1e-14);
5488 }
5489 let sum: f64 = rule.weights.iter().sum();
5490 assert_relative_eq!(sum, 2.0, epsilon = 1e-14, max_relative = 1e-14);
5491 }
5492
5493 #[test]
5494 fn test_cc_preference_prefers_moderate_central_case() {
5495 assert!(cloglog_should_prefer_cc(-0.2, 0.8, CLOGLOG_CC_TOL));
5496 }
5497
5498 #[test]
5499 fn test_cc_preference_prefers_moderately_large_case() {
5500 assert!(cloglog_should_prefer_cc(0.0, 2.0, CLOGLOG_CC_TOL));
5501 }
5502
5503 #[test]
5504 fn test_cc_preference_rejects_broad_case() {
5505 assert!(!cloglog_should_prefer_cc(0.0, 5.0, CLOGLOG_CC_TOL));
5506 }
5507
5508 #[test]
5509 fn test_matches_abramowitz_stegun_7_point_gauss_hermite_constants() {
5510 let known_nodes = [
5514 -2.651_961_356_835_233_4,
5515 -1.673_551_628_767_471_4,
5516 -0.816_287_882_858_964_7,
5517 0.0,
5518 0.816_287_882_858_964_7,
5519 1.673_551_628_767_471_4,
5520 2.651_961_356_835_233_4,
5521 ];
5522 let knownweights = [
5523 0.000_971_781_245_099_519_1,
5524 0.054_515_582_819_127_03,
5525 0.425_607_252_610_127_8,
5526 0.810_264_617_556_807_3,
5527 0.425_607_252_610_127_8,
5528 0.054_515_582_819_127_03,
5529 0.000_971_781_245_099_519_1,
5530 ];
5531
5532 let ctx = QuadratureContext::new();
5533 let gh = ctx.gauss_hermite();
5534 for i in 0..N_POINTS {
5535 assert_relative_eq!(gh.nodes[i], known_nodes[i], epsilon = 1e-12);
5536 assert_relative_eq!(gh.weights[i], knownweights[i], epsilon = 1e-12);
5537 }
5538 }
5539
5540 #[test]
5541 fn testzero_se_returns_mode() {
5542 let eta = 1.5;
5544 let se = 0.0;
5545 let ctx = QuadratureContext::new();
5546 let mean = logit_posterior_mean(&ctx, eta, se);
5547 let mode = sigmoid(eta);
5548 assert_relative_eq!(mean, mode, epsilon = 1e-10);
5549 }
5550
5551 #[test]
5552 fn test_symmetric_atzero() {
5553 let eta = 0.0;
5555 let se = 1.0;
5556 let ctx = QuadratureContext::new();
5557 let mean = logit_posterior_mean(&ctx, eta, se);
5558 assert_relative_eq!(mean, 0.5, epsilon = 0.01);
5560 }
5561
5562 #[test]
5563 fn test_shrinkage_at_extremes() {
5564 let eta = 3.0; let se = 1.0;
5567 let ctx = QuadratureContext::new();
5568 let mean = logit_posterior_mean(&ctx, eta, se);
5569 let mode = sigmoid(eta);
5570
5571 assert!(mean < mode, "Expected mean {} < mode {}", mean, mode);
5573 assert!(mean > 0.8, "Mean {} should still be high", mean);
5575 }
5576
5577 #[test]
5578 fn test_matches_monte_carlo() {
5579 let eta = 2.0;
5581 let se = 0.8;
5582
5583 let ctx = QuadratureContext::new();
5584 let quad_mean = logit_posterior_mean(&ctx, eta, se);
5585
5586 let n_samples = 100_000;
5588 let mut mc_sum = 0.0;
5589 let mut rng_state = 12345u64; for _ in 0..n_samples {
5591 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5593 let u1 = ((rng_state as f64) / (u64::MAX as f64)).max(1e-10); rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5595 let u2 = (rng_state as f64) / (u64::MAX as f64);
5596 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
5597 let eta_sample = eta + se * z;
5598 mc_sum += sigmoid(eta_sample);
5599 }
5600 let mc_mean = mc_sum / (n_samples as f64);
5601
5602 assert_relative_eq!(quad_mean, mc_mean, epsilon = 0.01);
5604 }
5605
5606 #[test]
5607 fn test_quadrature_integrates_x_squared() {
5608 let ctx = QuadratureContext::new();
5611 let gh = ctx.gauss_hermite();
5612 let mut sum = 0.0;
5613 for i in 0..N_POINTS {
5614 sum += gh.weights[i] * gh.nodes[i] * gh.nodes[i];
5615 }
5616 let expected = std::f64::consts::PI.sqrt() / 2.0;
5617 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5618 }
5619
5620 #[test]
5621 fn test_quadrature_integrates_x_fourth() {
5622 let ctx = QuadratureContext::new();
5625 let gh = ctx.gauss_hermite();
5626 let mut sum = 0.0;
5627 for i in 0..N_POINTS {
5628 let x = gh.nodes[i];
5629 sum += gh.weights[i] * x * x * x * x;
5630 }
5631 let expected = 3.0 * std::f64::consts::PI.sqrt() / 4.0;
5632 assert_relative_eq!(sum, expected, epsilon = 1e-10);
5633 }
5634
5635 #[test]
5636 fn test_moment_exactness_up_to_degree_13() {
5637 let ctx = QuadratureContext::new();
5638 let gh = ctx.gauss_hermite();
5639
5640 for degree in 0..=13usize {
5641 let approx: f64 = (0..N_POINTS)
5642 .map(|i| gh.weights[i] * gh.nodes[i].powi(degree as i32))
5643 .sum();
5644
5645 let expected = if degree % 2 == 1 {
5646 0.0
5647 } else {
5648 even_moment_exp_neg_x2(degree)
5649 };
5650
5651 let err = (approx - expected).abs();
5652 let rel_scale = approx.abs().max(expected.abs()).max(1.0);
5653 assert!(
5654 err <= 1e-10 || err / rel_scale <= 1e-10,
5655 "degree={} approx={} expected={} abs_err={}",
5656 degree,
5657 approx,
5658 expected,
5659 err
5660 );
5661 }
5662 }
5663
5664 #[test]
5665 fn test_integrated_sigmoid_matches_high_res_integral_random_pairs() {
5666 let ctx = QuadratureContext::new();
5667 let mut rng_state = 0x4d595df4d0f33173u64;
5668
5669 for _ in 0..20 {
5670 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5671 let u_eta = (rng_state as f64) / (u64::MAX as f64);
5672 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
5673 let u_se = (rng_state as f64) / (u64::MAX as f64);
5674
5675 let eta = -6.0 + 12.0 * u_eta;
5676 let se = 0.02 + 1.5 * u_se;
5677
5678 let ghq = logit_posterior_mean(&ctx, eta, se);
5679 let numeric = high_res_sigmoid_integral(eta, se);
5680 assert_relative_eq!(ghq, numeric, epsilon = 2e-3);
5681 }
5682 }
5683
5684 #[test]
5685 fn test_logit_posterior_derivative_remains_positive_in_positive_tail() {
5686 let eta = 20.0;
5687 let se = 0.0;
5688 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5689 .expect("logit posterior mean derivative should evaluate");
5690 assert!(dmu > 0.0);
5691 assert!(
5692 dmu < 1e-6,
5693 "positive-tail derivative should stay tiny but nonzero, got {dmu}"
5694 );
5695 }
5696
5697 #[test]
5698 fn test_logit_posterior_derivative_matches_central_difference() {
5699 let ctx = QuadratureContext::new();
5700 let eta = 1.7;
5701 let se = 0.9;
5702 let h = 1e-5;
5703
5704 let (_, dmu) = logit_posterior_meanwith_deriv(eta, se)
5705 .expect("logit posterior mean derivative should evaluate");
5706 let mu_plus = logit_posterior_mean(&ctx, eta + h, se);
5707 let mu_minus = logit_posterior_mean(&ctx, eta - h, se);
5708 let dmufd = (mu_plus - mu_minus) / (2.0 * h);
5709
5710 assert_eq!(dmu.signum(), dmufd.signum());
5711 assert_relative_eq!(dmu, dmufd, epsilon = 5e-6, max_relative = 2e-4);
5712 }
5713
5714 fn dense_sigmoid_normal_mean(mu: f64, sigma: f64) -> f64 {
5720 let a = -18.0_f64;
5721 let b = 18.0_f64;
5722 let n = 400_000usize; let h = (b - a) / n as f64;
5724 let integrand = |z: f64| -> f64 { sigmoid(mu + sigma * z) * normal_pdf(z) };
5725 let mut sum = integrand(a) + integrand(b);
5726 for i in 1..n {
5727 let z = a + (i as f64) * h;
5728 sum += if i % 2 == 0 { 2.0 } else { 4.0 } * integrand(z);
5729 }
5730 sum * h / 3.0
5731 }
5732
5733 #[test]
5734 fn test_logit_posterior_mean_exact_symmetry_identity() {
5735 let cases = [
5738 (-3.0, 0.5),
5739 (-1.2, 1.7),
5740 (0.0, 2.2),
5741 (2.3, 0.8),
5742 (3.0, 0.05),
5743 ];
5744 for (mu, sigma) in cases {
5745 let p = logit_posterior_mean_exact(mu, sigma);
5746 let q = logit_posterior_mean_exact(-mu, sigma);
5747 assert!(
5748 (p + q - 1.0).abs() < 1e-12,
5749 "symmetry broken at mu={mu} sigma={sigma}: p+q-1 = {:.3e}",
5750 p + q - 1.0
5751 );
5752 }
5753 }
5754
5755 #[test]
5756 fn test_logit_posterior_mean_exact_matches_high_res_integral() {
5757 let cases = [
5761 (-2.0, 0.4),
5762 (-0.7, 1.1),
5763 (0.8, 0.9),
5764 (2.4, 1.7),
5765 (3.0, 0.05),
5766 (3.0, 0.5),
5767 (-2.0, 2.0),
5768 (5.0, 3.0),
5769 ];
5770 for (mu, sigma) in cases {
5771 let exact = logit_posterior_mean_exact(mu, sigma);
5772 let numeric = dense_sigmoid_normal_mean(mu, sigma);
5773 assert!(
5774 (exact - numeric).abs() < 1e-10,
5775 "oracle ≠ dense reference at mu={mu} sigma={sigma}: \
5776 exact={exact:.13} ref={numeric:.13} err={:.3e}",
5777 (exact - numeric).abs()
5778 );
5779 }
5780 }
5781
5782 #[test]
5789 fn test_logit_posterior_mean_exact_no_truncation_bias_1459() {
5790 let table = [
5793 (1.0, 0.02),
5794 (1.0, 0.05),
5795 (1.0, 0.5),
5796 (1.0, 2.0),
5797 (3.0, 0.02),
5798 (3.0, 0.05),
5799 (3.0, 0.5),
5800 (3.0, 2.0),
5801 (-2.0, 0.02),
5802 (-2.0, 0.05),
5803 (-2.0, 0.5),
5804 (-2.0, 2.0),
5805 ];
5806 for (mu, sigma) in table {
5807 let exact = logit_posterior_mean_exact(mu, sigma);
5808 let reference = dense_sigmoid_normal_mean(mu, sigma);
5809 let err = (exact - reference).abs();
5810 assert!(
5811 err < 1e-10,
5812 "#1459 truncation bias resurfaced at mu={mu} sigma={sigma}: \
5813 err={err:.3e} (pre-fix bias here was ~{:.2e})",
5814 mu.abs() / (2.0 * std::f64::consts::PI.powi(2) * 4096.0)
5815 );
5816 }
5817
5818 let mu = 3.0;
5823 let errs: Vec<f64> = [0.05, 0.5, 2.0]
5824 .iter()
5825 .map(|&s| logit_posterior_mean_exact(mu, s) - dense_sigmoid_normal_mean(mu, s))
5826 .collect();
5827 for e in &errs {
5828 assert!(
5829 e.abs() < 1e-10,
5830 "residual {e:.3e} at mu=3 — old σ-independent plateau was 3.71e-5"
5831 );
5832 }
5833 }
5834
5835 #[test]
5843 fn test_faddeeva_weideman_matches_known_values() {
5844 let w0 = faddeeva_upper_halfplane(Complex { re: 0.0, im: 0.0 });
5846 assert!(
5847 (w0.re - 1.0).abs() < 1e-13 && w0.im.abs() < 1e-13,
5848 "w(0)={w0:?}"
5849 );
5850 let on_axis = [
5852 (0.1, 0.8964569799691268),
5853 (0.5, 0.6156903441929258),
5854 (1.0, 0.427583576155807),
5855 (2.0, 0.2553956763105058),
5856 (5.0, 0.11070463773306861),
5857 (9.0, 0.06230772403777468),
5858 ];
5859 for (y, want) in on_axis {
5860 let w = faddeeva_upper_halfplane(Complex { re: 0.0, im: y });
5861 assert!(
5862 (w.re - want).abs() < 1e-13 && w.im.abs() < 1e-13,
5863 "w(i·{y}): got {w:?}, want re={want}, err={:.2e}",
5864 (w.re - want).abs()
5865 );
5866 }
5867 let off_axis = [
5869 ((0.7, 1.3), (0.31327301971562715, 0.12443489420104513)),
5870 ((-1.5, 0.8), (0.21066359024766423, -0.27001624496296617)),
5871 ((3.0, 0.4), (0.030278754646989155, 0.1957320888774461)),
5872 ];
5873 for ((re, im), (wre, wim)) in off_axis {
5874 let w = faddeeva_upper_halfplane(Complex { re, im });
5875 assert!(
5876 (w.re - wre).abs() < 1e-13 && (w.im - wim).abs() < 1e-13,
5877 "w({re}+{im}i): got {w:?}, want ({wre},{wim})"
5878 );
5879 }
5880 let w = faddeeva_upper_halfplane(Complex { re: 3.0, im: 40.0 });
5884 assert!(
5885 (w.re - 0.01402158696172506).abs() < 1e-13
5886 && (w.im - 0.0010509664408184546).abs() < 1e-13,
5887 "tail value mismatch: w={w:?}"
5888 );
5889 }
5890
5891 #[test]
5892 fn test_integrated_logit_mean_close_to_exact_oracle() {
5893 let ctx = QuadratureContext::new();
5897 let cases = [(-3.0, 0.3), (-1.0, 0.8), (0.5, 1.2), (2.8, 1.0)];
5898 for (eta, se) in cases {
5899 let ghq = logit_posterior_mean(&ctx, eta, se);
5900 let exact = logit_posterior_mean_exact(eta, se);
5901 assert!(
5902 (ghq - exact).abs() < 1e-6,
5903 "production path drifts from oracle at eta={eta} se={se}: \
5904 ghq={ghq:.12} oracle={exact:.12} gap={:.3e}",
5905 (ghq - exact).abs()
5906 );
5907 }
5908 }
5909
5910 #[test]
5911 fn test_probit_posterior_mean_reduces_to_map_atzero_se() {
5912 let eta = 1.25;
5913 let p = probit_posterior_mean(eta, 0.0);
5914 let map = gam_math::probability::normal_cdf(eta);
5915 assert_relative_eq!(p, map, epsilon = 1e-12);
5916 }
5917
5918 #[test]
5919 fn test_probit_posterior_mean_shrinks_extremeswith_uncertainty() {
5920 let hi_eta = 3.0;
5921 let lo_eta = -3.0;
5922 let p_hi_map = probit_posterior_mean(hi_eta, 0.0);
5923 let p_hi_unc = probit_posterior_mean(hi_eta, 2.0);
5924 let p_lo_map = probit_posterior_mean(lo_eta, 0.0);
5925 let p_lo_unc = probit_posterior_mean(lo_eta, 2.0);
5926 assert!(p_hi_unc < p_hi_map);
5927 assert!(p_lo_unc > p_lo_map);
5928 }
5929
5930 #[test]
5931 fn test_survival_posterior_mean_is_bounded_and_shrinks_tail() {
5932 let ctx = QuadratureContext::new();
5933 let eta: f64 = 3.0;
5934 let map = (-(eta.exp())).exp();
5935 let pm = survival_posterior_mean(&ctx, eta, 1.5);
5936 assert!((0.0..=1.0).contains(&pm));
5937 assert!(pm > map);
5938 }
5939
5940 #[test]
5941 fn test_cloglog_and_survival_posterior_means_are_complements() {
5942 let ctx = QuadratureContext::new();
5943 let cases = [
5944 (-3.0, 0.0),
5945 (-0.2, 0.1),
5946 (0.4, 0.8),
5947 (2.0, 1.5),
5948 (10.0, 0.3),
5949 (0.0, 20.0),
5950 (10.0, 10.0),
5951 (-0.5, 100.0),
5952 ];
5953 for (eta, se) in cases {
5954 let clog = cloglog_posterior_mean(&ctx, eta, se);
5955 let surv = survival_posterior_mean(&ctx, eta, se);
5956 assert_relative_eq!(clog + surv, 1.0, epsilon = 2e-10, max_relative = 2e-10);
5957 }
5958 }
5959
5960 #[test]
5961 fn test_cloglog_and_survival_share_large_sigmaspecial_function_path() {
5962 let ctx = QuadratureContext::new();
5963 let eta = -0.2;
5964 let se = 0.8;
5965 let clog = cloglog_posterior_mean(&ctx, eta, se);
5966 let surv = survival_posterior_mean(&ctx, eta, se);
5967 let integrated =
5968 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, eta, se)
5969 .expect("cloglog integrated inverse-link moments should evaluate");
5970 assert_eq!(
5971 integrated.mode,
5972 IntegratedExpectationMode::ExactSpecialFunction
5973 );
5974 assert_relative_eq!(clog, integrated.mean, epsilon = 1e-12, max_relative = 1e-12);
5975 assert_relative_eq!(clog + surv, 1.0, epsilon = 1e-10, max_relative = 1e-10);
5976 }
5977
5978 #[test]
5979 fn test_cloglog_and_survival_posteriorvariances_match() {
5980 let ctx = QuadratureContext::new();
5981 let cases = [(-3.0, 0.0), (-0.2, 0.1), (0.4, 0.8), (2.0, 1.5)];
5982 for (eta, se) in cases {
5983 let (_, clogvar) = cloglog_posterior_meanvariance(&ctx, eta, se);
5984 let (_, survvar) = survival_posterior_meanvariance(&ctx, eta, se);
5985 assert_relative_eq!(clogvar, survvar, epsilon = 1e-12, max_relative = 1e-12);
5986 }
5987 }
5988
5989 #[test]
5990 fn test_survivalvariance_uses_exactsecond_moment_shift() {
5991 let ctx = QuadratureContext::new();
5992 let eta = -0.2;
5993 let se = 0.8;
5994 let (survival, _) = cloglog_survival_term_controlled(&ctx, eta, se);
5995 let (survival_sq, _) = cloglog_survivalsecond_moment_controlled(&ctx, eta, se);
5996 let (_, variance) = survival_posterior_meanvariance(&ctx, eta, se);
5997 assert_relative_eq!(
5998 variance,
5999 (survival_sq - survival * survival).max(0.0),
6000 epsilon = 1e-12,
6001 max_relative = 1e-12
6002 );
6003 }
6004
6005 #[test]
6006 fn test_lognormal_laplace_shift_matches_explicitmu_plus_logz() {
6007 let ctx = QuadratureContext::new();
6008 let mu = -0.2;
6009 let sigma = 0.8;
6010 let z = 2.0;
6011 let shifted = lognormal_laplace_term_controlled(&ctx, z, mu, sigma);
6012 let explicit = cloglog_survival_term_controlled(&ctx, mu + z.ln(), sigma);
6013 assert_eq!(shifted.1, explicit.1);
6014 assert_relative_eq!(shifted.0, explicit.0, epsilon = 1e-12, max_relative = 1e-12);
6015 }
6016
6017 #[test]
6018 fn test_integrated_dispatch_uses_closed_form_probit() {
6019 let ctx = QuadratureContext::new();
6020 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Probit, 0.7, 1.3)
6021 .expect("probit integrated inverse-link moments should evaluate");
6022 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
6023 let direct = probit_posterior_meanwith_deriv_exact(0.7, 1.3);
6024 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
6025 assert_relative_eq!(out.dmean_dmu, direct.dmean_dmu, epsilon = 1e-12);
6026 }
6027
6028 #[test]
6029 fn test_integrated_probit_jet_matches_closed_form_derivatives() {
6030 let ctx = QuadratureContext::new();
6031 let mu = 0.7;
6032 let sigma = 1.3;
6033 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Probit, mu, sigma)
6034 .expect("probit integrated inverse-link jet should evaluate");
6035 let s = (1.0 + sigma * sigma).sqrt();
6036 let z = mu / s;
6037 let pdf = gam_math::probability::normal_pdf(z);
6038 assert_relative_eq!(
6039 out.mean,
6040 gam_math::probability::normal_cdf(z),
6041 epsilon = 1e-12
6042 );
6043 assert_relative_eq!(out.d1, pdf / s, epsilon = 1e-12);
6044 assert_relative_eq!(out.d2, -z * pdf / (s * s), epsilon = 1e-12);
6045 assert_relative_eq!(out.d3, (z * z - 1.0) * pdf / (s * s * s), epsilon = 1e-12);
6046 }
6047
6048 #[test]
6049 fn test_integrated_logit_jet_matches_central_differences() {
6050 let ctx = QuadratureContext::new();
6063 let mu = 1.1;
6064 let sigma = 0.8;
6065 let out = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6066 .expect("logit integrated inverse-link jet should evaluate");
6067 assert!(matches!(
6068 out.mode,
6069 IntegratedExpectationMode::ExactSpecialFunction
6070 | IntegratedExpectationMode::QuadratureFallback
6071 ));
6072 let (ref_mean, ref_d1, ref_d2, ref_d3) = logit_reference_jet_highres_simpson(mu, sigma);
6073 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6074 assert_relative_eq!(out.d1, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6075 assert_relative_eq!(out.d2, ref_d2, epsilon = 1e-11, max_relative = 1e-10);
6076 assert_relative_eq!(out.d3, ref_d3, epsilon = 1e-11, max_relative = 1e-10);
6077 }
6078
6079 #[test]
6080 fn test_integrated_cloglog_jet_matches_central_differences() {
6081 let ctx = QuadratureContext::new();
6082 let mu = 0.4;
6083 let sigma = 0.6;
6084 let h = 1e-4;
6085 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
6086 .expect("cloglog integrated inverse-link jet should evaluate");
6087 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
6088 .expect("cloglog integrated inverse-link jet should evaluate");
6089 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
6090 .expect("cloglog integrated inverse-link jet should evaluate");
6091 let d1fd = (plus.mean - minus.mean) / (2.0 * h);
6092 let d2fd = (plus.d1 - minus.d1) / (2.0 * h);
6093 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
6094 assert_eq!(out.d1.signum(), d1fd.signum());
6095 assert_eq!(out.d2.signum(), d2fd.signum());
6096 assert_eq!(out.d3.signum(), d3fd.signum());
6097 assert_relative_eq!(out.d1, d1fd, epsilon = 2e-5, max_relative = 3e-4);
6098 assert_relative_eq!(out.d2, d2fd, epsilon = 4e-5, max_relative = 8e-4);
6099 assert_relative_eq!(out.d3, d3fd, epsilon = 8e-5, max_relative = 2e-3);
6100 }
6101
6102 #[test]
6103 fn test_integrated_cloglog_wide_sigma_d3_matches_simpson_and_d2_slope() {
6104 let ctx = QuadratureContext::new();
6105 let cases = [(0.0, 4.0), (-1.0, 4.0), (2.0, 3.0), (3.0, 3.0)];
6106 let h = 1e-4;
6107
6108 for (mu, sigma) in cases {
6109 let out = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu, sigma)
6110 .expect("wide-sigma cloglog integrated jet should evaluate");
6111 let reference = cloglog_reference_jet_highres_simpson(mu, sigma);
6112 let plus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu + h, sigma)
6113 .expect("wide-sigma cloglog integrated jet should evaluate");
6114 let minus = integrated_inverse_link_jet(&ctx, LinkFunction::CLogLog, mu - h, sigma)
6115 .expect("wide-sigma cloglog integrated jet should evaluate");
6116 let d3fd = (plus.d2 - minus.d2) / (2.0 * h);
6117
6118 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
6119 assert_relative_eq!(out.mean, reference.0, epsilon = 4e-8, max_relative = 4e-8);
6120 assert_relative_eq!(out.d1, reference.1, epsilon = 4e-8, max_relative = 4e-8);
6121 assert_relative_eq!(out.d2, reference.2, epsilon = 2e-9, max_relative = 2e-7);
6122 assert_relative_eq!(out.d3, reference.3, epsilon = 2e-9, max_relative = 2e-7);
6123 assert_relative_eq!(out.d3, d3fd, epsilon = 2e-7, max_relative = 4e-5);
6124 }
6125 }
6126
6127 #[test]
6128 fn test_latent_cloglog_jet5_matches_higher_order_central_differences() {
6129 let ctx = QuadratureContext::new();
6130 let mu = 0.35;
6131 let sigma = 0.7;
6132 let h = 2e-4;
6133
6134 let out = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu, sigma);
6135 let plus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu + h, sigma);
6136 let minus = latent_cloglog_inverse_link_jet5_controlled(&ctx, mu - h, sigma);
6137
6138 let d4fd = (plus.d3 - minus.d3) / (2.0 * h);
6139 let d5fd = (plus.d4 - minus.d4) / (2.0 * h);
6140
6141 assert_eq!(out.d4.signum(), d4fd.signum());
6142 assert_eq!(out.d5.signum(), d5fd.signum());
6143 assert_relative_eq!(out.d4, d4fd, epsilon = 2e-4, max_relative = 5e-3);
6144 assert_relative_eq!(out.d5, d5fd, epsilon = 6e-4, max_relative = 2e-2);
6145 }
6146
6147 #[test]
6148 fn test_logit_exact_derivative_matches_finite_difference() {
6149 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("controlled logit");
6159 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6160 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6161 assert!(out.dmean_dmu > 0.0);
6162 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6163 }
6164
6165 #[test]
6166 fn test_logit_small_sigma_returns_quadrature_truth_2623() {
6167 for &(mu, sigma) in &[
6173 (-3.0, 0.05),
6174 (-0.5, 0.10),
6175 (0.0, 0.20),
6176 (0.5, 0.24),
6177 (3.0, 0.15),
6178 ] {
6179 let out =
6180 logit_posterior_meanwith_deriv_controlled(mu, sigma).expect("controlled logit");
6181 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6182 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
6183 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-12, max_relative = 1e-11);
6184 assert_relative_eq!(
6185 out.dmean_dmu,
6186 ref_d1,
6187 epsilon = 1e-12,
6188 max_relative = 1e-11
6189 );
6190 }
6191 }
6192
6193 #[test]
6194 fn test_logit_exact_clamped_degenerate_branch_is_locally_flat() {
6195 let out = logit_posterior_meanwith_deriv_exact(-710.0, 0.0).expect("exact logit");
6196 let h = 1e-6;
6197 let plus = logit_posterior_meanwith_deriv_exact(-710.0 + h, 0.0)
6198 .expect("exact logit plus")
6199 .mean;
6200 let minus = logit_posterior_meanwith_deriv_exact(-710.0 - h, 0.0)
6201 .expect("exact logit minus")
6202 .mean;
6203 let fd = (plus - minus) / (2.0 * h);
6204 assert_eq!(fd, 0.0);
6205 assert_eq!(out.dmean_dmu, 0.0);
6206 }
6207
6208 fn simpson_integrate<F>(a: f64, b: f64, n_intervals: usize, f: F) -> f64
6209 where
6210 F: Fn(f64) -> f64,
6211 {
6212 assert_eq!(n_intervals % 2, 0, "Simpson integration requires an even n");
6213 let h = (b - a) / n_intervals as f64;
6214 let mut sum = f(a) + f(b);
6215 for i in 1..n_intervals {
6216 let x = a + i as f64 * h;
6217 let w = if i % 2 == 0 { 2.0 } else { 4.0 };
6218 sum += w * f(x);
6219 }
6220 sum * h / 3.0
6221 }
6222
6223 fn cloglog_reference_mean_and_derivative(mu: f64, sigma: f64) -> (f64, f64) {
6224 if sigma <= CLOGLOG_SIGMA_DEGENERATE {
6225 return (cloglog_mean_exact(mu), cloglog_mean_d1_exact(mu));
6226 }
6227
6228 let z_max = 12.0;
6232 let n_intervals = 4096;
6233 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
6234 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6235 let eta = mu + sigma * z;
6236 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_exact(eta)
6237 });
6238 let deriv = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6239 let eta = mu + sigma * z;
6240 inv_sqrt_2pi * (-0.5 * z * z).exp() * cloglog_mean_d1_exact(eta)
6241 });
6242 (mean, deriv)
6243 }
6244
6245 fn logit_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
6258 let z_max = 14.0;
6259 let n_intervals = 16384;
6260 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
6261 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
6262 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6263 let eta = mu + sigma * z;
6264 let (p, _, _, _) = component_point_jet(LinkComponent::Logit, eta);
6265 phi(z) * p
6266 });
6267 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6268 let eta = mu + sigma * z;
6269 let (_, p1, _, _) = component_point_jet(LinkComponent::Logit, eta);
6270 phi(z) * p1
6271 });
6272 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6273 let eta = mu + sigma * z;
6274 let (_, _, p2, _) = component_point_jet(LinkComponent::Logit, eta);
6275 phi(z) * p2
6276 });
6277 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6278 let eta = mu + sigma * z;
6279 let (_, _, _, p3) = component_point_jet(LinkComponent::Logit, eta);
6280 phi(z) * p3
6281 });
6282 (mean, d1, d2, d3)
6283 }
6284
6285 fn cloglog_reference_jet_highres_simpson(mu: f64, sigma: f64) -> (f64, f64, f64, f64) {
6286 let z_max = 14.0;
6287 let n_intervals = 16384;
6288 let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
6289 let phi = |z: f64| inv_sqrt_2pi * (-0.5 * z * z).exp();
6290 let mean = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6291 let eta = mu + sigma * z;
6292 let (g, _, _, _, _, _) = cloglog_point_jet5(eta);
6293 phi(z) * g
6294 });
6295 let d1 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6296 let eta = mu + sigma * z;
6297 let (_, g1, _, _, _, _) = cloglog_point_jet5(eta);
6298 phi(z) * g1
6299 });
6300 let d2 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6301 let eta = mu + sigma * z;
6302 let (_, _, g2, _, _, _) = cloglog_point_jet5(eta);
6303 phi(z) * g2
6304 });
6305 let d3 = simpson_integrate(-z_max, z_max, n_intervals, |z| {
6306 let eta = mu + sigma * z;
6307 let (_, _, _, g3, _, _) = cloglog_point_jet5(eta);
6308 phi(z) * g3
6309 });
6310 (mean, d1, d2, d3)
6311 }
6312
6313 #[test]
6314 fn test_cloglog_taylor_negative_tail_matches_mathematical_target() {
6315 let mu = -40.0;
6316 let sigma = 0.1;
6317 let out = cloglog_small_sigma_taylor(mu, sigma);
6318 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
6319
6320 assert!(
6321 out.dmean_dmu > 0.0,
6322 "negative-tail derivative should remain positive"
6323 );
6324 assert_relative_eq!(
6325 out.mean,
6326 expected_mean,
6327 epsilon = 1e-30,
6328 max_relative = 1e-12
6329 );
6330 assert_relative_eq!(
6331 out.dmean_dmu,
6332 expected_deriv,
6333 epsilon = 1e-30,
6334 max_relative = 1e-12
6335 );
6336 }
6337
6338 #[test]
6339 fn test_cloglog_degenerate_negative_tail_matches_pointwise_target() {
6340 let ctx = QuadratureContext::new();
6341 let mu = -40.0;
6342 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, 0.0);
6343
6344 assert!(
6345 out.dmean_dmu > 0.0,
6346 "degenerate negative-tail derivative should remain positive"
6347 );
6348 assert_relative_eq!(
6349 out.mean,
6350 cloglog_mean_exact(mu),
6351 epsilon = 1e-30,
6352 max_relative = 1e-15
6353 );
6354 assert_relative_eq!(
6355 out.dmean_dmu,
6356 cloglog_mean_d1_exact(mu),
6357 epsilon = 1e-30,
6358 max_relative = 1e-15
6359 );
6360 }
6361
6362 #[test]
6363 fn test_degenerate_probit_jet_is_exact_beyond_former_clamp() {
6364 let mu = -30.1;
6365 let probit = integrated_probit_jet(mu, 0.0);
6366 let pdf = gam_math::probability::normal_pdf(mu);
6367 assert!(
6368 pdf > 0.0,
6369 "test point must have a represented Gaussian tail"
6370 );
6371 assert_eq!(probit.mean, gam_math::probability::normal_cdf(mu));
6372 assert_eq!(probit.d1, pdf);
6373 assert_eq!(probit.d2, -mu * pdf);
6374 assert_eq!(probit.d3, (mu * mu - 1.0) * pdf);
6375
6376 let tail = (-710.0_f64).exp();
6392 assert!(
6393 tail > 0.0 && tail < f64::MIN_POSITIVE,
6394 "eta=-710 must sit in the subnormal tail, not underflow"
6395 );
6396 let logit = component_point_jet(LinkComponent::Logit, -710.0);
6397 assert_eq!(logit.0, tail);
6398 assert_eq!(logit.1, tail);
6399 assert_eq!(logit.2, tail);
6400 assert_eq!(logit.3, tail);
6401
6402 assert_eq!(
6404 (-750.0_f64).exp(),
6405 0.0,
6406 "eta=-750 must underflow f64 for this arm to mean anything"
6407 );
6408 let underflowed = component_point_jet(LinkComponent::Logit, -750.0);
6409 assert_eq!(underflowed.1, 0.0);
6410 assert_eq!(underflowed.2, 0.0);
6411 assert_eq!(underflowed.3, 0.0);
6412 }
6413
6414 #[test]
6415 fn test_degenerate_cloglog_component_jet_preserves_smooth_negative_tail() {
6416 let eta: f64 = -40.0;
6417 let t = eta.exp();
6418 let s = (-t).exp();
6419 let cloglog = component_point_jet(LinkComponent::CLogLog, eta);
6420 let expected_mean = -(-t).exp_m1();
6421 let expected_d1 = t * s;
6422 let expected_d2 = (t - t * t) * s;
6423 let expected_d3 = (t - 3.0 * t * t + t * t * t) * s;
6424
6425 assert!(cloglog.1 > 0.0, "negative-tail d1 should remain positive");
6426 assert_relative_eq!(
6427 cloglog.0,
6428 expected_mean,
6429 epsilon = 1e-30,
6430 max_relative = 1e-15
6431 );
6432 assert_relative_eq!(
6433 cloglog.1,
6434 expected_d1,
6435 epsilon = 1e-30,
6436 max_relative = 1e-15
6437 );
6438 assert_relative_eq!(
6439 cloglog.2,
6440 expected_d2,
6441 epsilon = 1e-30,
6442 max_relative = 1e-15
6443 );
6444 assert_relative_eq!(
6445 cloglog.3,
6446 expected_d3,
6447 epsilon = 1e-30,
6448 max_relative = 1e-15
6449 );
6450 }
6451
6452 #[test]
6453 fn test_zero_sigma_logit_and_cloglog_share_component_tail_jets() {
6454 let ctx = QuadratureContext::new();
6455 for (link, component, eta) in [
6456 (LinkFunction::Logit, LinkComponent::Logit, 50.0),
6457 (LinkFunction::CLogLog, LinkComponent::CLogLog, -50.0),
6458 ] {
6459 let integrated = integrated_inverse_link_jet(&ctx, link, eta, 0.0)
6460 .expect("degenerate integrated jet");
6461 let point = component_inverse_link_jet(component, eta);
6462 assert_eq!(integrated.mode, IntegratedExpectationMode::ExactClosedForm);
6463 assert_eq!(integrated.mean, point.mu);
6464 assert_eq!(integrated.d1, point.d1);
6465 assert_eq!(integrated.d2, point.d2);
6466 assert_eq!(integrated.d3, point.d3);
6467 }
6468 }
6469
6470 #[test]
6471 fn test_cloglog_controlled_matches_mathematical_target_on_small_sigma_grid() {
6472 let ctx = QuadratureContext::new();
6473 let cases = [
6477 (-30.0, 1e-10),
6478 (-30.0, 0.1),
6479 (-10.0, 0.24),
6480 (-3.0, 0.2),
6481 (0.0, 0.05),
6482 (0.4, 0.1),
6483 (3.0, 0.24),
6484 (10.0, 0.1),
6485 (30.0, 0.24),
6486 ];
6487
6488 for &(mu, sigma) in &cases {
6489 let approx = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
6490 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
6491 assert_relative_eq!(
6492 approx.mean,
6493 expected_mean,
6494 epsilon = 1e-12,
6495 max_relative = 2e-3
6496 );
6497 assert_relative_eq!(
6498 approx.dmean_dmu,
6499 expected_deriv,
6500 epsilon = 1e-12,
6501 max_relative = 4e-3
6502 );
6503 }
6504 }
6505
6506 #[test]
6507 fn test_cloglog_dispatch_uses_gamma_backend_for_large_sigma_central_regime() {
6508 let ctx = QuadratureContext::new();
6509 let out =
6510 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, -0.2, 0.8)
6511 .expect("cloglog integrated inverse-link moments should evaluate");
6512 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6513 assert!(out.mean.is_finite());
6514 assert!(out.dmean_dmu.is_finite());
6515 assert!(out.dmean_dmu >= 0.0);
6516 }
6517
6518 #[test]
6519 fn test_cloglog_dispatch_uses_large_sigma_asymptotic_without_ghq() {
6520 let ctx = QuadratureContext::new();
6521 let out =
6522 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::CLogLog, 0.0, 20.0)
6523 .expect("cloglog integrated inverse-link moments should evaluate");
6524 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6525 assert!(out.mean.is_finite());
6526 assert!(out.dmean_dmu.is_finite());
6527 assert!(out.dmean_dmu >= 0.0);
6528 }
6529
6530 #[test]
6531 fn test_cloglog_cc_matches_gamma_reference_on_central_case() {
6532 let ctx = QuadratureContext::new();
6533 let mu = -0.2;
6534 let sigma = 0.8;
6535 let cc = cloglog_survival_cc(&ctx, mu, sigma, CLOGLOG_CC_TOL).expect("cc backend");
6536 let gamma = cloglog_survival_gamma_reference(mu, sigma).expect("gamma backend");
6537 assert_relative_eq!(cc, gamma, epsilon = 5e-6, max_relative = 5e-6);
6538 }
6539
6540 #[test]
6541 fn test_cloglog_gamma_reference_matches_seeded_monte_carlo_small_case() {
6542 let mu = -0.2;
6543 let sigma = 0.8;
6544 let gamma =
6545 cloglog_posterior_meanwith_deriv_gamma_reference(mu, sigma).expect("gamma reference");
6546 let mut rng_state = 0x9e3779b97f4a7c15u64;
6547 let mut mean_mc = 0.0f64;
6548 let mut deriv_mc = 0.0f64;
6549 let n_samples = 300_000usize;
6550 for _ in 0..n_samples {
6551 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6552 let u1 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6553 rng_state = rng_state.wrapping_mul(6364136223846793005).wrapping_add(1);
6554 let u2 = ((rng_state as f64) / (u64::MAX as f64)).clamp(1e-12, 1.0 - 1e-12);
6555 let z = (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
6556 let eta = mu + sigma * z;
6557 mean_mc += cloglog_mean_exact(eta);
6558 deriv_mc += cloglog_mean_d1_exact(eta);
6559 }
6560 mean_mc /= n_samples as f64;
6561 deriv_mc /= n_samples as f64;
6562 assert_relative_eq!(gamma.mean, mean_mc, epsilon = 2e-3, max_relative = 2e-3);
6563 assert_relative_eq!(
6564 gamma.dmean_dmu,
6565 deriv_mc,
6566 epsilon = 2e-3,
6567 max_relative = 2e-3
6568 );
6569 }
6570
6571 #[test]
6572 fn test_logit_dispatch_uses_tail_asymptotic_outside_old_guard() {
6573 let ctx = QuadratureContext::new();
6574 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 35.0, 1.0)
6575 .expect("logit integrated inverse-link moments should evaluate");
6576 assert_eq!(out.mode, IntegratedExpectationMode::ControlledAsymptotic);
6577 assert!(out.mean.is_finite());
6578 assert!(out.dmean_dmu.is_finite());
6579 assert!(out.dmean_dmu >= 0.0);
6580 }
6581
6582 #[test]
6583 fn test_logit_dispatch_prefers_erfcx_in_moderate_regime() {
6584 let ctx = QuadratureContext::new();
6595 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 1.1, 0.8)
6596 .expect("logit integrated inverse-link moments should evaluate");
6597 assert!(matches!(
6598 out.mode,
6599 IntegratedExpectationMode::ExactSpecialFunction
6600 | IntegratedExpectationMode::QuadratureFallback
6601 ));
6602 assert!(out.mean.is_finite());
6603 assert!(out.dmean_dmu.is_finite());
6604 assert!(out.dmean_dmu >= 0.0);
6605 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6606 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6607 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6608 }
6609
6610 #[test]
6611 fn test_logit_dispatch_large_sigma_uses_accurate_quadrature_not_monahan() {
6612 let ctx = QuadratureContext::new();
6621 let out = integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, 0.5, 20.0)
6622 .expect("logit integrated inverse-link moments should evaluate");
6623 assert_eq!(out.mode, IntegratedExpectationMode::QuadratureFallback);
6624 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(0.5, 20.0);
6625 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6626 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6627 let kappa = (1.0 + std::f64::consts::PI * 20.0 * 20.0 / 8.0)
6630 .sqrt()
6631 .recip();
6632 let monahan_mean = gam_math::probability::normal_cdf(0.5 * kappa);
6633 assert!(
6634 (out.mean - monahan_mean).abs() > 1e-3,
6635 "dispatcher must not return the inaccurate Monahan mean {monahan_mean}; got {}",
6636 out.mean
6637 );
6638 }
6639
6640 #[test]
6641 fn test_logit_controlled_path_keeps_exact_backend_in_moderate_regime() {
6642 let out = logit_posterior_meanwith_deriv_controlled(1.1, 0.8).expect("logit controlled");
6652 assert!(matches!(
6653 out.mode,
6654 IntegratedExpectationMode::ExactSpecialFunction
6655 | IntegratedExpectationMode::QuadratureFallback
6656 ));
6657 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(1.1, 0.8);
6658 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-11, max_relative = 1e-10);
6659 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-11, max_relative = 1e-10);
6660 }
6661
6662 #[test]
6663 fn test_logit_dispatch_derivative_correct_at_mu_zero_small_sigma() {
6664 let ctx = QuadratureContext::new();
6673 for &(mu, sigma) in &[(0.0, 0.3), (0.0, 0.4), (0.0, 0.5)] {
6674 let out =
6675 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6676 .expect("logit integrated inverse-link moments should evaluate");
6677 assert_relative_eq!(out.mean, 0.5, epsilon = 1e-10);
6679 assert!(
6681 out.dmean_dmu <= 0.25 + 1e-9,
6682 "E[sigmoid'] must not exceed 0.25 at (μ={mu}, σ={sigma}); got {}",
6683 out.dmean_dmu
6684 );
6685 let (_, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6686 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-6);
6687 }
6688 }
6689
6690 #[test]
6691 fn test_logit_erfcx_exact_branch_is_self_certified() {
6692 for &(mu, sigma) in &[(8.0, 1.0), (10.0, 1.0), (15.0, 2.0)] {
6699 let out = logit_posterior_meanwith_deriv_exact(mu, sigma)
6700 .expect("erfcx branch should certify");
6701 assert_eq!(out.mode, IntegratedExpectationMode::ExactSpecialFunction);
6702 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6703 assert_relative_eq!(out.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-7);
6704 assert_relative_eq!(out.dmean_dmu, ref_d1, epsilon = 1e-9, max_relative = 1e-7);
6705 }
6706 assert!(
6710 logit_posterior_meanwith_deriv_exact(0.0, 0.3).is_err(),
6711 "erfcx branch must not claim ExactSpecialFunction when it cannot certify the derivative"
6712 );
6713 }
6714
6715 #[test]
6716 fn test_logit_integrated_derivative_is_even_in_mu() {
6717 let ctx = QuadratureContext::new();
6723 for &(mu, sigma) in &[(0.3, 0.3), (1.1, 0.8), (10.0, 1.0), (3.0, 3.0), (35.0, 1.0)] {
6724 let pos =
6725 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6726 .expect("logit moments (+μ)");
6727 let neg =
6728 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, -mu, sigma)
6729 .expect("logit moments (-μ)");
6730 assert_relative_eq!(
6731 pos.dmean_dmu,
6732 neg.dmean_dmu,
6733 epsilon = 1e-9,
6734 max_relative = 1e-7
6735 );
6736 assert_relative_eq!(
6738 neg.mean,
6739 1.0 - pos.mean,
6740 epsilon = 1e-9,
6741 max_relative = 1e-7
6742 );
6743 }
6744 }
6745
6746 #[test]
6747 fn test_logit_dmean_dmu_equals_fd_of_mean_across_regimes() {
6748 let ctx = QuadratureContext::new();
6763 let h = 1e-4;
6764 let cases = [
6765 (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), ];
6775 for &(mu, sigma) in &cases {
6776 let at = |m: f64| {
6777 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, m, sigma)
6778 .expect("logit moments")
6779 };
6780 let out = at(mu);
6781 let fd = (at(mu + h).mean - at(mu - h).mean) / (2.0 * h);
6782 assert!(
6783 (out.dmean_dmu - fd).abs() <= 1e-5,
6784 "dmean_dmu must equal d/dμ of mean at (μ={mu}, σ={sigma}): \
6785 returned {}, FD of mean {} (mode {:?})",
6786 out.dmean_dmu,
6787 fd,
6788 out.mode
6789 );
6790 assert!(
6794 out.dmean_dmu <= 0.25 + 1e-9 && out.dmean_dmu >= 0.0,
6795 "dmean_dmu out of [0, 0.25] at (μ={mu}, σ={sigma}): {}",
6796 out.dmean_dmu
6797 );
6798 }
6799 }
6800
6801 #[test]
6802 fn test_logit_scalar_matches_jet_at_large_sigma() {
6803 let ctx = QuadratureContext::new();
6809 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0)] {
6810 let scalar =
6811 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6812 .expect("scalar logit moments");
6813 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6814 .expect("jet logit moments");
6815 let (ref_mean, ref_d1, _, _) = logit_reference_jet_highres_simpson(mu, sigma);
6819 assert_relative_eq!(scalar.mean, ref_mean, epsilon = 1e-9, max_relative = 1e-8);
6820 assert_relative_eq!(
6821 scalar.dmean_dmu,
6822 ref_d1,
6823 epsilon = 1e-9,
6824 max_relative = 1e-8
6825 );
6826 assert_relative_eq!(scalar.mean, jet.mean, epsilon = 1e-12, max_relative = 1e-12);
6834 assert_relative_eq!(
6835 scalar.dmean_dmu,
6836 jet.d1,
6837 epsilon = 1e-12,
6838 max_relative = 1e-12
6839 );
6840 }
6841 }
6842
6843 #[test]
6844 fn test_logit_jet_accurate_at_wide_sigma() {
6845 let ctx = QuadratureContext::new();
6854 for &(mu, sigma) in &[(3.0, 3.0), (4.0, 4.0), (2.0, 5.0), (5.0, 5.0), (0.5, 20.0)] {
6855 let jet = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6856 .expect("wide-σ logit jet");
6857 let (rm, rd1, rd2, rd3) = logit_reference_jet_highres_simpson(mu, sigma);
6858 assert_relative_eq!(jet.mean, rm, epsilon = 1e-8, max_relative = 1e-7);
6859 assert_relative_eq!(jet.d1, rd1, epsilon = 1e-8, max_relative = 1e-6);
6860 assert_relative_eq!(jet.d2, rd2, epsilon = 1e-8, max_relative = 1e-6);
6861 assert_relative_eq!(jet.d3, rd3, epsilon = 1e-8, max_relative = 1e-6);
6862 let scalar =
6864 integrated_inverse_link_mean_and_derivative(&ctx, LinkFunction::Logit, mu, sigma)
6865 .expect("scalar logit moments");
6866 assert_relative_eq!(jet.d1, scalar.dmean_dmu, epsilon = 1e-12);
6867 assert_relative_eq!(jet.mean, scalar.mean, epsilon = 1e-12);
6868 }
6869 }
6870
6871 #[test]
6872 fn test_logit_jet_continuous_across_ghq_simpson_seam() {
6873 let ctx = QuadratureContext::new();
6881 let sigma = LOGIT_JET_GHQ_SIGMA_MAX;
6882 for mu in [-2.0, -0.5, 0.0, 0.7, 1.3, 3.0] {
6883 let ghq = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, mu, sigma)
6885 .expect("jet at seam (GHQ dispatch)");
6886 let simpson = logit_wide_sigma_jet(mu, sigma).expect("jet at seam (Simpson)");
6888 assert_relative_eq!(ghq.mean, simpson.mean, epsilon = 1e-9, max_relative = 1e-8);
6891 assert_relative_eq!(ghq.d1, simpson.d1, epsilon = 1e-9, max_relative = 1e-7);
6892 assert_relative_eq!(ghq.d2, simpson.d2, epsilon = 1e-9, max_relative = 1e-7);
6893 assert_relative_eq!(ghq.d3, simpson.d3, epsilon = 1e-8, max_relative = 1e-6);
6894 }
6895 }
6896
6897 #[test]
6898 fn test_logit_batch_uses_same_dispatchvalues() {
6899 let ctx = QuadratureContext::new();
6900 let eta = ndarray::array![-2.0, 0.0, 1.25, 35.0];
6901 let se = ndarray::array![0.1, 0.5, 1.0, 1.0];
6902 let batch_mean = logit_posterior_mean_batch(&ctx, &eta, &se)
6903 .expect("logit posterior mean batch should evaluate");
6904 let (batchmu, batch_dmu) = logit_posterior_meanwith_deriv_batch(&ctx, &eta, &se)
6905 .expect("logit posterior mean derivative batch should evaluate");
6906 for i in 0..eta.len() {
6907 let direct = integrated_inverse_link_mean_and_derivative(
6908 &ctx,
6909 LinkFunction::Logit,
6910 eta[i],
6911 se[i],
6912 )
6913 .expect("logit integrated inverse-link moments should evaluate");
6914 assert_relative_eq!(batch_mean[i], direct.mean, epsilon = 1e-12);
6915 assert_relative_eq!(batchmu[i], direct.mean, epsilon = 1e-12);
6916 assert_relative_eq!(batch_dmu[i], direct.dmean_dmu, epsilon = 1e-12);
6917 }
6918 }
6919
6920 #[test]
6921 fn exact_logit_small_se_branch_loses_tail_derivative() {
6922 let eta = 50.0_f64;
6923 let stable_z = (-eta).exp();
6924 let stable_dmu = stable_z / (1.0_f64 + stable_z).powi(2);
6925 assert!(stable_dmu > 0.0);
6926 let out = logit_posterior_meanwith_deriv_exact(eta, 0.0).expect("exact branch");
6927 let dmu = out.dmean_dmu;
6928 assert!(
6929 (dmu - stable_dmu).abs() < 1e-30,
6930 "exact logit small-se branch should use the stable derivative z/(1+z)^2 at eta={eta}; got {} vs {}",
6931 dmu,
6932 stable_dmu
6933 );
6934 }
6935
6936 #[test]
6937 fn integrated_family_moments_rejects_latent_cloglog_without_concrete_handler() {
6938 let ctx = QuadratureContext::new();
6944 let latent =
6945 gam_problem::types::LatentCLogLogState::new(0.4).expect("valid latent cloglog state");
6946 let spec =
6947 LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::LatentCLogLog(latent));
6948 let likelihood = GlmLikelihoodSpec::canonical(spec);
6949 let err = integrated_family_moments_jet(&ctx, &likelihood, 0.2, 0.5)
6950 .expect_err("latent cloglog moments should error in this dispatcher");
6951 assert!(format!("{err}").contains("LatentCLogLog"));
6952 }
6953
6954 #[test]
6955 fn integrated_family_moments_supports_stateful_sas() {
6956 let ctx = QuadratureContext::new();
6957 let sas = crate::mixture_link::state_from_sasspec(gam_problem::types::SasLinkSpec {
6958 initial_epsilon: 0.3,
6959 initial_log_delta: -0.2,
6960 })
6961 .expect("sas state should reconstruct from raw parameters");
6962 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Sas(sas));
6963 let likelihood = GlmLikelihoodSpec::canonical(spec);
6964 let out = integrated_family_moments_jet(&ctx, &likelihood, 0.2, 0.5)
6965 .expect("stateful SAS integrated moments should evaluate");
6966 assert!(out.mean.is_finite());
6967 assert!(out.d1.is_finite());
6968 assert!(out.d2.is_finite());
6969 assert!(out.d3.is_finite());
6970 assert!(out.mean > 0.0 && out.mean < 1.0);
6971 }
6972
6973 #[test]
6974 fn integrated_family_moments_supports_pure_probit_mixture() {
6975 let ctx = QuadratureContext::new();
6976 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6977 components: vec![gam_problem::types::LinkComponent::Probit],
6978 initial_rho: ndarray::Array1::<f64>::zeros(0),
6979 })
6980 .expect("single-component probit mixture state");
6981 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
6982 let likelihood = GlmLikelihoodSpec::canonical(spec);
6983 let out = integrated_family_moments_jet(&ctx, &likelihood, 0.7, 1.3)
6984 .expect("pure probit mixture integrated moments should evaluate");
6985 let exact = integrated_probit_jet(0.7, 1.3);
6986 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
6987 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
6988 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
6989 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
6990 assert_eq!(out.mode, IntegratedExpectationMode::ExactClosedForm);
6991 }
6992
6993 #[test]
6994 fn integrated_family_moments_supports_pure_logit_mixture() {
6995 let ctx = QuadratureContext::new();
6996 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
6997 components: vec![gam_problem::types::LinkComponent::Logit],
6998 initial_rho: ndarray::Array1::<f64>::zeros(0),
6999 })
7000 .expect("single-component logit mixture state");
7001 let spec = LikelihoodSpec::new(ResponseFamily::Binomial, InverseLink::Mixture(state));
7002 let likelihood = GlmLikelihoodSpec::canonical(spec);
7003 let out = integrated_family_moments_jet(&ctx, &likelihood, 1.1, 0.8)
7004 .expect("pure logit mixture integrated moments should evaluate");
7005 let exact = integrated_inverse_link_jet(&ctx, LinkFunction::Logit, 1.1, 0.8)
7006 .expect("canonical integrated logit jet");
7007 assert_relative_eq!(out.mean, exact.mean, epsilon = 1e-12);
7008 assert_relative_eq!(out.d1, exact.d1, epsilon = 1e-12);
7009 assert_relative_eq!(out.d2, exact.d2, epsilon = 1e-12);
7010 assert_relative_eq!(out.d3, exact.d3, epsilon = 1e-12);
7011 assert_eq!(out.mode, exact.mode);
7012 }
7013
7014 #[test]
7015 fn integrated_family_moments_supports_stateful_mixture() {
7016 let ctx = QuadratureContext::new();
7017 let state = crate::mixture_link::state_fromspec(&gam_problem::types::MixtureLinkSpec {
7018 components: vec![
7019 gam_problem::types::LinkComponent::Logit,
7020 gam_problem::types::LinkComponent::Probit,
7021 ],
7022 initial_rho: ndarray::array![0.35],
7023 })
7024 .expect("mixture state should reconstruct from rho");
7025 let spec = LikelihoodSpec::new(
7026 ResponseFamily::Binomial,
7027 InverseLink::Mixture(state.clone()),
7028 );
7029 let likelihood = GlmLikelihoodSpec::canonical(spec);
7030 let out = integrated_family_moments_jet(&ctx, &likelihood, 0.2, 0.5)
7031 .expect("stateful mixture integrated moments should evaluate");
7032 let direct = integrated_mixture_jet(&ctx, 0.2, 0.5, &state)
7033 .expect("direct integrated mixture jet should evaluate");
7034 assert_relative_eq!(out.mean, direct.mean, epsilon = 1e-12);
7035 assert_relative_eq!(out.d1, direct.d1, epsilon = 1e-12);
7036 assert_relative_eq!(out.d2, direct.d2, epsilon = 1e-12);
7037 assert_relative_eq!(out.d3, direct.d3, epsilon = 1e-12);
7038 assert_eq!(out.mode, direct.mode);
7039 }
7040
7041 #[test]
7042 fn integrated_family_moments_use_scale_dispersion_for_tweedie_and_gamma() {
7043 let ctx = QuadratureContext::new();
7047 let e = 0.3_f64;
7049 let se = 0.5_f64;
7050 let m = (e + 0.5 * se * se).exp();
7051
7052 let p = 1.5_f64;
7054 let phi = 2.0_f64;
7055 let tweedie = LikelihoodSpec::tweedie_log(p);
7056 let tweedie_likelihood = GlmLikelihoodSpec {
7057 spec: tweedie.clone(),
7058 scale: LikelihoodScaleMetadata::EstimatedTweediePhi { phi },
7059 };
7060 let out = integrated_family_moments_jet(&ctx, &tweedie_likelihood, e, se)
7061 .expect("tweedie integrated moments should evaluate");
7062 let expected = phi * m.powf(p);
7063 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
7064 assert_relative_eq!(out.variance / m.powf(p), phi, epsilon = 1e-12);
7066
7067 let shape = 4.0_f64;
7069 let gamma = LikelihoodSpec::gamma_log();
7070 let gamma_likelihood = GlmLikelihoodSpec {
7071 spec: gamma.clone(),
7072 scale: LikelihoodScaleMetadata::EstimatedGammaShape { shape },
7073 };
7074 let out = integrated_family_moments_jet(&ctx, &gamma_likelihood, e, se)
7075 .expect("gamma integrated moments should evaluate");
7076 let expected = m * m / shape;
7077 assert_relative_eq!(out.variance, expected, epsilon = 1e-12);
7078 assert_relative_eq!(out.variance / (m * m), 1.0 / shape, epsilon = 1e-12);
7080
7081 let poisson = LikelihoodSpec::poisson_log();
7083 let poisson_likelihood = GlmLikelihoodSpec::canonical(poisson);
7084 let out = integrated_family_moments_jet(&ctx, &poisson_likelihood, e, se)
7085 .expect("poisson integrated moments should evaluate");
7086 assert_relative_eq!(out.variance, m, epsilon = 1e-12);
7087
7088 let theta = 3.0_f64;
7090 let nb = LikelihoodSpec::negative_binomial_log(theta);
7091 let nb_likelihood = GlmLikelihoodSpec::canonical(nb);
7092 let out = integrated_family_moments_jet(&ctx, &nb_likelihood, e, se)
7093 .expect("negative-binomial integrated moments should evaluate");
7094 assert_relative_eq!(out.variance, m + m * m / theta, epsilon = 1e-12);
7095
7096 let missing_gamma = GlmLikelihoodSpec {
7098 spec: gamma,
7099 scale: LikelihoodScaleMetadata::Unspecified,
7100 };
7101 let err = integrated_family_moments_jet(&ctx, &missing_gamma, e, se)
7102 .expect_err("gamma without a shape in the scale metadata must error");
7103 assert!(
7104 format!("{err}").contains("GammaShape"),
7105 "unexpected error message: {err}"
7106 );
7107
7108 let missing_tweedie = GlmLikelihoodSpec {
7110 spec: tweedie,
7111 scale: LikelihoodScaleMetadata::Unspecified,
7112 };
7113 let err = integrated_family_moments_jet(&ctx, &missing_tweedie, e, se)
7114 .expect_err("tweedie without a φ in the scale metadata must error");
7115 assert!(
7116 format!("{err}").contains("EstimatedTweediePhi"),
7117 "unexpected error message: {err}"
7118 );
7119 }
7120
7121 #[test]
7124 fn cloglog_g_derivatives_at_zero() {
7125 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(0.0);
7126 let expected_g = 1.0 - (-1.0_f64).exp();
7128 assert_relative_eq!(g, expected_g, epsilon = 1e-14);
7129 let e_neg1 = (-1.0_f64).exp();
7131 assert_relative_eq!(g1, e_neg1, epsilon = 1e-14);
7132 assert_relative_eq!(g2, 0.0, epsilon = 1e-14);
7134 assert_relative_eq!(g3, -e_neg1, epsilon = 1e-14);
7136 assert_relative_eq!(g4, -e_neg1, epsilon = 1e-14);
7138 }
7139
7140 #[test]
7141 fn cloglog_g_derivatives_saturation() {
7142 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(50.0);
7144 assert_relative_eq!(g, 1.0, epsilon = 1e-10);
7145 assert_eq!(g1, 0.0);
7146 assert_eq!(g2, 0.0);
7147 assert_eq!(g3, 0.0);
7148 assert_eq!(g4, 0.0);
7149
7150 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(-50.0);
7152 let expected = (-50.0_f64).exp();
7153 assert_relative_eq!(g, expected, max_relative = 1e-10);
7154 assert_relative_eq!(g1, expected, max_relative = 1e-10);
7155 assert_relative_eq!(g2, expected, max_relative = 1e-10);
7157 assert_relative_eq!(g3, expected, max_relative = 1e-10);
7158 assert_relative_eq!(g4, expected, max_relative = 1e-10);
7159 }
7160
7161 #[test]
7162 fn cloglog_ghq_value_sigma_zero_matches_pointwise() {
7163 let ctx = QuadratureContext::new();
7164 for &mu in &[-2.0, -1.0, 0.0, 0.5, 1.5] {
7166 let val = cloglog_ghq_value(&ctx, mu, 0.0, 21);
7167 let (g, _, _, _, _) = cloglog_g_derivatives(mu);
7168 assert_relative_eq!(val, g, epsilon = 1e-14);
7169 }
7170 }
7171
7172 #[test]
7173 fn cloglog_ghq_value_bounded_zero_one() {
7174 let ctx = QuadratureContext::new();
7175 for &mu in &[-5.0, -2.0, 0.0, 1.0, 3.0, 10.0] {
7177 for &sigma in &[0.1, 0.5, 1.0, 2.0, 5.0] {
7178 let val = cloglog_ghq_value(&ctx, mu, sigma, 31);
7179 assert!((0.0..=1.0).contains(&val), "L({mu},{sigma}) = {val}");
7180 }
7181 }
7182 }
7183
7184 #[test]
7185 fn cloglog_ghq_derivatives_sigma_zero_matches_pointwise() {
7186 let ctx = QuadratureContext::new();
7187 let mu = 0.3;
7188 let d = cloglog_ghq_derivatives(&ctx, mu, 0.0, 21);
7189 let (g, g1, g2, g3, g4) = cloglog_g_derivatives(mu);
7190 assert_relative_eq!(d.l, g, epsilon = 1e-14);
7191 assert_relative_eq!(d.l_mu, g1, epsilon = 1e-14);
7192 assert_relative_eq!(d.l_mumu, g2, epsilon = 1e-14);
7193 assert_relative_eq!(d.l_mumumu, g3, epsilon = 1e-14);
7194 assert_relative_eq!(d.l_mumumumu, g4, epsilon = 1e-14);
7195
7196 assert_eq!(d.l_sigma, 0.0);
7198 assert_eq!(d.l_musigma, 0.0);
7199 assert_eq!(d.l_mumusigma, 0.0);
7200 assert_eq!(d.l_mumumusigma, 0.0);
7201 assert_eq!(d.l_sigmasigmasigma, 0.0);
7202 assert_eq!(d.l_musigmasigmasigma, 0.0);
7203
7204 assert_relative_eq!(d.l_sigmasigma, g2, epsilon = 1e-14);
7207 assert_relative_eq!(d.l_musigmasigma, g3, epsilon = 1e-14);
7208 assert_relative_eq!(d.l_mumusigmasigma, g4, epsilon = 1e-14);
7209 assert_relative_eq!(d.l_sigmasigmasigmasigma, 3.0 * g4, epsilon = 1e-14);
7210 }
7211
7212 #[test]
7213 fn cloglog_ghq_derivatives_finite_difference_mu() {
7214 let ctx = QuadratureContext::new();
7216 let mu = 0.5;
7217 let sigma = 0.8;
7218 let h = 1e-6;
7219 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
7220 let l_plus = cloglog_ghq_value(&ctx, mu + h, sigma, 31);
7221 let l_minus = cloglog_ghq_value(&ctx, mu - h, sigma, 31);
7222 let fd_mu = (l_plus - l_minus) / (2.0 * h);
7223 assert_relative_eq!(d.l_mu, fd_mu, epsilon = 1e-5);
7224
7225 let d_plus = cloglog_ghq_derivatives(&ctx, mu + h, sigma, 31);
7227 let d_minus = cloglog_ghq_derivatives(&ctx, mu - h, sigma, 31);
7228 let fd_mumu = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
7229 assert_relative_eq!(d.l_mumu, fd_mumu, epsilon = 1e-4);
7230 }
7231
7232 #[test]
7233 fn cloglog_ghq_derivatives_finite_difference_sigma() {
7234 let ctx = QuadratureContext::new();
7236 let mu = 0.2;
7237 let sigma = 1.0;
7238 let h = 1e-6;
7239 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
7240 let l_plus = cloglog_ghq_value(&ctx, mu, sigma + h, 31);
7241 let l_minus = cloglog_ghq_value(&ctx, mu, sigma - h, 31);
7242 let fd_sigma = (l_plus - l_minus) / (2.0 * h);
7243 assert_relative_eq!(d.l_sigma, fd_sigma, epsilon = 1e-5);
7244 }
7245
7246 #[test]
7247 fn cloglog_ghq_derivatives_finite_difference_cross() {
7248 let ctx = QuadratureContext::new();
7250 let mu = -0.5;
7251 let sigma = 0.6;
7252 let h = 1e-6;
7253 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 31);
7254 let d_plus = cloglog_ghq_derivatives(&ctx, mu, sigma + h, 31);
7255 let d_minus = cloglog_ghq_derivatives(&ctx, mu, sigma - h, 31);
7256 let fd_musigma = (d_plus.l_mu - d_minus.l_mu) / (2.0 * h);
7257 assert_relative_eq!(d.l_musigma, fd_musigma, epsilon = 1e-4);
7258 }
7259
7260 #[test]
7261 fn cloglog_ghq_l_mu_nonnegative() {
7262 let ctx = QuadratureContext::new();
7264 for &mu in &[-3.0, -1.0, 0.0, 1.0, 3.0] {
7265 for &sigma in &[0.1, 0.5, 1.0, 2.0] {
7266 let d = cloglog_ghq_derivatives(&ctx, mu, sigma, 21);
7267 assert!(
7268 d.l_mu >= -1e-14,
7269 "L_mu should be non-negative at mu={mu}, sigma={sigma}: got {}",
7270 d.l_mu
7271 );
7272 }
7273 }
7274 }
7275
7276 #[test]
7277 fn cloglog_ghq_adaptive_matches_explicit() {
7278 let ctx = QuadratureContext::new();
7279 let mu = 0.7;
7280 let sigma = 1.2;
7281 let adaptive = cloglog_ghq_derivatives_adaptive(&ctx, mu, sigma);
7282 let n = adaptive_point_count_from_sd(sigma);
7283 let explicit = cloglog_ghq_derivatives(&ctx, mu, sigma, n);
7284 assert_relative_eq!(adaptive.l, explicit.l, epsilon = 1e-15);
7285 assert_relative_eq!(adaptive.l_mu, explicit.l_mu, epsilon = 1e-15);
7286 assert_relative_eq!(adaptive.l_sigma, explicit.l_sigma, epsilon = 1e-15);
7287 assert_relative_eq!(adaptive.l_mumu, explicit.l_mumu, epsilon = 1e-15);
7288 }
7289
7290 #[test]
7291 fn cloglog_ghq_value_matches_mathematical_target_in_central_regime() {
7292 let ctx = QuadratureContext::new();
7293 for &mu in &[-1.0, 0.0, 0.5, 2.0] {
7294 for &sigma in &[0.1, 0.5, 1.0] {
7295 let ghq = cloglog_ghq_value(&ctx, mu, sigma, 51);
7296 let (expected_mean, _) = cloglog_reference_mean_and_derivative(mu, sigma);
7297 assert_relative_eq!(ghq, expected_mean, epsilon = 1e-12, max_relative = 2e-8);
7298 }
7299 }
7300 }
7301
7302 #[test]
7305 fn cloglog_negative_tail_mean_matches_exact_near_transition() {
7306 let eta: f64 = -30.0;
7310 let exact = {
7311 let ex = eta.exp();
7312 -(-ex).exp_m1()
7313 };
7314 let tail = cloglog_negative_tail_mean(eta);
7315 assert!(
7316 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7317 "tail mean at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7318 );
7319 }
7320
7321 #[inline]
7322 fn cloglog_negative_tail_derivative(eta: f64) -> f64 {
7323 if eta < -745.0 {
7325 0.0
7326 } else {
7327 let ex = safe_exp(eta);
7328 (ex * (-ex).exp()).max(0.0)
7329 }
7330 }
7331
7332 #[test]
7333 fn cloglog_negative_tail_derivative_matches_exact_near_transition() {
7334 let eta: f64 = -30.0;
7336 let ex = eta.exp();
7337 let exact = ex * (-ex).exp();
7338 let tail = cloglog_negative_tail_derivative(eta);
7339 assert!(
7340 (exact - tail).abs() < 1e-26 * exact.abs().max(1e-300),
7341 "tail derivative at η={eta}: exact={exact:.6e} tail={tail:.6e}"
7342 );
7343 }
7344
7345 #[test]
7346 fn cloglog_negative_tail_degenerate_branch_matches_target_near_transition() {
7347 let ctx = QuadratureContext::default();
7348 let sigma = 0.0;
7349 for &mu in &[-30.001, -30.0, -29.999] {
7350 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7351 assert_relative_eq!(
7352 out.mean,
7353 cloglog_mean_exact(mu),
7354 epsilon = 1e-28,
7355 max_relative = 1e-15
7356 );
7357 assert_relative_eq!(
7358 out.dmean_dmu,
7359 cloglog_mean_d1_exact(mu),
7360 epsilon = 1e-28,
7361 max_relative = 1e-15
7362 );
7363 }
7364 }
7365
7366 #[test]
7367 fn cloglog_negative_tail_small_sigma_branch_matches_target_near_transition() {
7368 let ctx = QuadratureContext::default();
7369 let sigma = 0.1;
7370 for &mu in &[-30.001, -30.0, -29.999] {
7371 let out = cloglog_posterior_meanwith_deriv_controlled(&ctx, mu, sigma);
7372 let (expected_mean, expected_deriv) = cloglog_reference_mean_and_derivative(mu, sigma);
7373 assert_relative_eq!(
7374 out.mean,
7375 expected_mean,
7376 epsilon = 1e-24,
7377 max_relative = 1e-10
7378 );
7379 assert_relative_eq!(
7380 out.dmean_dmu,
7381 expected_deriv,
7382 epsilon = 1e-24,
7383 max_relative = 1e-10
7384 );
7385 }
7386 }
7387
7388 fn ref_cholesky_heap(cov: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
7392 let n = cov.len();
7393 if n == 0 || cov.iter().any(|r| r.len() != n) {
7394 return None;
7395 }
7396 let mut base = cov.to_vec();
7397 for retry in 0..8 {
7398 let jitter = if retry == 0 {
7399 0.0
7400 } else {
7401 1e-12 * 10f64.powi(retry - 1)
7402 };
7403 if jitter > 0.0 {
7404 for i in 0..n {
7405 base[i][i] = cov[i][i] + jitter;
7406 }
7407 }
7408 let mut l = vec![vec![0.0_f64; n]; n];
7409 let mut ok = true;
7410 for i in 0..n {
7411 for j in 0..=i {
7412 let mut sum = base[i][j];
7413 for k in 0..j {
7414 sum -= l[i][k] * l[j][k];
7415 }
7416 if i == j {
7417 if !sum.is_finite() || sum <= 0.0 {
7418 ok = false;
7419 break;
7420 }
7421 l[i][j] = sum.sqrt();
7422 } else {
7423 l[i][j] = sum / l[j][j];
7424 }
7425 }
7426 if !ok {
7427 break;
7428 }
7429 }
7430 if ok {
7431 return Some(l);
7432 }
7433 }
7434 None
7435 }
7436
7437 #[test]
7438 fn cholesky_static_matches_heap_d2() {
7439 let cases: &[[[f64; 2]; 2]] = &[
7442 [[1.0, 0.0], [0.0, 1.0]],
7443 [[2.5, 0.3], [0.3, 0.75]],
7444 [[1.0, 0.9999], [0.9999, 1.0]],
7445 [[1e-10, 0.0], [0.0, 1e-10]],
7446 [[4.0, -1.5], [-1.5, 2.25]],
7447 ];
7448 for cov in cases {
7449 let stack = cholesky_static_with_jitter::<2>(cov).expect("stack cholesky");
7450 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7451 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7452 for i in 0..2 {
7453 for j in 0..2 {
7454 assert_eq!(
7455 stack[i][j].to_bits(),
7456 heap[i][j].to_bits(),
7457 "mismatch at ({i},{j}) for cov={cov:?}"
7458 );
7459 }
7460 }
7461 }
7462 }
7463
7464 #[test]
7465 fn cholesky_static_matches_heap_d3() {
7466 let cases: &[[[f64; 3]; 3]] = &[
7467 [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
7468 [[2.0, 0.5, 0.1], [0.5, 1.5, -0.2], [0.1, -0.2, 0.8]],
7469 [[4.0, 1.0, 0.5], [1.0, 3.0, 0.25], [0.5, 0.25, 2.0]],
7470 ];
7471 for cov in cases {
7472 let stack = cholesky_static_with_jitter::<3>(cov).expect("stack cholesky");
7473 let heap_in: Vec<Vec<f64>> = cov.iter().map(|r| r.to_vec()).collect();
7474 let heap = ref_cholesky_heap(&heap_in).expect("heap cholesky");
7475 for i in 0..3 {
7476 for j in 0..3 {
7477 assert_eq!(
7478 stack[i][j].to_bits(),
7479 heap[i][j].to_bits(),
7480 "mismatch at ({i},{j}) for cov={cov:?}"
7481 );
7482 }
7483 }
7484 }
7485 }
7486
7487 #[test]
7488 fn cholesky_static_d1() {
7489 let l = cholesky_static_with_jitter::<1>(&[[2.25]]).expect("d=1");
7490 assert_eq!(l[0][0], 1.5);
7491 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e-13]]).is_some());
7501 assert!(cholesky_static_with_jitter::<1>(&[[-1.0e3]]).is_none());
7504 }
7505}
7506
7507#[cfg(test)]
7508mod log_survival_panel_2714_tests {
7509 use super::{
7510 LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER, LOG_SURVIVAL_PANEL_MAX_NODES,
7511 LOG_SURVIVAL_PANEL_MIN_NODES, LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION, LogSurvivalBranch,
7512 QuadratureContext, log_survival_jet, log_survival_panel,
7513 };
7514
7515 const LOG_SURVIVAL_REFERENCE: &[(f64, f64, f64)] = &[
7524 (-30.0, 0.05, -9.369_327_311_241_221e-14),
7525 (-20.0, 0.002, -2.061_157_744_749_916_5e-9),
7526 (-20.0, 2.0, -1.522_997_352_870_166_6e-8),
7527 (-8.0, 0.15, -3.392_565_812_965_201_6e-4),
7528 (-8.0, 8.0, -1.981_213_453_495_813_4e-1),
7529 (-3.0, 0.02, -4.979_653_073_925_083_3e-2),
7530 (-3.0, 1.0, -7.722_156_176_750_611e-2),
7531 (-1.0, 0.005, -3.678_823_479_542_769_3e-1),
7532 (0.0, 0.05, -9.999_992_216_748_937e-1),
7533 (0.0, 20.0, -7.163_521_358_485_562e-1),
7534 (1.0, 0.15, -2.668_055_285_814_174_4e0),
7535 (1.8, 0.15, -5.742_720_404_435_083e0),
7536 (1.8, 0.5, -4.257_254_029_226_836e0),
7537 (3.2, 0.005, -2.452_531_812_111_085_3e1),
7538 (3.2, 0.15, -2.014_631_867_969_775_7e1),
7539 (3.2, 4.0, -1.691_825_149_278_477e0),
7540 (5.0, 0.02, -1.442_777_351_752_632_2e2),
7541 (5.0, 1.0, -1.128_111_355_370_178_3e1),
7542 (8.0, 0.002, -2.963_400_229_027_156e3),
7543 (8.0, 0.05, -1.113_594_437_164_973_2e3),
7544 (8.0, 0.5, -7.097_988_851_759_84e1),
7545 (8.0, 60.0, -8.137_859_923_376_214e-1),
7546 (12.0, 0.002, -1.289_824_339_085_843_7e5),
7547 (12.0, 0.15, -1.181_338_131_091_633_2e3),
7548 (12.0, 1.0, -5.819_669_561_598_102e1),
7549 (20.0, 0.05, -3.135_653_820_340_577_3e4),
7550 (-50.0, 8.0, -7.328_008_793_232_462e-10),
7551 (-18.0, 0.001, -1.522_998_735_970_428_8e-8),
7552 ];
7553
7554 #[test]
7575 fn log_survival_matches_a_high_precision_reference_2714() {
7576 let ctx = QuadratureContext::new();
7577 let mut worst = 0.0_f64;
7578 let mut worst_at = (f64::NAN, f64::NAN);
7579 for &(mu, sigma, expected) in LOG_SURVIVAL_REFERENCE {
7580 let got = log_survival_jet(&ctx, mu, sigma, 0).log_survival;
7581 let error = (got - expected).abs() / expected.abs().max(1.0);
7582 if error > worst {
7583 worst = error;
7584 worst_at = (mu, sigma);
7585 }
7586 assert!(
7587 error <= 1.0e-13,
7588 "#2714: ln S({mu}, {sigma}) = {got:.17e} against reference \
7589 {expected:.17e} (relative {error:.3e}); one log-space panel has \
7590 to hold the f64 floor at every (mu, sigma), because a routed \
7591 surface whose error is a step function of (mu, sigma) is the \
7592 defect this replaced"
7593 );
7594 }
7595 assert!(
7598 worst > 0.0,
7599 "#2714: every reference row reproduced bit-exactly, which means the \
7600 table is not testing the quadrature"
7601 );
7602 let (mu, sigma) = worst_at;
7603 println!("[2714] worst relative ln S error {worst:.3e} at (mu={mu}, sigma={sigma})");
7604 }
7605
7606 #[test]
7616 fn log_survival_tower_is_the_derivative_of_the_value_2714() {
7617 let ctx = QuadratureContext::new();
7618 let mut worst = 0.0_f64;
7619 for &(mu, sigma) in &[
7620 (-8.0, 0.15),
7621 (-3.0, 0.5),
7622 (0.0, 1.0),
7623 (1.8, 0.15),
7624 (3.2, 0.15),
7625 (3.2, 2.0),
7626 (5.0, 1.0),
7627 (8.0, 4.0),
7628 (0.0, 8.0),
7629 (-3.0, 20.0),
7630 ] {
7631 let jet = log_survival_jet(&ctx, mu, sigma, 1);
7632 let first = jet.scaled_mu_derivatives[1];
7633 let analytic = first.sign * (first.log_abs - jet.log_survival).exp() / sigma;
7635 let step = (f64::EPSILON.cbrt()) * (1.0 + mu.abs());
7636 let up = log_survival_jet(&ctx, mu + step, sigma, 0).log_survival;
7637 let down = log_survival_jet(&ctx, mu - step, sigma, 0).log_survival;
7638 let numeric = (up - down) / (2.0 * step);
7639 let error = (analytic - numeric).abs() / analytic.abs().max(1.0e-3);
7640 worst = worst.max(error);
7641 assert!(
7642 error <= 1.0e-7,
7643 "#2714: d/dmu ln S at (mu={mu}, sigma={sigma}) is {analytic:.12e} \
7644 analytically and {numeric:.12e} by central difference of the \
7645 shipped value (relative {error:.3e}). A tower that is not the \
7646 derivative of the value it ships with is what collapses the \
7647 joint-Newton trust region."
7648 );
7649 }
7650 println!("[2714] worst |analytic - FD| / |analytic| = {worst:.3e}");
7651 }
7652
7653 #[test]
7670 fn log_survival_tower_is_independent_of_the_requested_order_2714() {
7671 let ctx = QuadratureContext::new();
7672 let mut compared = 0usize;
7673 for &(mu, sigma) in &[
7674 (-20.0, 0.002),
7675 (-8.0, 0.15),
7676 (-3.0, 0.5),
7677 (0.0, 1.0),
7678 (1.8, 0.15),
7679 (3.2, 0.15),
7680 (3.2, 2.0),
7681 (5.0, 1.0),
7682 (8.0, 4.0),
7683 (0.0, 8.0),
7684 (12.0, 0.02),
7685 (-3.0, 20.0),
7686 ] {
7687 let reference = log_survival_jet(&ctx, mu, sigma, LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER);
7688 for requested in 1..=LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER {
7689 let jet = log_survival_jet(&ctx, mu, sigma, requested);
7690 for order in 1..=requested {
7691 let got = jet.scaled_mu_derivatives[order];
7692 let want = reference.scaled_mu_derivatives[order];
7693 assert_eq!(
7694 (got.log_abs.to_bits(), got.sign.to_bits()),
7695 (want.log_abs.to_bits(), want.sign.to_bits()),
7696 "#2714: sigma^{order} d^{order} S/d mu^{order} at (mu={mu}, \
7697 sigma={sigma}) is {:.17e} (sign {}) when {requested} orders \
7698 are requested and {:.17e} (sign {}) when \
7699 {LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER} are. The tower must \
7700 not depend on how much of it the caller wants.",
7701 got.log_abs, got.sign, want.log_abs, want.sign
7702 );
7703 compared += 1;
7704 }
7705 }
7706 }
7707 assert!(
7710 compared >= 12 * (LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER
7711 * (LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER + 1)
7712 / 2),
7713 "#2714: only {compared} tower entries were compared"
7714 );
7715 println!("[2714] tower order-independence: {compared} entries bit-identical");
7716 }
7717
7718 #[test]
7726 fn log_survival_tower_certification_is_a_prefix_2714() {
7727 let ctx = QuadratureContext::new();
7728 let mut saw_truncation = false;
7729 for &(mu, sigma, ..) in LOG_SURVIVAL_TOWER_REFERENCE {
7730 let jet = log_survival_jet(&ctx, mu, sigma, LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER);
7731 let Some(prefix) = jet.certified_prefix_order(LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER)
7732 else {
7733 continue;
7734 };
7735 for requested in 0..=LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER {
7737 assert_eq!(
7738 jet.certified_prefix_order(requested),
7739 Some(prefix.min(requested)),
7740 "#2714: the certified prefix at (mu={mu}, sigma={sigma}) is \
7741 {prefix} but reading it against a request of {requested} \
7742 disagrees — the prefix must be a property of the tower, \
7743 truncated by the request and not decided by it"
7744 );
7745 assert_eq!(
7749 jet.certified_scaled_mu_derivatives(requested).is_some(),
7750 requested <= prefix,
7751 "#2714: the whole-tower and prefix readings disagree at \
7752 (mu={mu}, sigma={sigma}, requested={requested})"
7753 );
7754 }
7755 if prefix < LOG_SURVIVAL_MAX_MU_DERIVATIVE_ORDER {
7756 saw_truncation = true;
7757 }
7758 }
7759 assert!(
7763 saw_truncation,
7764 "#2714: no reference row produced a truncated certified prefix, so \
7765 this test did not exercise the behaviour it exists for"
7766 );
7767 }
7768
7769 const LOG_SURVIVAL_TOWER_REFERENCE: &[(f64, f64, [f64; 5], [f64; 5])] = &[
7783 (-20.0, 0.002, [-2.0611577447499165e-9, -26.214606100483358, -32.429214200966715, -38.643822303511239, -44.858430410178095], [1.0, -1.0, -1.0, -1.0, -1.0]),
7784 (-8.0, 0.15, [-0.00033925658129652018, -9.8862169612318004, -11.783683981303588, -13.681498274811226, -15.580007954561913], [1.0, -1.0, -1.0, -1.0, -1.0]),
7785 (-8.0, 8.0, [-0.19812134534958133, -1.3528211702547069, -1.452147542328205, -3.240571743395306, -0.71024921821004234], [1.0, -1.0, -1.0, 1.0, 1.0]),
7786 (-3.0, 0.02, [-0.049796530739250832, -6.9616394585654159, -10.92476204736755, -14.944640460549261, -19.104083600790619], [1.0, -1.0, -1.0, -1.0, -1.0]),
7787 (-3.0, 0.5, [-0.055971877835218903, -3.6398579829672956, -4.4066491238893504, -5.2575250298329575, -6.3338354346007994], [1.0, -1.0, -1.0, -1.0, -1.0]),
7788 (-3.0, 20.0, [-0.60124121574363753, -0.9283108111246211, -3.0420009372801215, -0.94719003039898553, -1.9515751519096107], [1.0, -1.0, -1.0, 1.0, 1.0]),
7789 (-1.0, 0.005, [-0.36788234795427693, -6.666196411633682, -12.423205395414024, -20.715002885443461, -22.768239058212055], [1.0, -1.0, -1.0, -1.0, 1.0]),
7790 (0.0, 0.05, [-0.99999922167489372, -3.9969814880950871, -13.681704026824476, -9.9897069128682089, -12.994203213800263], [1.0, -1.0, 1.0, 1.0, 1.0]),
7791 (0.0, 1.0, [-0.96297240050030377, -1.3514828821346528, -3.4776689251068834, -2.0721333282379262, -4.077937287249247], [1.0, -1.0, 1.0, 1.0, -1.0]),
7792 (0.0, 8.0, [-0.75101079589424262, -0.93383523622026295, -3.6184199118831959, -0.96301212066547202, -2.5659122387330754], [1.0, -1.0, 1.0, 1.0, -1.0]),
7793 (1.0, 0.15, [-2.6680552858141743, -3.6128621614846425, -5.0108306306587141, -8.9368697272493613, -7.5728553818569553], [1.0, -1.0, 1.0, -1.0, -1.0]),
7794 (1.8, 0.15, [-5.7427204044350835, -5.9515476490674363, -6.340239096657947, -7.0035313750765859, -8.252090685678664], [1.0, -1.0, 1.0, -1.0, 1.0]),
7795 (1.8, 0.5, [-4.2572540292268353, -3.8392878811891043, -3.6246169834513785, -3.7636817578612526, -5.1584862006168726], [1.0, -1.0, 1.0, -1.0, 1.0]),
7796 (3.2, 0.005, [-24.525318121110852, -26.624235940415906, -28.764769508877318, -30.95061498411088, -33.186272517159154], [1.0, -1.0, 1.0, -1.0, 1.0]),
7797 (3.2, 0.15, [-20.146318679697758, -19.21570557965608, -18.328812156432285, -17.490308387278464, -16.706077981177801], [1.0, -1.0, 1.0, -1.0, 1.0]),
7798 (3.2, 2.0, [-2.9809436264637528, -2.3631698967287561, -1.9865287099827521, -2.1204144743340792, -2.9045009851726548], [1.0, -1.0, 1.0, -1.0, -1.0]),
7799 (3.2, 4.0, [-1.691825149278477, -1.3630396002384465, -1.5160548511506286, -2.9955755434621241, -0.80260762836904534], [1.0, -1.0, 1.0, 1.0, -1.0]),
7800 (5.0, 1.0, [-11.281113553701783, -9.9522520507276806, -8.6794206852298702, -7.4721966945023752, -6.3443308739944037], [1.0, -1.0, 1.0, -1.0, 1.0]),
7801 (8.0, 0.05, [-1113.5944371649733, -1110.1523256471932, -1106.7108385495186, -1103.2699768904941, -1099.8297416923265], [1.0, -1.0, 1.0, -1.0, 1.0]),
7802 (8.0, 0.5, [-70.979888517598404, -68.673129437259724, -66.374650578770375, -64.084648631557826, -61.803330805347407], [1.0, -1.0, 1.0, -1.0, 1.0]),
7803 (8.0, 4.0, [-3.9201018255002301, -3.0676533191912935, -2.3815942308869494, -1.9700438387194749, -2.2785027569363032], [1.0, -1.0, 1.0, -1.0, 1.0]),
7804 (8.0, 60.0, [-0.8137859923376214, -0.92937929368933511, -2.8751091738859909, -0.95047056774657423, -1.7838093439866154], [1.0, -1.0, 1.0, 1.0, -1.0]),
7805 (12.0, 0.002, [-128982.43390858436, -128977.07394828311, -128971.7139945779, -128966.35404746883, -128960.99410695599], [1.0, -1.0, 1.0, -1.0, 1.0]),
7806 (12.0, 1.0, [-58.196695615981018, -55.917604210289579, -53.648028066976042, -51.388233991426026, -49.138505716496612], [1.0, -1.0, 1.0, -1.0, 1.0]),
7807 (-50.0, 8.0, [-7.3280087932324622e-10, -19.248520695018656, -17.486423137639009, -15.749919102771783, -14.041433160139017], [1.0, -1.0, -1.0, -1.0, -1.0]),
7808 (20.0, 0.05, [-31356.538203405772, -31351.094833958202, -31345.651481726192, -31340.208146710605, -31334.764828912308], [1.0, -1.0, 1.0, -1.0, 1.0]),
7809 ];
7810
7811 #[test]
7826 fn log_survival_tower_matches_a_high_precision_reference_2714() {
7827 let ctx = QuadratureContext::new();
7828 let mut worst_certified = 0.0_f64;
7829 let mut worst_certified_at = (f64::NAN, f64::NAN, 0usize);
7830 let mut certified_rows = 0usize;
7831 let mut refused_rows = 0usize;
7832 for &(mu, sigma, expected_log, expected_sign) in LOG_SURVIVAL_TOWER_REFERENCE {
7833 let jet = log_survival_jet(&ctx, mu, sigma, 4);
7834 let certified = jet.certified_scaled_mu_derivatives(4).is_some();
7835 if certified {
7836 certified_rows += 1;
7837 } else {
7838 refused_rows += 1;
7839 }
7840 if !certified {
7841 continue;
7850 }
7851 for order in 0..=4 {
7852 let entry = jet.scaled_mu_derivatives[order];
7853 assert_eq!(
7854 entry.sign, expected_sign[order],
7855 "#2714: sigma^{order} d^{order} S/d mu^{order} at (mu={mu}, \
7856 sigma={sigma}) is certified and has sign {} against the \
7857 high-precision {}. A sign error in the tower is not an \
7858 accuracy question — it reverses the curvature every \
7859 latent-survival row hands the joint Newton.",
7860 entry.sign, expected_sign[order]
7861 );
7862 let error = (entry.log_abs - expected_log[order]).abs();
7869 let bar = 1.0e-13 * expected_log[order].abs().max(1.0);
7870 assert!(
7871 error <= bar,
7872 "#2714: sigma^{order} d^{order} S/d mu^{order} at (mu={mu}, \
7873 sigma={sigma}) is CERTIFIED by its own measured \
7874 cancellation ({:.3} nats) and yet lands {error:.3e} from \
7875 the reference against a bar of {bar:.3e}. Certification \
7876 means 'this entry is at the working floor'; if it is not, \
7877 the gate that replaced the sigma >= 8 constant is weaker \
7878 than the constant was.",
7879 entry.log_cancellation
7880 );
7881 let relative = error / expected_log[order].abs().max(1.0);
7882 if relative > worst_certified {
7883 worst_certified = relative;
7884 worst_certified_at = (mu, sigma, order);
7885 }
7886 }
7887 }
7888 assert!(
7892 certified_rows > 0 && refused_rows > 0,
7893 "#2714: the tower reference must straddle the certification gate — \
7894 got {certified_rows} certified and {refused_rows} refused rows"
7895 );
7896 let (mu, sigma, order) = worst_certified_at;
7897 println!(
7898 "[2714] tower: {certified_rows} certified / {refused_rows} refused; \
7899 worst certified error {worst_certified:.3e} at (mu={mu}, sigma={sigma}, j={order})"
7900 );
7901 }
7902
7903 #[test]
7920 fn log_survival_tower_cancellation_bounds_its_own_error_2714() {
7921 const ERROR_MODEL_SLACK: f64 = 128.0;
7928 let ctx = QuadratureContext::new();
7935 let mut worst_ratio = 0.0_f64;
7936 let mut worst_at = (f64::NAN, f64::NAN, 0usize);
7937 let mut observed_max_cancellation = 0.0_f64;
7938 for &(mu, sigma, expected_log, _) in LOG_SURVIVAL_TOWER_REFERENCE {
7939 let jet = log_survival_jet(&ctx, mu, sigma, 4);
7940 for order in 0..=4 {
7941 let entry = jet.scaled_mu_derivatives[order];
7942 if !entry.log_cancellation.is_finite() {
7943 continue;
7944 }
7945 observed_max_cancellation =
7946 observed_max_cancellation.max(entry.log_cancellation);
7947 let predicted = ERROR_MODEL_SLACK
7948 * f64::EPSILON
7949 * entry.log_cancellation.exp().max(entry.log_abs.abs().max(1.0));
7950 let error = (entry.log_abs - expected_log[order]).abs();
7951 let ratio = error / predicted;
7952 if ratio > worst_ratio {
7953 worst_ratio = ratio;
7954 worst_at = (mu, sigma, order);
7955 }
7956 assert!(
7957 error <= predicted,
7958 "#2714: the tower's own cancellation bar is not honest at \
7959 (mu={mu}, sigma={sigma}, j={order}): it reports \
7960 {:.4} nats of cancellation, which the signed-log-sum-exp \
7961 error model prices at {predicted:.3e}, and the achieved \
7962 error against the high-precision reference is \
7963 {error:.3e}. The admission bar \
7964 LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION is derived FROM \
7965 that model, so a model that under-predicts makes the bar \
7966 meaningless.",
7967 entry.log_cancellation
7968 );
7969 }
7970 }
7971 assert!(
7972 observed_max_cancellation > LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION,
7973 "#2714: no row of the reference reaches the admission bar \
7974 ({LOG_SURVIVAL_TOWER_MAX_LOG_CANCELLATION} nats), so this test \
7975 never grades the model in the regime the bar exists for — worst \
7976 observed cancellation was {observed_max_cancellation:.3}"
7977 );
7978 let (mu, sigma, order) = worst_at;
7979 println!(
7980 "[2714] cancellation bar: worst achieved/predicted = {worst_ratio:.3e} \
7981 at (mu={mu}, sigma={sigma}, j={order}); max cancellation \
7982 {observed_max_cancellation:.3} nats"
7983 );
7984 }
7985
7986 #[test]
8002 fn log_survival_has_no_routing_step_in_it_2714() {
8003 let ctx = QuadratureContext::new();
8004 let step = 1.0e-3_f64;
8005 for &sigma in &[0.05_f64, 0.15, 0.5, 1.0, 4.0, 8.0, 20.0] {
8006 let mut worst = 0.0_f64;
8007 let mut worst_mu = f64::NAN;
8008 let value = |mu: f64| log_survival_jet(&ctx, mu, sigma, 0).log_survival;
8009 let mut index = -100_i32;
8010 while index <= 100 {
8011 let mu = f64::from(index) * 0.1;
8012 let mid = value(mu);
8013 let fine =
8014 (value(mu - step) - 2.0 * mid + value(mu + step)) / (step * step);
8015 let coarse = (value(mu - 2.0 * step) - 2.0 * mid + value(mu + 2.0 * step))
8016 / (4.0 * step * step);
8017 let disagreement = (fine - coarse).abs() / (1.0 + fine.abs());
8018 if disagreement > worst {
8019 worst = disagreement;
8020 worst_mu = mu;
8021 }
8022 index += 1;
8023 }
8024 assert!(
8025 worst <= 2.0e-6,
8026 "#2714: the h and 2h second differences of ln S(., sigma={sigma}) \
8027 disagree by {worst:.3e} (relative) near mu={worst_mu}, i.e. a jump \
8028 of about {:.3e} in the value. The surface must be smooth across \
8029 every internal routing decision, because that is precisely what \
8030 the joint-Newton accept test differentiates.",
8031 worst * step * step / 0.75
8032 );
8033 println!("[2714] sigma={sigma}: worst Richardson disagreement {worst:.3e}");
8034 }
8035 }
8036
8037 #[test]
8051 fn log_survival_panel_node_count_is_live_and_bounded_2714() {
8052 let mut counts = Vec::new();
8053 for &sigma in &[1.0e-4_f64, 1.0e-2, 0.15, 1.0, 8.0, 60.0, 1.0e3] {
8054 let panel = log_survival_panel(LogSurvivalBranch::Survival, 1.8, sigma, 0);
8055 assert!(
8056 panel.nodes >= LOG_SURVIVAL_PANEL_MIN_NODES
8057 && panel.nodes <= LOG_SURVIVAL_PANEL_MAX_NODES,
8058 "#2714: node count {} at sigma={sigma} is outside \
8059 [{LOG_SURVIVAL_PANEL_MIN_NODES}, {LOG_SURVIVAL_PANEL_MAX_NODES}]",
8060 panel.nodes
8061 );
8062 assert!(
8063 !panel.nodes.is_multiple_of(2),
8064 "#2714: the Clenshaw-Curtis panel needs an odd node count so the \
8065 grid is symmetric and gap-free (sigma={sigma} gave {})",
8066 panel.nodes
8067 );
8068 assert!(
8069 panel.z_lo < panel.z_hi,
8070 "#2714: degenerate panel at sigma={sigma}"
8071 );
8072 counts.push(panel.nodes);
8073 }
8074 let (low, high) = (
8075 *counts.iter().min().expect("non-empty"),
8076 *counts.iter().max().expect("non-empty"),
8077 );
8078 assert!(
8079 high >= 4 * low,
8080 "#2714: the node ladder must actually MOVE with sigma — it spans \
8081 {low}..{high} over sigma in [1e-4, 1e3], and a ladder that does not \
8082 move is the #2469 inertness that made the escape-hatch quadrature \
8083 unable to resolve its own transition"
8084 );
8085 }
8086
8087 #[test]
8091 fn log_survival_panel_brackets_its_own_peak_2714() {
8092 for &(mu, sigma) in &[
8093 (-20.0, 0.002),
8094 (-3.0, 0.15),
8095 (0.0, 1.0),
8096 (3.2, 0.15),
8097 (12.0, 0.002),
8098 (0.0, 60.0),
8099 ] {
8100 for branch in [LogSurvivalBranch::Survival, LogSurvivalBranch::Complement] {
8101 let panel = log_survival_panel(branch, mu, sigma, 0);
8102 let mid = 0.5 * (panel.z_lo + panel.z_hi);
8103 let at_mid = branch.log_integrand(mu, sigma, mid);
8104 let at_lo = branch.log_integrand(mu, sigma, panel.z_lo);
8105 let at_hi = branch.log_integrand(mu, sigma, panel.z_hi);
8106 assert!(
8107 at_lo < at_mid && at_hi < at_mid,
8108 "#2714: {branch:?} panel at (mu={mu}, sigma={sigma}) does not \
8109 bracket its own maximum: L(lo)={at_lo:.6e}, L(mid)={at_mid:.6e}, \
8110 L(hi)={at_hi:.6e}"
8111 );
8112 }
8113 }
8114 }
8115}