bayesmith.evidence
Source-derived signatures and contracts. See the reference index for coverage and conventions.
Deprecated compatibility namespace. Use bayesmith.marginal. This namespace must not be reused for graph-level EvidenceTask semantics.
Module contract
Deprecated alias for :mod:`bayesmith.marginal`. Retires at 1.0.
This subpackage was called ``evidence`` through 0.4.0. The name claimed
something the package does not do: the Bayesian evidence is
``p(d) = INT p(theta) PROD_i L_i(theta) d theta``, a number obtained by
integrating the parameters OUT, and nothing here computes or consumes one.
The package computes it elsewhere: an ``EvidenceTask`` answered by
:mod:`bayesmith.dispatch.evidence` for a graph that is one exact
linear-Gaussian block, and compared across two results by
:mod:`bayesmith.evaluation.evidence`. What the modules under here actually
build is each dataset's MARGINAL LIKELIHOOD
``L_i(theta)``, with its own nuisances integrated away: a function of the
parameters, not a number. Hence ``marginal``.
The rename is not cosmetic bookkeeping. ``evidence`` is a word this package
uses correctly in fifteen other places, in its ordinary English sense of
"the grounds for a verdict" -- including in ``errors.py``, which is the only
eagerly imported module and so the first bayesmith prose a user reads. A
subpackage holding the term made every one of those ambiguous, and no
disclaimer written inside the subpackage can reach the prose outside it.
**Everything still works through this path**, including deep imports::
from bayesmith.evidence import SqrtInfo # fine, warns
from bayesmith.evidence.compress import compress # fine, warns
import bayesmith.evidence.chain # fine, warns
Deep imports are aliased through ``sys.modules`` rather than left to a
module-level ``__getattr__``, because a ``__getattr__`` shim does NOT support
them -- ``from pkg.old.kernel import helper`` raises ``ModuleNotFoundError``
against one. That distinction is load-bearing here: of the seventeen names
0.4.0 published from this subpackage, the intersection with bayesmith's
top-level ``__all__`` was EMPTY, so every one of them was reachable only by
the deep path this file has to keep working.SqrtInfo
Defined in bayesmith.marginal.sqrtinfo, line 56.
class SqrtInfo(eqx.Module)
factor: jax.Array
target: jax.Array
offset: jax.Array
names: tuple[str, ...] = eqx.field(static=True)
shapes: tuple[tuple[int, ...], ...] = eqx.field(static=True)``log L(x) = -0.5 ||R x - z||^2 + offset`` over a flat, named vector.
Attributes:
factor: ``(r, n)`` array ``R``. ``r < n`` means the term constrains
only an ``r``-dimensional subspace -- the normal case for a single
epoch, not an error.
target: ``(r,)`` array ``z``.
offset: scalar; every part of the log-density that does not depend on
the latents -- the residual chi-square about the storage origin
and the (masked) Gaussian normalisation.
names: the latents this term is over, in the order they are ravelled.
shapes: each latent's shape, in the same order. ``()`` for a scalar.SqrtInfo.width
width(self) -> intNumber of columns the named latents ravel to.
SqrtInfo.ravel
ravel(self, values: dict[str, Any]) -> jax.ArrayFlatten ``{name: array}`` into this term's column order.SqrtInfo.log_prob
log_prob(self, values: dict[str, Any]) -> jax.ArrayThe log-density this term encodes, at the given latent values.
SqrtInfo.information
information(self) -> jax.Array``F = R^T R`` -- the information this term carries. May legitimately be singular: one epoch usually constrains only a subspace, and only the campaign total plus the prior need be positive definite. Named ``information`` rather than ``fisher`` because it is the observed information of a stored term, not a Fisher matrix of a model -- ``exact/fisher.py`` owns that word here.
SqrtInfo.null
null(cls, names: tuple[str, ...], shapes: tuple[tuple[int, ...], ...]) -> SqrtInfoA term that says nothing -- the identity of :meth:`combine`. Square rather than zero-row so an accumulator's pytree keeps a fixed treedef across a whole campaign, which is what stops ``jit`` retracing once per epoch.
SqrtInfo.combine
combine(cls, first: SqrtInfo, second: SqrtInfo) -> SqrtInfoThe term whose log-density is the sum of the two given ones.
Stack the augmented factors and re-triangularise. Writing
``y = [x; -1]`` so that ``[R | z] y = R x - z``, the stacked product
has the same norm as its triangular factor because ``Q`` has
orthonormal columns::
||R_a x - z_a||^2 + ||R_b x - z_b||^2 = ||R_tot x - z_tot||^2 + rho^2
``rho`` is the corner of the triangular factor -- the part of the two
residuals that no single quadratic form in ``x`` can express. It is a
constant, so it belongs in the offset; dropping it leaves every
combined term wrong by an amount that grows with the campaign and is
invisible in the posterior's SHAPE, which is why the tests here
compare absolute log-densities.marginalise
Defined in bayesmith.marginal.sqrtinfo, line 290.
marginalise(info: SqrtInfo, block: Sequence[str]) -> SqrtInfoIntegrate named latents out of a square-root information form, exactly.
Permute the block's columns first, re-triangularise, and drop the leading
rows and columns. That drop **is** the Schur complement, and the Gaussian
integral over the block contributes exactly
``+ (n_block/2) log(2 pi) - sum log|R_bb,ii| - 0.5 rho^2``
and nothing else. In particular it does **not** contribute the block's own
prior normalisation: whoever appended the prior rows owns
``-sum(log(std)) - (n/2) log(2 pi)``, and the two ``2 pi`` halves cancel
while ``sum(log(std))`` has nothing to cancel against. rheplicant shipped
that second term missing once -- 1.07 nats for three nuisances at
``std=0.7``, 27.47 for twenty-five at ``std=3``, and **exactly zero at**
``std=1``, which is how a probe built on unit priors passed. A constant is
invisible in a posterior's shape, so the tests here compare ABSOLUTE
log-densities and use a non-unit prior.
Marginalising **every** name is legal and returns a zero-width term whose
``log_prob({})`` is the marginal likelihood; marginalising **nothing** is
likewise legal and is the identity on the log-density.
Args:
info: the joint form, prior rows already appended by the caller.
block: which names to integrate out.
Returns:
A term over ``info``'s remaining names, in their original relative
order, whose log-density is the integral of ``info``'s over ``block``.
Raises:
StructureError: if a name is repeated or not in ``info``; if the
re-triangularisation is not finite; or if the block is not
constrained -- an unconstrained direction makes the integral
divergent, and finite arithmetic returns a large plausible number
for it rather than an infinity anyone would notice.marginalise_arrays
Defined in bayesmith.marginal.sqrtinfo, line 230.
marginalise_arrays(factor: jax.Array, target: jax.Array, offset: jax.Array, n_block: int) -> tuple[jax.Array, jax.Array, jax.Array, jax.Array]The Schur complement in square-root form, with no Python control flow.
The block is the **leading** ``n_block`` columns; permuting is the
caller's job, because a caller that already knows its layout should not
pay for a name lookup once per epoch inside a ``lax.scan``.
This exists because :func:`marginalise` cannot be traced: it concretises
twice, for the pivot scale and for the dtype's floor. A filter that
evaluates this arithmetic inside a likelihood -- under ``scan``, and
DIFFERENTIATED with respect to a transition's own parameters -- needs the
half that traces. ``grad`` is the half that matters: a correlation time
inferred rather than pinned is differentiated on every leapfrog step.
**The refusal is not weakened, it is moved -- and this function cannot
make it.** What :func:`marginalise` catches is a block that does not
constrain itself, which makes the integral divergent. Under a trace that
judgement is unavailable: it needs a comparison against a value, and there
is no value. What this hands back instead is the evidence -- ``pivots``,
as data -- so an eager caller can judge them and a traced one can carry
them out.
Args:
factor: ``(r, n)`` array ``R``, block columns first.
target: ``(r,)`` array ``z``.
offset: scalar; the constant this term already carries.
n_block: how many leading columns to integrate out. ``0`` is legal and
is the identity on the density -- the re-triangularisation folds
any excess rows into the corner and the offset absorbs it.
Returns:
``(factor, target, offset, pivots)`` -- the retained form, the offset
with the Gaussian integral's constant folded in, and ``|diag(R)|`` of
the re-triangularisation so a checked caller can test it.compress
Defined in bayesmith.marginal.compress, line 67.
compress(design: Mapping[str, jax.Array], data: jax.Array, precision: Any, shapes: Mapping[str, tuple[int, ...]], *, offset_prediction: jax.Array | None=None) -> SqrtInfoOne epoch's likelihood as a sufficient statistic over the named latents.
Args:
design: ``{latent: (n_data, n_i)}`` blocks. **Column order follows
this mapping's iteration order**, which becomes the term's latent
order -- the same convention rheplicant's compressor uses, and the
one :meth:`~bayesmith.marginal.sqrtinfo.SqrtInfo.combine` needs
two terms to agree on.
data: ``(n_data,)`` this epoch's observations.
precision: the epoch's ``N^-1``, from
:func:`~bayesmith.exact.gaussian.precision_at` or built directly.
shapes: ``{latent: shape}``, so the term can ravel a value dict back.
offset_prediction: the part of the prediction that does not depend on
the latents, subtracted from ``data`` before whitening. ``None``
means zero.
Returns:
A :class:`~bayesmith.marginal.sqrtinfo.SqrtInfo` whose ``log_prob``
is this epoch's exact Gaussian log-likelihood -- **normalisation
included**, because that constant is what a campaign's evidence is
made of and dropping it is invisible in every posterior.
Raises:
StructureError: if a design block's rows do not match the data; if a
block's columns do not match its declared shape; or if the
covariance is correlated AND reports an unobserved sample, which
has no exact meaning -- the observed submatrix of a stationary
covariance is not stationary.compress_epoch
Defined in bayesmith.marginal.compress, line 318.
compress_epoch(global_design: Mapping[str, jax.Array], data: jax.Array, precision: Any, global_shapes: Mapping[str, tuple[int, ...]], *, nuisance_design: Mapping[str, jax.Array] | None=None, nuisance_shapes: Mapping[str, tuple[int, ...]] | None=None, nuisance_prior_std: Mapping[str, Any] | None=None, nuisance_prior_mean: Mapping[str, Any] | None=None, offset_prediction: jax.Array | None=None) -> SqrtInfoOne epoch, with its own nuisances integrated out, over the globals.
The streaming analysis in one call: compress the epoch over BOTH sets,
append the nuisances' prior rows, and marginalise them. What comes back is
a term over the globals alone, and folding those across a campaign is
:meth:`~bayesmith.marginal.sqrtinfo.SqrtInfo.combine`.
**The prior rows go in before the marginalisation, not after**, and that
is not a convenience. A per-epoch nuisance is integrated exactly once, so
a prior that arrives later has nowhere to be applied -- and without one
the block need not constrain itself, which makes the integral divergent.
:func:`~bayesmith.marginal.sqrtinfo.marginalise` refuses that case rather
than returning the large plausible number finite arithmetic would give.
Args:
global_design: ``{latent: (n_data, n_i)}`` for the latents that
SURVIVE the epoch.
data: ``(n_data,)`` this epoch's observations.
precision: the epoch's ``N^-1``.
global_shapes: ``{latent: shape}`` for the survivors.
nuisance_design: ``{latent: (n_data, n_i)}`` for the latents
integrated away inside the epoch. ``None`` means there are none,
and this reduces to :func:`compress`.
nuisance_shapes, nuisance_prior_std, nuisance_prior_mean: the
nuisances' shapes, prior widths and prior centres.
``nuisance_prior_mean`` defaults to zero.
offset_prediction: as for :func:`compress`.
Returns:
A term over ``global_design``'s latents whose ``log_prob`` is the
epoch's marginal log-likelihood -- exactly
``log N(d | A_g x_g + A_n m_n + c, N + A_n S_n A_n^T)``, which is the
independent formula the tests compare against.
Raises:
StructureError: from :func:`compress` or from
:func:`~bayesmith.marginal.sqrtinfo.marginalise`; and if
nuisances are named without a prior, which the marginalisation
has no way to make convergent.epoch_joint
Defined in bayesmith.marginal.compress, line 268.
epoch_joint(global_design: Mapping[str, jax.Array], data: jax.Array, precision: Any, global_shapes: Mapping[str, tuple[int, ...]], *, nuisance_design: Mapping[str, jax.Array] | None=None, nuisance_shapes: Mapping[str, tuple[int, ...]] | None=None, nuisance_prior_std: Mapping[str, Any] | None=None, nuisance_prior_mean: Mapping[str, Any] | None=None, offset_prediction: jax.Array | None=None) -> tuple[SqrtInfo, tuple[str, ...]]The epoch's term over BOTH sets, prior rows appended -- before integrating.
Split out of :func:`compress_epoch` so the assembly exists once and the
two marginalisers can both use it. ``compress_epoch`` passes the result
to :func:`~bayesmith.marginal.sqrtinfo.marginalise`, which checks and
cannot be traced; a campaign passes it to
:func:`~bayesmith.marginal.sqrtinfo.marginalise_arrays`, which traces and
hands back ``pivots`` for the caller to check on every epoch at once.
The same split as ``marginalise``/``marginalise_arrays`` one level down,
and for the same reason.
**The nuisances are the LEADING columns**, which is what lets a traced
caller skip the name permutation entirely.
Returns:
``(joint term, nuisance names)``.nuisance_prior
Defined in bayesmith.marginal.compress, line 203.
nuisance_prior(names: tuple[str, ...], shapes: Mapping[str, tuple[int, ...]], prior_std: Mapping[str, Any], prior_mean: Mapping[str, Any] | None, over: tuple[str, ...], over_shapes: Mapping[str, tuple[int, ...]]) -> SqrtInfo``N(x_n | m_n, S_n)`` as a ``[R | z]`` term over the FULL column set. Rows ``I/std`` in the nuisance columns, zero everywhere else, target ``m_n/std``, and the offset the prior's own normalisation: ``-sum(log std) - (n/2) log 2 pi``. **That offset is the caller's, and this is the caller.** :func:`~bayesmith.marginal.sqrtinfo.marginalise` contributes ``+(n/2) log 2 pi - sum log|R_bb,ii| - 1/2 rho^2`` and deliberately not the prior's normalisation -- the two ``2 pi`` halves cancel and ``sum(log std)`` has nothing to cancel against, which is how rheplicant shipped it missing and why a unit-prior fixture could not see it. Emitted over the full column set rather than over the nuisances alone so it can be :meth:`~bayesmith.marginal.sqrtinfo.SqrtInfo.combine`d with the epoch's own term, which is what puts the prior rows in BEFORE the marginalisation rather than after.
observed_mask
Defined in bayesmith.marginal.compress, line 49.
observed_mask(precision: Any) -> jax.Array | NoneWhich samples this covariance says were taken, or ``None``. ``None`` when the question does not apply -- a covariance with no per-sample sigma has no per-sample "was this taken", and :func:`compress` refuses to guess one. Otherwise a boolean array, ``True`` where ``sigma`` is finite. Reads through :func:`~bayesmith.exact.precision.per_sample_sigma`, which is the one place in the package that asks which implementation a ``Precision`` is, rather than adding a second.
Factorization
Defined in bayesmith.marginal.factorize, line 86.
class Factorization(eqx.Module)
epoch_plate: str = eqx.field(static=True)
survivors: tuple[str, ...] = eqx.field(static=True)
per_epoch: tuple[str, ...] = eqx.field(static=True)A campaign's latents, split by the extent of data they are constant over.
Attributes:
epoch_plate: the plate whose axis indexes epochs.
survivors: latents outside that plate -- constant over the whole
campaign, so they are what a streamed term is over.
per_epoch: latents inside it -- integrated away within their own
epoch, and therefore never accumulated.factorize
Defined in bayesmith.marginal.factorize, line 184.
factorize(graph: Graph, epoch_plate: str, *, at: dict[str, Any] | None=None, rtol: float=LEAK_RTOL) -> FactorizationDerive the campaign's factorization from the graph, and check it.
Plate membership gives the partition; :func:`epoch_leakage` and
:func:`~bayesmith.exact.linearity.check_linearity` decide whether it is
TRUE. Both are refusals rather than diagnostics, for measured reasons:
* a per-epoch latent that reaches its neighbours' data would be integrated
once per epoch instead of once, and the wrong answer is finite and
plausible;
* a map that is not affine in the survivors and per-epoch latents JOINTLY
makes the fold's design -- one jacobian, taken at zeros -- describe the
point rather than the model. Measured on `mu = g * n`: the jacobian
at zeros is identically zero, so the leakage probe read 0.0 and the
folded term was a finite constant independent of the data.
Affinity is checked BEFORE leakage, because the leakage score is one of
those jacobians at one point -- on a non-affine map it is meaningless.
Args:
graph: the model.
epoch_plate: the plate whose axis indexes epochs.
at: values for the latents this check holds fixed. Defaults to zeros,
which is enough because the leakage of an affine model does not
depend on the point -- and a model that is not affine in its
per-epoch latents cannot be folded at all.
rtol: relative leakage above which a latent is refused.
Raises:
StructureError: if the plate is not in the graph; if no latent
survives the epoch, which would leave nothing to accumulate; if
the observed predictions are not jointly affine in the survivors
and per-epoch latents; or if a per-epoch latent is not
epoch-local.epoch_leakage
Defined in bayesmith.marginal.factorize, line 102.
epoch_leakage(graph: Graph, name: str, epoch_plate: str, at: dict[str, Any]) -> floatHow far outside its own epoch a per-epoch latent reaches, relatively.
``0.0`` means the latent touches only its own epoch's observations, which
is what makes integrating it once per epoch the same as integrating it
once. Anything else means it does not.
The measurement is one forward-mode jacobian of the observed predictions
with respect to the latent -- the same map :func:`compress` will build a
design from -- so it asks about the model the fold will actually see,
not about a declaration.
Returns:
``max |off-epoch| / max |any|`` over every observed node, or ``0.0``
if the latent moves nothing at all. Relative to the latent's OWN
largest sensitivity, so a latent the model barely uses is judged on
the same scale as one it leans on -- an absolute floor would wave
through a leak in a weakly-used nuisance and refuse a well-used one
for roundoff.compress_campaign
Defined in bayesmith.marginal.campaign, line 273.
compress_campaign(graph: Graph, epoch_plate: str, *, at: dict[str, Any] | None=None, factorization: Factorization | None=None) -> SqrtInfoFold every epoch into one term over the survivors. Holds one epoch's term at a time, which is the point of the exercise: the raw data and the forward evaluation are gone by the end, and what remains is a sufficient statistic whose size does not grow with the campaign.
epoch_terms
Defined in bayesmith.marginal.campaign, line 47.
epoch_terms(graph: Graph, epoch_plate: str, *, at: dict[str, Any] | None=None, factorization: Factorization | None=None, stacked_only: bool=False) -> list[SqrtInfo]One ``[R | z]`` term per epoch, over the survivors.
Each epoch's own per-epoch latents are integrated out inside it, at the
prior the graph declares for them -- so what comes back is a list of terms
all over the same latents, ready to fold.
Args:
graph: the model. Exactly one observed node, plated on ``epoch_plate``.
epoch_plate: the plate whose axis indexes epochs.
at: values for latents this holds fixed. Defaults to zeros, which is
what an affine model's design does not depend on.
factorization: a partition already derived and checked. ``None``
derives one, **including the leakage check** -- pass one only if
you have already run that.
stacked_only: return the terms as STACKED arrays rather than a list,
so :func:`compress_campaign` can fold them under ``lax.scan``
instead of unpacking E of them only to re-stack. Internal; the
list is the shape a caller wants.
Raises:
StructureError: from :func:`factorize`; or if the graph does not have
exactly one observed node plated on ``epoch_plate``, which is what
makes "one epoch's data" a thing that can be sliced.epoch_observation
Defined in bayesmith.marginal.campaign, line 296.
epoch_observation(graph: Graph, epoch_plate: str) -> strThe one observed node whose data IS this plate's epochs, or a refusal. PUBLIC because a caller deciding *whether* to fold needs the same structural verdict the fold itself makes, and cannot get it from ``factorize`` alone -- see :func:`~bayesmith.dispatch.streaming.streaming_route`, which was the package's only cross-module private import until this name lost its underscore.
coherent_mode
Defined in bayesmith.marginal.diagnostics, line 68.
coherent_mode(terms: Sequence[SqrtInfo], values: Mapping[str, Any]) -> dict[str, float]Is there a fault common to every epoch? -- the detectable half.
A common mode contributes no variance, so it is invisible to split-half,
to leave-one-out and to the posterior width. What it moves is the MEAN of
the per-epoch residual, and a mean over ``N`` epochs is resolved at
``sqrt(N)`` -- which is why this reports a z-score rather than a
magnitude.
``scatter`` is reported beside the mean and is the part to read second: a
mean-level fault leaves it near 1, while an under-estimated noise model
raises both together. The two failing differently is what tells them
apart.
**Not the same statistic as :func:`template_modes`, and the pair is worth
reading together** (ledger D45). This one evaluates ``||R x - z||^2`` at a
point YOU choose and can be re-asked at any value of the survivors; that
one reads what was recorded when the raw data still existed -- the
residual left after each epoch's own best fit, which no later call can
recompute -- and it carries the named template projections this one has
no way to know about. They coincide only at the per-epoch minimiser, which
is not the point a campaign evaluates at. This one also reads its dof off
the FIRST term's row count, so a campaign whose epochs differ in flagging
wants ``template_modes``, which sums each epoch's own.
Args:
terms: the campaign's stored per-epoch terms.
values: the point to evaluate the residual at -- normally the
campaign's own posterior mean.
Returns:
``{"n_epochs", "chi2_mean", "chi2_dof", "chi2_z", "scatter"}``.
Raises:
StructureError: if the campaign is empty, or if the epochs do not
agree about which latents they are over -- a mean over terms that
describe different parameters is not a mean of anything.epoch_chi_square
Defined in bayesmith.marginal.diagnostics, line 45.
epoch_chi_square(terms: Sequence[SqrtInfo], values: Mapping[str, Any]) -> np.ndarray``||R x - z||^2`` per epoch, at one point over the survivors. The residual each stored term still carries -- computed from the term alone, which is what makes it available a thousand epochs after the data was discarded.
refuse_undeclared_coherent_error
Defined in bayesmith.marginal.diagnostics, line 142.
refuse_undeclared_coherent_error(declared: bool, *, what: str='a coherent error shared by every epoch') -> NoneRefuse to report on the half no statistic can see.
**There is no number to return here, and that is the point.** The in-span
half of a coherent error is absorbed into the survivors identically in
every epoch: it leaves no residual, so :func:`coherent_mode` reads noise,
a held-out epoch reads noise, split-half agrees to roundoff -- and the
answer is biased, without limit, as the campaign grows. Measured in this
package's own tests: an in-span injection leaves ``chi2_z`` inside its
noise band while displacing the estimate by exactly the injected amount.
So the honest interface is a refusal keyed on what the analyst DECLARES,
not a report keyed on what the numbers show. A campaign that has not
declared this class is not clean; it is unexamined, and the two must not
read the same.
Args:
declared: whether the analyst has stated that this error class has
been bounded by some means OUTSIDE the campaign's own residuals --
an end-to-end simulation, a redundant instrument, an independent
measurement of the shared component.
what: names the class, for the message.
Raises:
StructureError: unless ``declared``.ResidualSummary
Defined in bayesmith.marginal.compress, line 407.
class ResidualSummary(NamedTuple)
chi2: float
dof: int
reduced_chi2: float
template_names: tuple[str, ...]
projections: AnyA hundred bytes saying what one epoch's residual shows, after the fit.
The design this implements asks for "a fixed-size (~100 byte) per-epoch
residual summary: chi-square, DOF, and the residual projected onto a
handful of NAMED systematic templates", so that "a common mode then shows
as a nonzero MEAN at sqrt(N) significance with correct scatter". That
sentence is quoted rather than cited: it lived in a numbered list item of
an rheplicant design note which is untracked there, so no reader of this
repository could reach it, and the "Section 9.3" this docstring used to
say resolved to nothing -- while THIS repository has since grown an
unrelated `### 9.3`, so grepping the number returns a confident wrong
answer (R9).
Recorded at compression time because that is the last moment the raw data
exists. A campaign keeps these and discards everything else, so a fault
found a thousand epochs later is found in these numbers or not at all.
Attributes:
chi2: the whitened residual left AFTER the epoch's own best fit -- the
part no value of the latents could have absorbed.
dof: observed samples minus the rank of what the epoch fitted.
reduced_chi2: ``chi2 / dof``, or ``nan`` when ``dof`` is zero. An
epoch whose design saturates its data has no residual to speak of,
and a zero here would read as a perfect fit.
template_names: the named systematic shapes, in the stored order.
projections: ``(len(template_names),)`` -- each a standard normal
under the null -- or ``None`` when no templates were given. Stored
positionally, which is why :func:`epoch_residuals` refuses a
campaign whose epochs name them differently.residual_summary
Defined in bayesmith.marginal.compress, line 446.
residual_summary(fitted: SqrtInfo, precision: Any, *, templates: Mapping[str, jax.Array] | None=None) -> ResidualSummaryOne epoch's out-of-span residual and its named template projections.
**``fitted`` must be the term over EVERY column the epoch fits, nuisances
included, and it must not have been marginalised or QR'd yet** -- so it is
what :func:`compress` returns for the union of the global and nuisance
designs, whose ``factor`` is the whitened design and whose ``target`` is
the whitened residual.
Passing the GLOBAL block alone is wrong wherever nuisances are used, and
wrong in a way no shape catches: the nuisance's contribution stays in the
residual, so the chi-square is inflated by whatever the nuisance explained
while the dof is over-counted by its rank. Upstream measured that over
4000 clean epochs of a two-global design with a three-column nuisance --
including the nuisance gives ``dof = 3`` and a mean chi-square of 3.0000,
excluding it gives ``dof = 6`` against 28.44, which is a nine-sigma
per-epoch detection of nothing at all.
The residual is taken after the epoch's own best fit, so what is left is
the OUT-OF-SPAN half of any error. That is deliberate and it is also the
limit: the in-span half is absorbed into the latents identically in every
epoch and leaves nothing here, which is why a systematic floor is a
declaration rather than a measurement.
Templates are projected onto the part of themselves that is also out of
span. One lying inside the design's column space projects to exactly
``0.0`` and says so, rather than reporting a small number that reads like
a null result -- and that reading of zero is exclusive, because a
non-finite or wrongly shaped template is refused by name first.
**"Inside the span" is a RELATIVE test, and it has to be.** What is left
of an in-span template after projection is roundoff, not zero -- measured
at ``6.0e-07`` of its own norm in float32 for a template that is literally
a design column -- and dividing by that norm returns an arbitrary unit
vector's dot with the residual, measured ``-0.2517``. The cut is
``sqrt(eps)`` of the arithmetic in hand. It is the same FORMULA
:func:`~bayesmith.exact.reduced_basis.numerical_rank` uses and a different
RULE: that one is generous because a quadratic form squares the
conditioning, this one because declaring a template in-span makes it
quieter, and quieter is the safe direction for a detection statistic.
The projector is the flat-prior one: everything the columns COULD explain
is removed. That is the conservative direction -- it can only make a
template quieter -- and it is what makes the null exact rather than
dependent on whatever prior the marginalisation used.
Args:
fitted: the pre-marginalisation term over every fitted column.
precision: the epoch's ``N^-1``, so a template given in the model's own
units is whitened the same way the design was. The SAME object
:func:`compress` was given.
templates: ``{name: (n_data,)}`` named systematic shapes, in model
units. ``None`` for an epoch that names none.
Returns:
A :class:`ResidualSummary`.
Raises:
StructureError: if a template is not finite where the epoch observed,
or does not have one entry per sample.epoch_residuals
Defined in bayesmith.marginal.diagnostics, line 192.
epoch_residuals(summaries: Sequence[ResidualSummary]) -> tuple[dict[str, Any], ...]The per-epoch table, in the order the summaries were given.
Per-epoch half of the design quoted in
:class:`~bayesmith.marginal.compress.ResidualSummary` (R9).
Order is preserved rather than sorted by anything. For a chained campaign
the archive order IS the time order, and a diagnostic that reordered it
would make a drift look like scatter.
**Distinct from :func:`coherent_mode`, and the difference is the input.**
That one evaluates ``||R x - z||^2`` at a point you choose, from the
stored term, and can be re-asked at any value of the survivors. This one
reads what :func:`~bayesmith.marginal.compress.residual_summary` recorded
when the raw data still existed: the residual left after the epoch's OWN
best fit, which no later call can recompute. The two coincide only at the
per-epoch minimiser, which is not the point a campaign evaluates at.
Args:
summaries: one per epoch, from
:func:`~bayesmith.marginal.compress.residual_summary`.
Returns:
One dict per epoch with ``"chi2"``, ``"dof"``, ``"reduced_chi2"`` and
``"templates"``, the last being ``{name: projection}``.
Raises:
StructureError: for an empty campaign, or one whose epochs name
different templates -- see :func:`refuse_mixed_templates`.refuse_mixed_templates
Defined in bayesmith.marginal.diagnostics, line 240.
refuse_mixed_templates(summaries: Sequence[ResidualSummary]) -> tuple[str, ...]The one template list this campaign shares, or a refusal.
The projections are stored as a bare array whose meaning is POSITIONAL,
and the summary is read by name across a whole campaign. Two epochs naming
``("gain_ripple", "ground_pickup")`` and ``("ground_pickup",
"gain_ripple")`` have compatible shapes and incompatible contents, so
averaging them positionally reports each template's mean as a mixture of
both -- a finite, plausible number about nothing.template_modes
Defined in bayesmith.marginal.diagnostics, line 271.
template_modes(summaries: Sequence[ResidualSummary]) -> dict[str, Any]Is there a fault shaped like one of the named templates?
Named-template half of the design quoted in
:class:`~bayesmith.marginal.compress.ResidualSummary` (R9), counterpart to
:func:`coherent_mode`'s unnamed one. Neither substitutes for the other: a
chi-square z needs no guess about what the fault looks like and therefore
cannot say what it IS, while a named template says what it is and is
silent when the guess was wrong.
Under the null each projection is a standard normal in its own right --
that is what the whitening buys -- so the campaign mean is a z-score at
``sqrt(N)`` with no calibration step.
``scatter`` is reported beside every mean and is the part to read second.
A mean-level fault leaves it at 1.0 while an under-estimated noise model
raises both together, and the two failing differently is what tells them
apart. It is the population spread about the sample mean (``ddof=0``), so
a shifted mean does not inflate it.
**The chi-square z here uses each epoch's OWN dof**, summed, rather than
the first epoch's. A campaign whose nights differ in flagging has a
different null per night, and it is the SUM that is chi-square distributed:
``Var(sum chi2_k) = 2 sum k``. That reduces to the textbook form for a
uniform campaign and stays exact for a ragged one, which
:func:`coherent_mode` -- reading a row count off the first term -- cannot
do.
Returns:
``{"n_epochs", "chi2_mean", "chi2_dof", "chi2_z", "templates"}``,
where ``chi2_dof`` is the shared dof or ``None`` for a ragged
campaign, and ``templates`` is ``{name: {"mean", "scatter", "z"}}``.
Raises:
StructureError: for an empty campaign or mixed template lists.held_out_z
Defined in bayesmith.marginal.diagnostics, line 372.
held_out_z(terms: Sequence[SqrtInfo], prior_fisher: Any, prior_mean: Any=None) -> tuple[dict[str, float], ...]How surprising is each epoch to the rest of the campaign?
For a linear-Gaussian model this is exact and needs no simulation. Write
the leave-one-out posterior ``N(mu_-e, Sigma_-e)`` and the epoch's own
factor ``[R_e | z_e]``. Then ``z_e = R_e theta + eps`` with unit-covariance
``eps`` -- that is what the square-root form MEANS -- and ``mu_-e - theta``
is independent of ``eps``, so::
m = R_e mu_-e - z_e ~ N(0, I + R_e Sigma_-e R_e^T)
and ``m^T (I + V)^-1 m`` is chi-square on ``rank(R_e)`` degrees of
freedom. The returned ``z`` is that, standardised.
**Computed by subtracting one summand, never by downdating the
accumulator.** A QR accumulation cannot be un-summed stably. The campaign
total is formed once in ``(F, b)`` form and one epoch's contribution
subtracted per row, which is ``O(N)`` rather than ``O(N^2)`` and loses at
most ``log10(N)`` of float64's sixteen digits -- affordable because this
is offline, and safe because it removes one positive-semidefinite summand
out of N rather than a triangular factor out of its own product.
**What it can see, and what it cannot.** A single rogue epoch stands out
by a wide margin, and so does a campaign whose epochs genuinely differ.
But where every epoch carries the SAME design -- the realistic case -- a
coherent error's in-span half shifts ``z_e`` and ``mu_-e`` by amounts that
cancel in ``m``, and the clean and biased campaigns return the same
scores. ``tests/marginal/test_held_out.py`` measures both directions
rather than asserting them. Read this beside
:func:`refuse_undeclared_coherent_error`, never instead of it.
Args:
terms: the campaign's stored terms, in any order -- this statistic is
exchangeable even when the campaign is not, because each epoch is
scored against all the others.
prior_fisher: ``F_prior`` over the same latents in the same column
order. **Required, not optional**: ``Sigma_-e`` is singular at
small N without one.
prior_mean: the prior's mean over the raveled columns. Zero if absent.
Returns:
One dict per epoch, in the order given, with ``"chi2"``, ``"dof"`` and
``"z"``.
Raises:
StructureError: for an empty campaign, epochs over different latents, a
prior over a different number of columns, a singular leave-one-out
information, or a score that comes out non-finite.shrinkage_power
Defined in bayesmith.marginal.diagnostics, line 551.
shrinkage_power(sigmas: Mapping[int, Any]) -> floatThe fitted exponent of ``sigma_N ~ N^p``. **A sanity check, not a test.**
Returned as a bare float only from here; :func:`shrinkage_report` is what
a caller should print, because it carries the caveat in the same object as
the number.
For a Gaussian model ``sigma_N = (sum_e F_e + F_prior)^-1/2`` does not read
the data. So this quantity is **data-independent**, ``p = -0.5`` holds by
construction, and a uniform rescaling ``F_e -> (1+c) F_e`` cannot move it
at all -- the intercept is free. A diagnostic that could be fooled by
injecting a shared systematic and watching this number is therefore no
diagnostic, which is the whole reason
:func:`refuse_undeclared_coherent_error` exists.
Args:
sigmas: ``{campaign size: posterior widths}``. A scalar width is
accepted and treated as a one-element array. Every size must
report the same number of widths, in the same order.
Returns:
The ordinary-least-squares slope of ``log sigma`` on ``log N``, pooled
over the widths.
Raises:
StructureError: for fewer than two campaign sizes, a non-positive
size, a ragged table, or a sigma that is not finite and strictly
positive.shrinkage_report
Defined in bayesmith.marginal.diagnostics, line 585.
shrinkage_report(sigmas: Mapping[int, Any]) -> dict[str, Any]:func:`shrinkage_power` with its limits attached to it. ``detects_coherent_bias`` is ``False``, always, and it is a FIELD rather than a sentence in a docstring because a number and its caveat travel together or they do not travel. A deterministic error shared across epochs contributes no variance at all, so the posterior width is the same array element for element on a clean campaign and on a badly biased one, and this power is the same to every digit either way.
systematic_floor
Defined in bayesmith.marginal.diagnostics, line 647.
systematic_floor(term: SqrtInfo, prior_fisher: Any, floors: Mapping[str, Any], n_epochs: int) -> dict[str, dict[str, Any]]Has this campaign out-run its own calibration?
The floor is the declared width of a shared calibration product -- one
solution, one beam model, one flag table serving every epoch --
**projected into the latents' units by the analyst**. It is a declaration
and not a measurement, and that is forced rather than lazy: the in-span
half of a coherent error biases the survivors identically in every epoch
and leaves no residual anywhere, so it passes per-epoch chi-square,
split-half and leave-one-out, and NO statistic computed from the stored
terms can recover it.
What the campaign does know is its own width, and that width falls as
``N^-1/2`` while the shared product's does not fall at all. So the whole
content of this function is: when does one pass under the other.
``crossing_epoch`` is **computed from the observed width**, not quoted:
``sigma_N`` is extrapolated as ``sigma_N sqrt(N / N')`` and solved for
``sigma_N' = floor``, giving ``N' = ceil(N (sigma_N / floor)^2)``. That
extrapolation ignores the prior's share, which shrinks as N grows, so it
predicts the crossing marginally early -- the conservative direction.
Args:
term: the campaign's folded term over the survivors.
prior_fisher: ``F_prior`` over the same columns, added to the term's
own information before inverting.
floors: ``{latent: declared width in that latent's units}``. Every
entry must name a latent the term is over and be finite and
strictly positive.
n_epochs: how many epochs the term was folded from. Used only for the
extrapolation.
Returns:
One entry per named latent, each a dict of ``"sigma"``, ``"floor"``,
``"below_floor"``, ``"crossing_epoch"`` and ``"direction"``.
Raises:
StructureError: for a campaign of no epochs, a floor naming a latent
the term is not over, a floor that is not a finite positive width,
a prior of the wrong shape, or information that is not positive
definite.tightest_direction
Defined in bayesmith.marginal.diagnostics, line 609.
tightest_direction(block: Any) -> tuple[float, Any]``(width, unit direction)`` of the narrowest direction of a covariance block. **The tightest DIRECTION, not the tightest coordinate**, and the difference is a basis rotation wide. For any correlated posterior the smallest eigen-direction is below every diagonal entry, so a refusal keyed on ``min(diag(block))`` watches a quantity that is not the first to go under. ``tests/marginal/test_systematic_floor.py`` measures the gap on a near-collinear campaign, where the coordinate widths sit an order of magnitude ABOVE a floor that the tightest direction is well below. The direction is returned because a bare number is half the value: "your error bar is too tight" is not actionable, "your error bar on ``0.707 x0 + 0.707 x1`` is too tight" is. Its sign is fixed by making the largest-magnitude component positive, so the same posterior reports the same vector rather than one that flips with LAPACK's mood. ``(nan, None)`` for a block that is not finite: NaN cannot be eigendecomposed, and a poisoned campaign must report ``nan`` and be refused by the NaN-safe comparison in :func:`systematic_floor` rather than raise a linear-algebra error here.