Population inference (ch. 4): the vpop summary-likelihood fit - #88
Merged
Conversation
The Chapter 4 guide states the hierarchical VPop design: the three objects (anchored prior pi, proposal pi-tilde, inferred population (mu, omega)), the omega eigenbasis, the center-vs-spread target split (across_patient feeds mu+omega+reachability, center_only feeds mu only), the multi-round re-sim TSNPE loop, reporting under pi via the self-normalized reweight, and the full-path SBC gate. Wires it into the docs index (Ch. 4 row + status). vpop/proposal.py adds the generic spine piece: EigenbasisPopulation (the log-normal population/proposal in the spread eigenbasis, with an analytic log_prob that plugs into the importance reweight), widen_on_identified (build pi-tilde by widening sigma_u on the identified directions), and reachable_accept_fn (the data-anchored, observable-space truncation). 13 tests.
numpy's default quantile method ('linear', Hyndman-Fan type 7) shrinks the
sample IQR at the n published QSP targets carry: 0.788 of the population IQR
at n=6, 0.848 at n=9, 0.891 at n=12 (measured on a Gaussian population and
reproduced to two decimals on a 2000-patient PDAC cloud at each target's real
n, so it is a property of the estimator and n, not of any population shape).
Any estimator comparing that summary against a POPULATION quantile inherits
the shrinkage as a 15-25% LOW bias on sigma_u_hat at n=6-12, which is the
failure ch. 4's matched footing exists to prevent. Pin 'normal_unbiased' for
the observed anchor and the cohort summary alike (residual bias ~1%), as a
module constant rather than a per-call argument: the requirement is that both
sides use the SAME estimator, and a parameter would let them diverge.
Second fix: resolvability is per-anchor and does not make the anchor VECTOR
Gaussian. Nine deciles from eleven patients all pass the old filter; measured
D^2/k against the asymptotic covariance is 0.56 at n=6-8 and 0.75 at n=9-12
(1.00 is exact), and no bandwidth repairs it. _resolvable_grid now also
enforces n >= 2k, falling back to a symmetric 3-anchor subset with ties broken
outward, so a small-n target keeps spread anchors instead of collapsing to the
median. Three anchors measured calibrated at every spacing from (.25,.5,.75)
out to (.1,.5,.9).
Both constants carry their measurements in the docstrings.
… to take (ch. 4b) Ch. 4 keeps the model (two covariances, provenance split, matched footing, discrepancy layer, variance budget) and is demoted to the escalation path; its amortized NPE estimator is heavier than the data can exercise. Ch. 4b keeps that model, leads with the misspecification diagnosis the fixed-cloud ESS was asking for, and fits with an explicit summary likelihood sampled by NUTS. Ch. 4b's likelihood section is rewritten around two measured results: - The asymptotic quantile law holds at every n in the PDAC target set, but only in log space and only at three anchors below n=20. Escalation case 1 closes. The measured D^2/k table ships with the claim. - The law is now stated in its general copula form, because targets measured on the SAME patients have correlated sampling error and the per-target form was wrong for 33 of 49 PDAC observables (largest group: ten observables at n=10 from one GVAX+nivo arm). Factorizing over them counts ten correlated readouts as ten independent ones. eta_s does not cover it: a study offset is a shared systematic shift, this is shared sampling noise, and omitting the block inflates tau_eta_hat so finite-n noise reads as between-study heterogeneity. The model is restated per cohort_id. Also: differentiability recorded as a build requirement rather than a hope (NN emulator decided); the viability filter's discontinuity relocated to where it actually bites (variable cohort membership under sorting, not the normalizer) with the smooth-weight fix; build order gains the block-covariance validation as a gate on the fit; vpop/diagnostics.py corrected from "to build" to mostly built.
Ran the check the previous commit's build order calls for: cohorts drawn as rows of the cloud, so one study's targets are read off the same patients, and the measured study covariance compared against the copula formula. The formula is right. Across-observable entries reproduce to within a few percent (median measured/formula ratio 1.0 over eight multi-target studies). The magnitude claim in the previous commit was wrong and is corrected. It said factorizing counts ten correlated readouts as ten independent ones, up to an order of magnitude of overstated precision. Measured, the model's own cross-observable dependence is weak (median |rho| = 0.21 across the ten GVAX+nivo readouts, anchor-vector correlation 0.14 against 0.63 within an observable), so independence costs about four points of coverage, not an order of magnitude. What the dependence IS concentrated in is a handful of near-duplicate pairs: (cd8_pct, cd8gzmb_pct) at rho = 0.990 and (cd8_fc, cd8gzmb_fc) at rho = 0.994, with anchor correlations of 0.63 and 0.67. Those are the same pairs duplicate_observables already flags as a model-structure finding, so the covariance correction and the misspecification triage point at the same objects. Also corrected: the density plug-in is weaker than the earlier table showed. That table was measured against a smoothed population, which flatters the sparsity estimator. Against the raw cloud, h=0.035 makes the asserted variance ~18% too large on two heavy-tailed fold-change marginals while staying within 2% elsewhere, so the bandwidth wants a per-observable rule rather than one declared constant.
Checked the near-duplicate correlations the previous commit called a model-structure finding. They are not one. GZMB+ CD8 cells are a subset of CD8 cells, and the model's per-patient GZMB fraction does vary: log-sd 0.31 to 0.41. It is small next to the level, which varies with log-sd 2.4 to 2.6. With log X_gzmb = log X_cd8 + log r and the two roughly independent, rho = (1 + (sigma_r/sigma_L)^2)^-1/2, which is 0.99 at sigma_L = 2.55 and reproduces the measurement. So rho is manufactured by the width of the population it is evaluated on, and the cloud is a generous proposal, not the fitted population: log-sd 2.55 for cd8_pct_baseline against 1.03 in the published values. At the observed width the same decomposition gives rho ~ 0.93. Right way round for the implementation, since C_jl is evaluated at the current phi and so lands on the fitted population, but it makes the measured block sizes an upper bound. The screen does find one genuine duplicate: treg_fraction_cd4 and treg_fraction_cd4_hiraoka2006 are byte-identical columns, one quantity under two target names. The lesson for duplicate_observables is that a Spearman screen on levels conflates "same variable" with "same scale, different fraction". Also recorded: all 21 target rows carrying a values list store it SORTED, so per-patient pairing is gone. The covariance block is unaffected (it needs only the model's copula), but whether the model's cross-observable dependence is CORRECT is not answerable from the current target data.
The chapter specified the design but not the machinery, so a build starting from it would make three wrong assumptions. "Infer all P components of mu and log sigma_u" needed a hard qualifier. 227 of 271 parameters have a pool log-sd of 0.001 -- the simulations never varied them -- so only ~44 directions carry signal, and measured against the layered omega prior only 17% of the population's nominal log-variance sits on parameters the emulator has ever seen move. Along the rest the likelihood is not flat, it is undefined, and a prior-only report there is an artifact of the simulation design rather than of the data. Widening the pool is the fix; reporting them as inferred is not. The surrogate-error term was presented as a tidy variance contribution without saying how big it is. Measured on the new NN emulator (held-out R^2 0.984), the small-n targets where sigma_u is identified are sampling-limited (eps/se 0.56) and the larger-n targets are surrogate-limited (eps/se 2.6), so 4767 nominal patients carry the information of 835. Do not quote a target's n as its weight. Viability is 38%, not a rounding error, so Z(phi) is large and strongly phi-dependent and its saturation bites inside the range where sigma_u is read. Records that RestrictionClassifier.score already returns the unthresholded probability the smooth-weight fix needs, and that its gradient is piecewise-constant. Also records the within_cohort trap (any budget CSV predating the omega-layer fix is not comparable), and PDAC_IL10_tumor_tissue_Herremans2023 as a broken output extraction rather than a modelling problem.
Two places the simulation design was being reported as a modelling result. Part I listed five causes of a low prevalence-weighting ESS. It was missing the one currently measured to be present: the cloud is narrow by construction. vary_policy pins a non-varying parameter by driving its sigma to ~0, conflating omega_j = 0 with mu_j known (pdac-build 30967d1), and 227 of 271 parameters have a pool log-sd of 0.001, so the cloud carries 17% of the population variance the omega layer specifies. A tilt asked to reproduce the observed SPREAD from a cloud that concentrated has to strain and pays in ESS. It goes in as cause 3, ahead of support deficiency, because it is diagnostically distinct from it (a pinned-narrow cloud can still cover every bin, so support_deficient stays clean while the ESS collapses) and because it is one number to check. The end-to-end triage gains it as step 2, flagged as the exception in both directions: cheapest to check, most expensive to fix, since the fix is a re-simulation rather than an afternoon. The variance budget gains a fourth row, "not representable": the share of tr(Gamma_omega) on parameters the pool never varied. It was being accounted as "asserted", which reads as a prior choice. The two are not the same claim. Asserted means the data were silent along that direction, which is correct Bayesian behaviour and is what a prior is for. Not representable means we never asked, and no cloud, emulator, or fit built on that pool can say anything about it. At 83% it is the largest line in the budget and the most over-readable, so it ships as its own row with omega_supported_frac beside it. Neither is reachable from the estimator side, which is the point of recording them here rather than in a fitter's docstring.
The measurements behind targets/anchors.py:QUANTILE_METHOD and MIN_PATIENTS_PER_ANCHOR, and behind ch. 4b's copula covariance block, existed only as throwaway scripts. They need rerunning whenever the simulation pool changes, so they belong in the repo. Two checks. The diagonal one asks, per observable at that target's real n, whether the asymptotic quantile law describes the finite-n sampling spread (Mahalanobis D^2/k against the asserted covariance, plus the finite-n IQR bias that lands on sigma_u_hat). The block one draws cohorts as ROWS of the cloud so a study's targets share patients, and compares the copula covariance against the independence assumption. Parameterised via argparse with no project paths baked in: cloud npz, targets CSV, and the column/key names are all options. Reproduces the session results on the current PDAC cloud: d2/k of 0.97 at n=6-8 and 0.96 at n=9-12 on the quartile grid, IQR ratio 1.00 to 1.04, and the cross-target blocks including the two near-duplicate pairs at rho 0.99 where independence visibly degrades (d2k 0.54 copula vs 0.83 independent).
The closed-form alternative to the amortized estimator. The observed data are per-target quantile anchors at each target's real n; their joint asymptotic sampling law is available in closed form, so write it down and sample phi with NUTS instead of training a density estimator on a hyper-proposal. vpop/summary_likelihood.py Cov(Q_j(p_a), Q_l(p_b)) = [C_jl - p_a p_b] / (n f_j f_l), with the within-observable branch exact (min(p_a,p_b)) and the cross-observable branch a Gaussian copula. bvn_cdf is a 32-node Gauss-Legendre quadrature of the d/drho Phi_2 = phi_2 identity under the r=sin(t) substitution: differentiable in all three arguments, 8e-9 against scipy. StudyBlock stacks every target sharing a cohort_id, so 33 of 49 PDAC observables stop being counted as independent observations. sigma_c and epsilon stay diagonal (neither is a property of the patient sample). The mean model pushes a common-random-number population through a caller- supplied differentiable emulator and reads population quantiles off the resulting cloud. One forward, one sort, gathers into a precomputed query table. Two things measured while building it, both of which changed the design: freeze_covariance_at. Sigma scales with the population width (f ~ 1/sigma), so a live covariance's -0.5 logdet term rewards a narrow population on its own. Fitting real n-cohort quantiles from a known Gaussian population, mean sigma_hat/sigma_true: 0.878/0.941/0.963 live at n=20/40/120 against 1.002/0.987/0.979 frozen. Same low bias on the population spread, at the same published n, that QUANTILE_METHOD removes on the summary side. Freezing is the estimating-equation reading in which Sigma is a weight, and it is what makes ch. 4b's own claim about the density plug-in literally true. The plug-in's phi_0 barely matters, so there is no loop. On a curved forward map at n=20, a phi_0 twice too wide moves phi_hat by 0.14-0.32 of one sampling sd and leaves the spread unbiased. Refitting at phi_hat is a diagnostic, not a required pass. (Measured on a linear map first, which cannot see this: scaling sigma_u scales every f_j identically, so Sigma moves by a constant and the argmin does not move at all.) vpop/population_fit.py phi = [mu | log sigma_u | eta_raw | delta_raw | log tau_eta | log tau_delta], one unconstrained vector, offsets non-centered. All P components of mu and log sigma_u are inferred; the prior answers where the data are silent, and spread_movement is the readout that replaces K selection. Pyro NUTS against a potential_fn, chosen over NumPyro because the emulator is a torch module and the whole phi -> anchors map has to be one autodiff graph.
Measured against the real PDAC target set through the ch. 4b summary
likelihood: 16 of 65 observables produced non-finite anchors and dropped
silently out of the log-scale fit. All 16 were the ci95 branch, and the cause
was not sign but a missing center. They carry target_value = NaN with a
perfectly good positive interval, so the Gaussian fallback computed
NaN + z*sigma. 14 are qualitative mechanistic constraints ("tumour does not
regress: fold change 1 to 5"), which state a RANGE and no point estimate.
State the branches by what a target actually supplies rather than requiring all
three inputs:
center + both bounds positive -> two-sided log-sd about the median (unchanged)
no center, both bounds positive -> geometric midpoint sqrt(lo*hi); a 95%
interval implies its own center under the same distribution used to
expand it, and the recovered 2.5/97.5 anchors reproduce the bounds
center positive, one bound usable -> one-sided log-sd. A CI reported as
touching zero was otherwise thrown onto the Gaussian branch, where it
returns negative quantiles for a strictly positive quantity
A signed quantity keeps the linear Gaussian and may legitimately have negative
anchors; those cannot feed a log-scale likelihood and the caller must exclude
them, which is now said in the docstring rather than discovered.
On the PDAC set this takes the ch. 4b fit from 49 usable targets in 24 study
blocks (147 anchor rows) to 65 in 27 blocks (163 rows).
States the model the vpop path implements: log-normal patients with a fixed correlation matrix and free marginal widths, cohort quantiles as the data, and an anchor covariance held fixed during the fit. Every term carries a status, so the six inferred quantities are separable from the inputs at a glance. Records five things the code does not say out loud: - Corr(Omega) = R is an input taken from stage 1, not a fitted quantity, and the copula behind it covers 65 of 271 parameters and states independence elsewhere - eta and delta are exactly aliased, because no observable appears in two cohorts, so 163 anchor rows overstate the information by a wide margin - kappa, a per-observable width correction, is what stops a narrow simulator being reported as diverse patients; without it only omega can absorb the gap - the emulator error slot exists in the likelihood and nothing populates it - the n_biological gate on the omega prior is load-bearing: without it, 145 epistemic posterior widths would enter as between-patient spread Also records why the eigenbasis was dropped. It rotated the spread into directions the measurements resolve, which moved Corr(Omega) by a rank-K update and left sigma_u with no per-parameter meaning.
The observation block is the cohort, not the target. F is joint over the readouts a cohort reports, the bootstrap resamples whole patients, and V_c carries the dependence between readouts of one cohort instead of assuming it away. Independence across cohorts replaces it as the standing assumption. Also: E_c defined on the statistics rather than the readouts, with its diagonal fallback and what that asserts; a flat fit section that holds omega at omega_0 rather than zero, so the bootstrap and the eligibility rule keep working and the scale rows become a held-out width check; and reported uncertainties entering through a D_c rho_c D_c split, with a rule for whether a lone reported uncertainty becomes a D_c entry or a scale statistic.
Adds examples/qsp_scale_mechanism.py, a synthetic 179-species / 271-parameter sparse reaction network sized from cpp/qsp/ode/QSP_enum.h and parameters/pdac_priors.csv, with 5 nested trial arms and readouts in the kinds vpop_marginal_targets.csv reports. examples/toy_population_fit.py now runs against either it or the 3-species toy via --size small|medium|full. Model: - tau_beta fixed at 0.15 rather than sampled. The hyperprior sat over 2 free directions, bought no pooling, and was the largest source of divergences (24 -> 0 in the flat fit). - The eq:disc affine map pivots about each study's own level, so kappa is a spread term within a cohort. Its lever on the location rows falls from 84 sd to 0.7 sd, and b goes from a 2.2 sd miss to within 1 sd. This makes the measurement bias study-dependent, at no cost in parameters. - The flat fit is an optimiser plus a Gauss-Newton covariance rather than NUTS, and pins b. On noiseless data the flat fit recovers b exactly, but with real rows its only channel is the mean-median gap and 19 parameters against 25 location rows leaves it 3 sd off. Pinning costs nothing measurable across 6 data seeds and removes a nuisance mu_hat had to trade against. Sampler: - NUTS is seeded with G^-1 from the Gauss-Newton information at the plug-in instead of adapting a dense mass matrix, which at 210+ coordinates it cannot estimate from warmup. Trees terminate instead of saturating at 1023 steps; the small problem's two fits go from 12.5 min to 9 s. Diagnostics: - conditioning_report: posterior anisotropy before any MCMC. Shows the a block carrying 95% of the Fisher information. - phi0_sensitivity: rebuilds V_c across Sigma_1 to test whether the plug-in matters at all. - print_pivot_offsets: kappa's residual location lever, per row. - Build-time warnings for singular V_c blocks and for |dbar| exceeding a cohort's sampling noise. - Emulator loss weighted by measured readout sensitivity, so capacity is not spent on species no readout reads.
Adds examples/run_toy_population_full.sbatch and a --no-progress flag so the NUTS progress bar does not fill a scheduler log with carriage returns. The job does not install the package. Both example modules import only jax, numpyro, diffrax, numpy and scipy, while qsp-inference itself depends on maple over git, which a compute node generally cannot fetch. uv supplies its own Python so the job does not depend on the cluster's default interpreter. No account or path is baked in: the allocation is passed to sbatch and the job runs from SLURM_SUBMIT_DIR.
Under a scheduler the progress bar goes to stderr and stdout is a file, so the log was silent through the ~18 minutes of setup before NUTS starts, which is indistinguishable from a hung job. Phase markers carry a wall clock and elapsed minutes, so each stage can be costed from the log alone. stderr now goes to its own file, which means the progress bar can be kept rather than silenced: .out is readable prose, .err carries a live ETA.
Both are pure functions of the design and the seed and neither depends on the data, yet at QSP scale they are the bulk of the wall clock: the emulator is ~102,000 ODE solves plus 10,000 Adam steps, and eq:Ec is 225,000 solves per build, done at least twice a run. Keyed by a content hash over the exact array bytes of every input, so a hit means the inputs are identical rather than similar. The V_c key includes a hash of the trained emulator weights themselves rather than the config that produced them, so a cached V_c can never be paired with a different surrogate. --no-cache forces a rebuild; --cache-dir moves it.
The XLA CPU backend does not automatically thread a vmapped diffrax solve, so the emulator build ran at 4.4 cores of 24 and took 15.5 minutes on a batch node. Host devices are the only axis JAX parallelises well on CPU, so the true-ODE cloud is now pmapped across them, with padding for batches that do not divide evenly and a vmap fallback for small ones. Verified bit-identical to vmap. Only the true ODE is wrapped. The emulator's cloud is called inside the NUTS gradient, where pmap overhead and nested tracing would both cost more than they save. The job script now asks for one host device per core and one thread per device, and runs 4 chains so that sampling occupies 4 devices rather than 2.
…endezvous At 32 devices on 32 cores the emulator build died with an all-reduce rendezvous timeout: 'Expected 32 threads to join the rendezvous, but only 29 arrived on time'. XLA's CPU collectives need every device thread co-scheduled inside 40 seconds, which one device per core does not reliably achieve. Eight devices leaves roughly four cores of headroom each. Throughput is lower than the theoretical maximum, but the emulator and V_c are cached, so their build is a one-time cost.
print_width_shortfall still reached for samples['b'] unconditionally. Pinning b removed that key, so the full-size run died with KeyError: 'b' immediately after printing the flat recovery table. summarise_recovery and print_readout_effects were guarded; this consumer was missed. With b pinned the shortfall is predicted at kappa = 1, which is exactly what the flat fit assumed.
JAX dispatches asynchronously, so mcmc.run returns once the work is queued. The full-size run printed 'NUTS done in 17.5s' while the chains were still at 0/1000 and sampled for a long time afterwards, which made a running job look finished and a phase marker report a completion that had not happened. Blocks on get_samples before stopping the clock.
The mass matrix was built from the Gauss-Newton information at the prior centre, where a = b = 0. That is not where the sampler lives: a ends up pinned to about 1/849 of its prior width. At full size the symptom was a straggler -- three chains finished 1000 iterations in ~70 minutes at ~4 s/it while the fourth managed 150 at 31 s/it, and numpyro writes nothing until every chain finishes. A frozen metric offers no per-chain recovery, so a chain that lands where the plug-in curvature fits badly stays badly conditioned. population_map locates the mode by Adam and the metric is linearised there. print_dbar_absorption tests the draft's claim that gamma_r absorbs the emulator's systematic offset. gamma = Z a is indexed by readout while dbar is per row, so only the component in the column space can be taken. Measured on the small problem, 0.39 to 0.41 of dbar lies in span(gamma, kappa) and the remaining 0.92 is absorbed by nothing and lands on mu. The claim does not hold.
Both full-size runs produced exactly one chain doing 24-35 s/it while its siblings did 4-8, and a different chain each time. Linearising the metric at the mode tightened the fast chains from an 11x spread to 1.13x but did not remove the straggler, which points at initialisation rather than curvature: init_to_median disperses chains across a 572-dimensional prior, and a frozen metric tuned for the mode does not fit where a distant chain starts. The MAP is already computed for the metric, so reusing it as the init is free. Cost: chains no longer start over-dispersed, so r_hat loses power to detect a missed mode. Accepted deliberately -- numpyro writes nothing until every chain finishes, so one straggler currently costs the whole run.
A source that prints one uncertainty can tell you how precisely it pinned its own statistic, or how spread out its patients were. Not both. It now always goes to width, because with M empty the scale rows are all omega has, while sampling noise comes free from the bootstrap on every row. So U leaves V_c and becomes an "se" row of T_hat: the log sampling standard error of a median, predicted at phi. eq:Vsplit collapses to V_c = V_boot + E_c, since scales and correlations now both come from the bootstrap and the old split reconstructed V_boot exactly. Estimator is hybrid, closed form above n_c = 20 and a frozen-design bootstrap below it, where the asymptotics are poor. Cost runs with n_c, so bootstrapping only the small cohorts measured 1.12x the gradient at QSP scale against 2.7x for bootstrapping everywhere. The resample design is stored as uniforms rather than indices because the flat fit runs a smaller cloud, and an out-of-range JAX gather clamps rather than fails. The same estimator generates the data and predicts it, so the toy tests the mechanism rather than estimator mismatch. V_c now rests entirely on the predicted cloud, so print_width_gate runs before the fit. It reports the width shortfall at phi_0 against the prior sd of s + b_1, which is the reachable question; a bare shortfall is expected whenever s* > 0 and says nothing on its own. Drops the "synthesized" concept entirely. Reconstruction error is a different quantity from a reported spread and belongs in V_c as an additive term, not as a scale override; it is not built. V_c's cache key gains the cohort stats. Keyed on names alone, adding an se row would have silently paired the new table with a stale V_c. Measured on the small problem: two se rows tighten (s+b_1)/sqrt2 by about 6% and leave (s-b_1)/sqrt2 untouched, which is the predicted behaviour, since kappa scales a standard error exactly as it scales an IQR.
The old test asked whether the mechanism CAN reproduce a column of Z's row effect, by projecting onto col(d tau / d(mu, omega)) and reading the residual. Once the mechanism spans row space that is answered yes by arithmetic: at P >> K the residual is zero for every column whatever Z is, so passing it says nothing. At full scale it returned 0.000 across the board. Ask the price instead. Among mechanism moves reproducing the effect, take the smallest in prior units, delta* = argmin delta^T Sigma^-1 delta s.t. J_mech delta = W e_j, which in prior-standardised coordinates is the minimum-norm solution, so ||delta*|| reads directly in prior standard deviations. Large means only an implausible excursion imitates that column and it is identified against the mechanism; small means confounded. The target is a one-prior-sd move on the coefficient, so the ratio is scale free. Reports the unreachable fraction beside it. Before saturation the mechanism cannot span row space and part of a column is unreachable outright, which is a stronger verdict than an expensive one. Effective rank and the retained singular value range are printed too, since the minimum-norm solution loads on the weakest retained direction and the cost should be readable as a consequence rather than taken on trust. conditioning_report now hands back its standardised blocks so this reuses that Jacobian. At QSP scale it is a forward-mode pass over 572 coordinates and the slowest step in the script.
Cohorts that share patients were asserted independent, with merging them named as the fix. Merging is not generally available: a trial contributing 16 evaluable baseline biopsies against day-21 arms of 9 and 10 has neither a containment nor a common n. So eq:obs is now block diagonal by block, where a block is the connected component of cohorts joined either by a row built from both or by shared people. The two relations need opposite draws, disjoint arms resampled independently and shared patients resampled once from the block's set, so which one joined a pair is kept and not just the partition. V, E and the flat likelihood follow the same index; a one-cohort block recovers what was written before. What remains asserted is that no two publications overlap. The within-block correlation is predicted at phi_0 rather than measured, since a study printing order statistics at two timepoints does not print their joint, and where per-patient values survive digitisation V_boot can come from the reported patients instead.
resampling.py turns covariance_blocks' partition into a draw recipe. A counted PatientBlock becomes one joint draw whose strata say which patients each member cohort sees; every other cohort draws alone. Uncounted blocks cannot be drawn and are reported instead, so a caller cannot take V without the caveat. covariance.py reindexes the toy's bootstrap_V / emulator_E / assemble_V / subset_V from cohort to block. bootstrap_V takes rows_fn(cohort_id, indices) so the cloud and the readouts stay with the caller. On the PDAC registry: 24 blocks over 29 cohorts, li2022 drawing 19 patients jointly at the declared 6/10/0 overlaps. Median rows off that draw correlate 0.43 between baseline and arm A (6 shared), 0.74 between baseline and arm B (10 shared, arm B inside baseline) and 0.007 between the arms (none). emulator_E takes precomputed surrogate-minus-simulator differences. At QSP scale those must come from held-out pool rows, so it waits on the new emulator.
eq:smoothq in the mass form: Beta(kappa, n-kappa+1) integrated over each cloud member's slice of [0,1], which is exact and sums to one without normalising, rather than the density sampled at the midpoint. tau is E[what the cohort printed], not the population quantile. Verified against simulated 16-samples for normal, Gumbel and t(3) clouds at p = 0.25/0.5/0.75, so the shape-independence of the order-statistic law is exercised rather than assumed. The order statistic an estimator picks is a convention nobody records: type 7 (R, numpy) targets rank 0.30 for a lower quartile at n=9 where type 6 (SPSS) targets 0.25, some 0.4 sampling SE apart. Declared per call, defaulting to type 7.
Scaling column j of L_R' by omega_j commutes with the row sum, so z (L_R' * omega) is (z L_R') * omega and the matmul is phi-free. It was being recomputed every gradient: (N,P) x (P,P) at P=271. Cuts 2.59e10 flops to 2.00e10 but only 490 ms to 468 ms, because the step is bound by traffic over the (N,P) array rather than by the arithmetic on it, and that traffic is unchanged. Kept for being exact and free, not for the wall clock.
Data drawn at phi* are recoverable in principle: z is frozen and shared between the draw and the fit, so tau's Monte Carlo displacement cancels. Adam rather than a quasi-Newton method, since most directions are flat and a curvature step along them is fit to rounding. The potential comes from initialize_model, so the latent sites and their Jacobians are numpyro's business and not this module's. On the corpus at N=2000 it converges over 550 parameters and returns mu_raw and u_raw to the prior mode, which is what the 83xP Jacobian says should happen. a is recovered. s and b_1 come back with equal and opposite errors, but the likelihood alone is flat along s and steep against b_1, so that is the s prior shrinking a 2.2-sd draw and b_1 absorbing it, not the aliased sum it resembles. eq:betaprior centres beta, so a single declared species gives beta = 0 identically: a parameter with no effect, and a flat direction under NUTS. Refused alongside the every-width-measured case.
The pool was iid. Sobol goes in the copula's independent latent space, before the Cholesky, which is the only place that leaves the marginals exact; on the live 271-parameter composite prior it takes the mean per-coordinate chi2 from 29.4 to 3.1. The sampler is in the pool fingerprint because it decides which rows come out, with 'iid' contributing nothing so pools cached before the field existed are not orphaned to describe a draw that has not changed. vpop/emulator.py ports the trained net to JAX: NUTS needs gradients of the whole phi to cohort-quantiles map, and a torch module cannot supply them. Dropout is identity at eval, so the port is the weight matrices and two normalisations, and check_against_torch closes it because a misread layer stride returns plausible numbers.
Two bits of wiring the primitives do not cover. row_specs, block_draw_plan, quantile_mass_table and build_h_fn already exist, and composing them is a few lines at the call site, so they stay there. The scenario table does not, because reference-vs-readout ordering is where the fold-change bug lived. Keyed by (arm, time) rather than by cohort: cohorts sharing an arm read the same trajectory and are separated by eq:disc and the weighting, not by the mechanism. On the live corpus that is 7 scenarios over 5 arms rather than 30, and each of the four fold changes takes its reference from its own arm's day 0. prior_cholesky derives L_R from the PriorSpec the pool was drawn from rather than accepting it alongside. The surrogate learns the manifold that correlation induces, so a Mechanism carrying a different one evaluates it off-manifold with every number still plausible.
tumor_doubling_time reads each patient's first and last timepoint, and the last is patient-dependent, so no fixed readout time expresses it and h_r cannot compose it. The reduce computes it where the trajectory still exists; this is how it reaches h_r. Not through g_fn's block. That block is one fixed width per scenario, so a column only one arm emits would have to be padded everywhere else, and padding costs the other arms their training rows: emitting it for all five dropped 21% of the urelumab arm, selected for being responders, since a doubling time is undefined when the tumour does not grow. Those are the patients the fit is about. So it travels beside the block. build_extra_fn resolves each (scenario, name) against the arm that serves it and refuses one that does not, h_fn declares which pairs it needs, and readout_cloud evaluates exactly those. It is not scaled by beta, which is right for a log-ratio of one species at two times: a per-species bias cancels identically and beta could never have reached it. The third argument is passed only when there is one, so an h_r with no precomposed readouts keeps the signature it had.
min and max are order statistics 1 and n, so the Beta kernel already had them; they were unsupported only because nothing could consume them. Fitting them gives stromal_fraction its only width, at the cost of a row reading the 0.993 quantile. row_specs gains exclude, for a printed number that is not evidence: the corpus carries a ci95 pair built by assuming CV = 1.0, and fitting it would return the assumption to omega as if a study had measured it. It raises on an exclusion matching nothing, so a later corpus edit cannot quietly re-admit one. restrict_plans drops cohorts reporting no row from a block's row order while leaving its strata alone: a dropped member of a counted block still occupies its place, and removing it there would change the overlap its co-members draw.
h_r returns logs and eq:disc acts there, which is what makes gamma a multiplicative assay bias. The row functionals were inheriting that scale, so the likelihood compared a log-scale prediction against a natural-scale T-hat. Mean, sd, se and iqr do not commute with exp, so the exp belongs between eq:disc and the functional rather than applied to tau afterwards. The width gate found it on its first run: mean shortfall went from +1.319 log units (2.26 prior sd, unreachable) to -0.223 (0.38 prior sd), and max |z| from 122 to 7.7. Location rows go from a median T-hat/tau of 13.5, computable for 48 of 121 rows because a log is often negative, to 0.851 over all 121. kappa now leaves only the POPULATION median at the pivot. A reported median is an expectation over n draws, and in the source's units the map is a power rather than affine, so the expectation stops commuting with it. The gap is the median's own sampling spread and closes with n.
Splits tau_from_readouts out of tau_all so the same rows can be computed from a cloud the simulator produced, where there is no phi to push through the emulator. Both sides read the same resampled patients, so the cloud's Monte Carlo fluctuation is common to the two and differences out, leaving the surrogate's own error carried to the reported statistics. Measuring it on the readouts and expanding would not do: a quantile row and a spread row do not feel a per-readout error alike.
pin_discrepancy, pin_u and pin_aux make the model falsifiable. With 18 free discrepancy parameters a good fit says little, since most of a mismatch can be absorbed; with them off it has nowhere to hide and shows as residual structure that can be read and attributed. pin_u is declared rather than defaulted either way: pinning it asserts the widths are proportional, leaving it free adds 271 parameters that no scale row informs. Neither is neutral. The cost is that real assay bias then lands on mu, so this diagnoses and does not ship.
NUMPY_QUANTILE_METHOD took five Hyndman-Fan names and QUANTILE_CONVENTIONS implemented three, so a target declaring type2 or type8 loaded clean and raised KeyError inside the sampler. Nothing declares one today, which is why it never fired. type8 is continuous. type2 is not, but the caller already mixes the two order statistics around h, so the averaging case is h = np + 1/2 at integer np and ceil(np) otherwise. Asserts the two registries are equal so they cannot drift again, and pins h against hand-checked Hyndman-Fan values.
eq:Ec split into dbar and E_B, and eq:gcorr added: an affine alpha + Gamma (theta - mu_0) ridge-fitted over n_g solves at phi_0, in the transform the net was trained in, with g_c replacing g throughout including inside eq:Ec. E_B alone understates V_B by whatever dbar is, and nothing in the fit registers the difference. The two pieces do not scale together: narrowing the cloud shrinks V^boot and the fluctuation but leaves dbar where it was, so the emulator's share grows as omega falls. Gamma also enters d tau / d mu, so the correction reaches the Jacobian and not only the level. Freezing alpha is the riskier freeze: a stale V_B gives a row the wrong width, which eq:ratio tests, while a stale alpha moves tau and is paid for in mu.
population-generative-model.tex carried 1048 uncommitted lines and generative_model.md / data-and-model-sketch.* were never tracked. Committing the working state so the next commit's deletion is recoverable rather than final.
The rival drafts described the same model with an older discrepancy layer (estimated tau_eta / tau_delta / tau_kappa against the live eq:abprior's fixed sigma_a = sigma_b = 0.5), so reading the wrong one made the implementation look non-compliant. Recoverable from the previous commit. Four sections have no model-draft counterpart and are worth harvesting if wanted: what the emulator must cover, how omega_0 is built, the correction alias, and what data to gather next.
A sigma with no resolvable n leaves the prior untouched, but the level was set from whether a sigma was offered rather than from whether it changed anything. In pdac that is 178 offered against 9 resolvable, so the provenance reported a whole prior as data-informed when the data never reached it. Values are unchanged: shrink_toward_prior already returns prior_omega when n is None. This is the label only.
A cut parameter is left unvisited by the component BFS so several components can claim it, and only the upstream one is meant to reach the priors. The aggregation resolved that by dict order instead, which took the discarded downstream view for 3 of the 7 cuts. TGFb_50_CAF_act was carrying the NPE fit submodel_config.yaml cuts for SBC bias. Resolution keys on the component's own target list, from the cache's target_yamls, not on the union over the params it carries: a merged component holds targets belonging to its other parameters and matches an upstream it never used. Passing cascade_cuts without that list now raises rather than silently passing everything through. run_comparison no longer discovers submodel_config.yaml for itself. It sets the component partition through cascade_cut_params, so a path the function guessed could differ from the caller's and re-partition a run with no symptom. Test fixtures moved off MomentSpread and the moments dict, which predate the flat-statistics schema; two of them had been asserting against None.
Marginals were selected from {normal, lognormal, gamma, invgamma} by AIC, and
then load_composite_prior_log sampled the winner and refitted a normal to its
logs. Every marginal in the composite prior was therefore Gaussian whatever
stage 1 found: the prior was a multivariate lognormal by construction, and a
copula exists precisely to keep the marginals and the dependence separate.
The family also had no say in whether the result was usable. AIC picked an
inverse gamma of shape 0.906 for IFNg_50_CXCL9, which has no finite mean, and
1.965 for k_APC_mig_ID_LN, no finite variance. Both reached the prior.
Marginals are now the posterior's own quantiles on log theta, indexed by normal
score so that continuing linearly past the grid is exactly a Gaussian tail. The
tail slope is regressed over the outer band rather than read off the endpoint
derivative, which is the noisiest knot: on a known lognormal that gave 0.54
against a true 0.80, 3.4 nats of error in the upper tail.
On this corpus 61% of the fitted marginals failed Anderson-Darling and no
candidate family passed for 113 of them. Against the samples, KS falls from a
median of 0.010 to 0.005, and from 51 of 189 beyond the 5% critical value to
none. sd_1 moves by at most 6%, so the widths downstream are unaffected; what
changes is shape, tails, and log_prob.
The copula is fitted from rank pseudo-observations. Through the fitted CDFs a
poor marginal leaves u non-uniform and the correlation absorbs that error as if
it were dependence.
_all_normal is gone. It was an automatic fork: all-normal marginals took an
exact standardization, anything else sent every marginal through cdf -> Phi^-1
with u clamped, capping |z| at 5.61. One non-normal marginal silently degraded
the other 270, and which log_prob ran was decided by the data. Exactness is now
per marginal, and on the real prior all 271 reach z directly.
Grids live in a .npz beside the YAML. Inline they made the file 12x larger and
5x slower to parse, and every parameter repeated the same z. The YAML carries a
sha256 of the sidecar, checked on load, because the two are one artifact in two
files and nothing else stops them drifting apart.
…them Twelve fields carried defaults, so a caller could assert tau_s, sigma_a, n_beta or pin_discrepancy by not mentioning them. A rubric number asserted by a dataclass is indistinguishable in the posterior from one somebody chose, and five of the pdac fit's numbers were entering exactly that way. pin_discrepancy=False is not the neutral setting it looked like: it turns the discrepancy layer on, which is a modelling decision, so it is stated too. tau_omega_measured is the one exception. It decides nothing while `measured` is empty, so it is required exactly when a width is declared measured, and __post_init__ raises rather than supplying a width for a claim about how far a measurement is trusted.
comparison.py built each posterior-predictive draw from the MCMC samples and then overwrote every declared-nuisance parameter with a fresh prior draw. An observable depending only on nuisance parameters therefore had a prior predictive reported as its posterior predictive, understating coverage on 24 of the 43 components that looked broken. ppc_audit reads that block back as per-observable evidence and attaches no verdict. Thresholds for "badly fitting" did not survive testing: the binary reachability flag rested on min/max of 200 draws, and 9 of 42 flagged observables sat inside the noise band of that order statistic. What did hold up is the sensitivity ratio, prior predictive width over CSV prior width, which separates an observable theta cannot move from a merely tight prior. refit_check decides instead, by fitting a target set with and without a change over an identical isolated directory. resolve_target_set carries the closure the stage DAG demands: it walks the entire cascade cut list and rejects any upstream it cannot place, whether or not that cut's parameter is in the run.
pin_u let a caller assert the answer twice: tau_u already says how far u may move, and a flag on top of it made zero and free the only choices. Wide prints a width profile whose individuation is a prior draw; pinned stops a direction the scale rows do constrain from moving at all. Small and free is neither, and needs no rank cutoff defended. The comment at the sample site carries what the posterior cannot: that tau_u is derived from the role table's narrowest gap so u cannot carry a parameter out of the role it was assigned, that the normal prior shrinks uniformly and that is a choice against a heavy tail, and that a component sitting at tau_u is the expected reading rather than a finding.
omega said how wide a parameter varies between patients and implicitly asserted that it varies lognormally. For a quantity bounded on (0,1) that assertion is false at every width: Emax_Cy_Treg at median 0.80 puts 26% of the population above complete inhibition at omega 0.35, and the only omega that fixes it is one denying the quantity varies at all. The overrides CSV gains a `scale` column, orthogonal to the role: the role says how wide, the scale says on what. eq:crn is a Gaussian copula, so the margin is separable from the correlation and only patient_cloud changes. mu stays on the log scale for every parameter, so eq:muprior and the stage-1 copula need no reinterpretation and z = 0 still returns exp(mu) exactly. Logit-normal rather than Beta because omega keeps its meaning as a scale on an unbounded latent, so the omega_0 role table, s, u and the tau_u derivation all carry over. Beta's dispersion is not a log-sd and its median has no closed form. The guard in patient_cloud is load-bearing. Clipping exp(mu) at the bound keeps the value finite but leaves 1/(1 - m) ~ 1e8 in the branch jnp.where discards, and that is enough to make d/dmu nan; the base has to be masked before the nonlinearity, not clamped at the singularity. load_omega_overrides also now refuses a repeated parameter name. It returns a dict, so a second row for the same name used to win silently and drop whatever rationale the first one carried.
The gate read the rows SCALE_STATS names, so a source printing an iqr was measured and a source printing q25 and q75 was not, though both report the same spread of the same patients. That made the width budget a sample of editorial style rather than of the corpus. paired_width_rows differences a target's widest symmetric quantile pair into one width and enters it beside the printed se, sd and iqr rows. One row per target, so q25/q50/q75 is not read as three overlapping widths; a target that printed an iqr outright is skipped, since the two would be the same number twice. The derived row needs no new data and no new assumption: both quantiles belong to one cohort, so their contrast variance is already exact in V_B, and the delta method carries it to the log scale the printed rows use. Non-monotone quantiles raise rather than pass, separately for the observed pair (a corpus error) and the predicted pair (a bug, since expected order statistics of one sorted cloud are monotone in p).
…ract main's PR #86 added the test with a bare gamma marginal and cross-checked that it still admits a log fit. This branch changed what that means: since the prior carries the posterior's own shape, a gamma without an empirical_log grid is not a log-space marginal and copula_prior refuses it by name rather than silently refitting it to a family. So the cross-check was asserting a contract that no longer exists. The rounding assertions are the subject of the test and are untouched; the downstream half now guards that a tiny scale reaches the format check as a live value rather than a zero, which is the regression the test was written for.
copula_prior and inference.importance import torch at module scope, and torch is the optional [sbi] extra. Both test modules imported them at module scope too, so collection failed outright in an environment without it rather than skipping, and the whole run aborted. Same guard main already uses for the copula_prior cross-check.
…y arithmetic order_statistic_mass built its edges as cumsum(w) / sum(w). Those are cumulative probabilities and belong in [0, 1], but cumsum and sum reduce in different orders, so the last edge can land an ulp above 1 -- and betainc is nan just outside the interval rather than saturating at it. Whether the ulp appears depends on the platform's summation order, which is why this passed on arm64 and returned nan on x86 CI, taking the whole Beta kernel with it. Clipping the edges is the fix rather than nudging the endpoints, since the bound is what the quantity means and not a tolerance.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The population-inference stack for chapter 4: the random-effects target
theta_i ~ F(theta | mu, omega)fit against cohort summaries, plus the guide that documents it.Nearly all of this is new surface —
vpop/predict.py,fit.py,rows.py,gates.pyand their tests,docs/,examples/— and cannot regress anything main ships. The commits that do touch existing modules were split out into #87 and should land first.Main pieces
eq:postas a NumPyro model:mu,s/u_raw,a/b,beta_raw,log_R, and nothing else.n_c-sample.V_Bper block, resampled by the registry's blocks; Beta order-statistic kernel for reported quantiles.mustays on the log scale, so eq:muprior is untouched.omega's patternuis shrunk by a prior rather than pinned by a flag.Note for review: merging main in surfaced one semantic conflict. PR #86's tiny-scale rounding test cross-checked that a bare gamma marginal still admits a log fit; on this branch a gamma without an
empirical_loggrid is refused by name instead of being silently refitted. The rounding assertions are untouched and the downstream half now asserts the current contract — worth a look, since it is a real behaviour change rather than a test fix.897 unit tests pass.