Struct LatentZConditionalCalibration 

pub struct LatentZConditionalCalibration {
    pub mean_coeffs: Vec<f64>,
    pub var_coeffs: Vec<f64>,
    pub basis_ncols: usize,
    pub var_floor: f64,
    pub homoskedastic_var: f64,
    pub post_mean: f64,
    pub post_sd: f64,
    pub theta1_cov: Array2<f64>,
}

Conditional location-scale calibration of the latent score (#905).

The marginal-slope Auto trigger’s pooled-z gate (KS / skewness / kurtosis + the rank inverse-normal transform) only inspects the marginal law of z. A conditional shift E[z | C] = m(C) ≠ 0 — the allele-frequency-driven grouping mean shift — passes the marginal gate while leaving z | C off-center, so the slope contribution b(C)·m(C) leaks into the influence channel q. Rank-INT provably cannot fix this: no transform T depending only on the marginal F_Z can enforce E[T(Z) | C] ≡ const for all joint laws.

The unique Fisher-orthogonal location-scale correction (for the Gaussian working metric the closed-form probit kernel assumes) is ζ = (z − m(C)) / √v(C), where m(C) = E[z|C] and v(C) = Var(z|C) are estimated by weighted ridge regression of z (and its squared residual) on the marginal-index span a(C) = [1 | X_marginal]. The corrected ζ is conditionally centered (and homoskedastic when the variance block is active) by construction, so the b(C)·m(C) leakage vanishes. Matching the first two conditional moments does not by itself make ζ standard normal (a two-point residual law survives location-scale correction unchanged in shape), so build_latent_measure_with_geometry re-checks the calibrated sample against the standard-normal adequacy gate and retains an empirical latent measure for the residual distribution when that re-check fails; only a passing ζ uses the closed-form standard-normal kernel. Persisted so prediction rebuilds a(C) from the (reproducible) marginal design and applies the identical map to incoming z.

Fields

mean_coeffs: Vec<f64>

Coefficients for the conditional mean m(C) = β_m·[1 | a(C)] over the basis [1 | marginal-design row]. Length 1 + basis_ncols (leading entry is the intercept).

var_coeffs: Vec<f64>

Coefficients for the conditional variance v(C) = max(β_v·[1 | a(C)], var_floor). Length 1 + basis_ncols, or empty when the conditional-variance block of the Rao gate was not significant (mean-only correction); then v(C) ≡ homoskedastic_var.

basis_ncols: usize

Number of marginal-design columns in the basis (excludes the leading intercept). The predict-time marginal design must present exactly this many columns.

var_floor: f64

Floor on the fitted conditional variance, in the (normalized) latent-score scale (= AUTO_Z_CONDITIONAL_VAR_FLOOR_FRAC · the global weighted variance of the training score).

homoskedastic_var: f64

The homoskedastic conditional variance v(C) ≡ Var(z | C), used when the Breusch-Pagan stage did not fire and v is therefore constant in C.

This is the residual variance of the conditional-mean regression, Σ w (z − m̂(C))² / Σ w, and NOT the global (marginal) variance of z — which is what it used to hold, and what its old on-disk name still says (gam#2768).

The distinction is the whole correction. ζ = (z − m(C))/√v is supposed to be conditionally standard normal, because the marginal-slope identity that makes q the MARGINAL index, E_ζ[Φ(q√(1+b²) + bζ)] = Φ(q), holds only at Var(ζ|C) = 1; at Var(ζ|C) = v it becomes Φ(q√(1+b²)/√(1+b²v)), a multiplicative distortion of every marginal coefficient. Dividing by the marginal variance guaranteed v ≠ 1: with z standardised, 1 = Var(m(C)) + E[Var(z|C)], so the residual variance is 1 − R² and is strictly below the marginal variance whenever the gate fires at all. The bug was therefore not a corner case — it was on every fired gate, and it grew with exactly the conditional structure the correction exists to remove. At R² = 0.25 it left sd(ζ) = 0.87 against the post_sd ≈ 1 this struct documents, distorted the marginal coefficients by ~4%, and (worse) made the calibrated residual fail the standard-normal adequacy re-check on the SD clause alone at any appreciable n — sending BMS to the empirical measure and, per gam#2718, withholding the covariance.

Kept under the on-disk name global_var so a model saved before the fix still deserializes AND still applies the map it was fitted with: the stored number is whatever that fit divided by, and predict must reproduce the fit, not the current formula.

post_mean: f64

Weighted mean of the calibrated training sample (sanity-check, ≈ 0).

post_sd: f64

Weighted SD of the calibrated training sample (sanity-check, ≈ 1).

theta1_cov: Array2<f64>

Joint first-stage (generated-regressor) sandwich covariance of θ₁ = (mean_coeffs, var_coeffs), shape dim θ₁ × dim θ₁.

Replaced a stored PAIR of per-stage sandwiches that Self::theta1_covariance assembled block-diagonally – i.e. that asserted Cov(mean_coeffs, var_coeffs) = 0. The stages are not independent: stage B regresses the squared MEAN residual on the same basis, so the stacked bread is block lower-triangular and the meat has a cross-block proportional to the residual’s third moment. Both vanish under a Gaussian residual, and neither vanishes on the branch this covariance serves (gam#2484). Built by stacked_first_stage_sandwich_cov; the two retired fields were its diagonal blocks.

Fit-time only: predict applies the map from mean_coeffs/var_coeffs and never reads their uncertainty. Verified rather than assumed – the only production consumer of this matrix is the Murphy-Topel assembly in block_specs.rs, which runs at fit.

#[serde(default)] so a model saved before gam#2484 – carrying the two retired per-stage blocks and no joint – still deserializes on the predict path. An empty default is NOT a silent zero: it is refused at the single point of consumption (see Self::generated_regressor_correction), because a covariance that quietly vanishes is exactly the failure this issue exists to prevent.

Implementations

impl LatentZConditionalCalibration

pub fn apply( &self, z: ArrayView1<'_, f64>, a_block: ArrayView2<'_, f64>, ) -> Result<Array1<f64>, String>

Apply ζ = (z − m(C))/√v(C) to a batch. a_block is the marginal design (n × basis_ncols); z is the (normalized) latent score. Used at both training and predict time, so the map is identical.

pub fn theta1_dim(&self) -> usize

Dimension of the first-stage parameter vector θ₁ = (mean_coeffs, var_coeffs) whose estimation uncertainty the generated-regressor correction propagates. Equals len(mean_coeffs) when the variance block is inactive, otherwise len(mean_coeffs) + len(var_coeffs).

pub fn zeta_theta1_jacobian_row( &self, z: f64, a_row: ArrayView1<'_, f64>, ) -> Vec<f64>

Per-row sensitivity ∂ζ_i/∂θ₁ of the calibrated score to the first-stage calibration coefficients, stacked as [∂ζ/∂mean_coeffs | ∂ζ/∂var_coeffs] (length Self::theta1_dim). With ζ = (z − m(C))/√v(C), A_i = [1 | a(C_i)], m = A_iᵀ·mean_coeffs, v = A_iᵀ·var_coeffs:

∂ζ/∂m = −1/√v, ∂ζ/∂v = −(z − m)/(2 v^{3/2}) = −ζ/(2v),

and by the chain rule through the affine basis ∂ζ/∂mean_coeffs = (∂ζ/∂m)·A_i, ∂ζ/∂var_coeffs = (∂ζ/∂v)·A_i. The variance block contributes only when var_coeffs is active AND the fitted v(C_i) is above the floor (a floored row has ∂v/∂var_coeffs = 0 in the applied map). z is the (normalized) raw latent score at this row.

pub fn theta1_covariance(&self) -> Array2<f64>

Joint first-stage covariance V₁ of θ₁ = (mean_coeffs, var_coeffs) of θ₁, ordered to match Self::zeta_theta1_jacobian_row. The two stages are fit on (asymptotically) uncorrelated estimating equations (the mean score Σ w û A and the Breusch–Pagan variance score Σ w (û² − v) A are orthogonal under the Gaussian working model), so the joint first-stage covariance is block-diagonal to first order — the same approximation the Rao gate above uses.

pub fn generated_regressor_term( &self, hbeta_inv_g: ArrayView2<'_, f64>, ) -> Array2<f64>

Murphy–Topel generated-regressor correction term for the second-stage slope covariance. Given the second-stage information H_β (the penalized joint Hessian of the slope fit, whose inverse is the naive V_β) and the cross-derivative G = ∂(score_β)/∂θ₁ (p_β × dim θ₁), the corrected covariance is

V_β = V_β^naive + (H_β⁻¹ G) V₁ (H_β⁻¹ G)ᵀ.

This returns the additive rank-dim θ₁ term (H_β⁻¹ G) V₁ (H_β⁻¹ G)ᵀ given the already-formed hbeta_inv_g = H_β⁻¹ G (p_β × dim θ₁). The caller forms G by accumulating the per-row slope-score sensitivity to ζ_i times Self::zeta_theta1_jacobian_row (chain rule ∂score_β/∂θ₁ = Σ_i (∂score_β/∂ζ_i) (∂ζ_i/∂θ₁)).

pub fn generated_regressor_correction( &self, score_zeta_sensitivity: ArrayView2<'_, f64>, z: ArrayView1<'_, f64>, a_block: ArrayView2<'_, f64>, vb: ArrayView2<'_, f64>, ) -> Result<Array2<f64>, String>

Assemble the full Murphy–Topel generated-regressor correction (Vb·G)·V₁·(Vb·G)ᵀ for the second-stage slope covariance, given the ONE engine-side quantity it cannot reconstruct post-fit: the per-row reduced-frame slope-score sensitivity to the calibrated score, s_i = ∂score_β,i/∂ζ_i (a p_β-vector in the joint flat-β reduced frame solved_fit.beta_covariance() lives in). With score_β,i = ∂ℓ_i/∂β, s_i = ∂²ℓ_i/∂β∂ζ_i = J_iᵀ·(∂²ℓ_i/∂η_i∂ζ_i) is the mixed (β, ζ) second derivative of the warped row kernel contracted through the slope design Jacobian J_i — exactly the #932 RowProgram z-jet channel (z is already a row-program input; one extra mixed (β, z) jet channel reads off ∂²ℓ/∂β∂z). It must be evaluated at the converged β̂ in the SAME reduced frame as vb.

Everything else is built here from the stored first-stage quantities and the second-stage fit, dissolving the post-fit-reconstruction blocker:

• G = Σ_i s_i · (∂ζ_i/∂θ₁)ᵀ (p_β × dim θ₁), the chain-rule outer product accumulated row-by-row with ∂ζ_i/∂θ₁ = Self::zeta_theta1_jacobian_row(z_i, a_row_i) (exact-zero on floored rows, so floored rows contribute nothing — G’s support is the gate-fired rows);
• Vb·G = vb·G since the naive second-stage covariance vb IS H_β⁻¹ (the coordinator’s H_β⁻¹ G = Vb.dot(G));
• the term (Vb·G)·V₁·(Vb·G)ᵀ via Self::generated_regressor_term.

score_zeta_sensitivity is n × p_β (row i = s_i); z is the per-row normalized latent score (n); a_block is the marginal design n × basis_ncols whose rows feed zeta_theta1_jacobian_row; vb is the naive reduced-frame slope covariance n_β × n_β. The returned term is PSD (a congruence of the PSD V₁), so adding it to vb makes the corrected slope SE strictly ≥ the naive SE whenever the gate fires (G ≠ 0) and exactly equal when every row is floored (G = 0).

Trait Implementations

impl Clone for LatentZConditionalCalibration

fn clone(&self) -> LatentZConditionalCalibration

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for LatentZConditionalCalibration

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<'de> Deserialize<'de> for LatentZConditionalCalibration

fn deserialize<__D>(__deserializer: __D) -> Result<Self, __D::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl PartialEq for LatentZConditionalCalibration

fn eq(&self, other: &LatentZConditionalCalibration) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl Serialize for LatentZConditionalCalibration

fn serialize<__S>(&self, __serializer: __S) -> Result<__S::Ok, __S::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.

impl StructuralPartialEq for LatentZConditionalCalibration