bayesmith.marginal.compress
Source-derived signatures and contracts. See the reference index for coverage and conventions.
Module contract
One epoch of a linear-Gaussian model, compressed into a ``[R | z]`` term.
The bridge between B9's noise interface and B11's kernel. For
d ~ N(A x + c, N)
the log-likelihood is exactly a square-root information term::
R = N^-1/2 A z = N^-1/2 (d - c)
log p(d | x) = -1/2 ||R x - z||^2 - 1/2 log det (2 pi N)
so :func:`compress` is a whitening and a normalisation, and everything after
it is :mod:`bayesmith.marginal.sqrtinfo`'s arithmetic.
**This module owns the "unobserved sample" concept, and that is B11's first
design decision made rather than assumed.** ``sigma = inf`` means a sample was
not taken. :meth:`~bayesmith.exact.precision.Precision.whiten` already gives
it weight zero with no special case, so the QUADRATIC half needs nothing --
but ``log_normalizer()`` is ``+inf`` there, and it is right to be: a sample
with infinite variance has no density, which is the honest answer to the
question a ``Precision`` is asked. Reading it as ``0`` is a statement that the
sample is UNOBSERVED, which is a modelling concept this layer has and the
interface does not. So the mask lives here, and ``precision.py`` keeps a
normaliser that is never silently wrong.
**Masking is a DIAGONAL concept, and that is measured rather than assumed.**
For a stationary covariance the observed submatrix is neither circulant nor a
subset of the spectrum -- on a 6-point kernel with one sample dropped, its
log-determinant is ``-0.7084`` while the closest subset sum of log-eigenvalues
is ``0.47`` nats away. So an unobserved sample inside a correlated epoch is
refused rather than approximated; see :func:`compress`.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.
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.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.
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)``.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.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.