fit_mppca_samples() fits a mixture of probabilistic principal component
analysers (Tipping and Bishop, 1999): each component covariance is a rank-q
term plus equal noise in every direction. It can be fitted when a component
holds fewer points than there are variables, where a full-covariance fit is
singular, and returns an ordinary fitted mixture.from_cf() fits a Gaussian mixture to a distribution given by its
characteristic function, for targets with no usable density such as stable
laws. The fit reports its error on frequencies
it did not use and warns when the target has heavy tails.gmm_cf() evaluates the characteristic function of a mixture.gmm_cf_distance() returns a lower and an upper bound for the Kolmogorov
distance between a characteristic-function target and a mixture along a
direction, computed by numerical inversion and without sampling. The upper
bound needs a target whose projected characteristic function is integrable.from_cf_elliptical() fits normal components with a common centre and
proportional covariances to an elliptical characteristic function, such as
a multivariate Cauchy or an elliptical stable law. It needs no starting
value and matches the distribution out to about its 0.999 quantile, even
with heavy tails.from_cf_sources() fits a linear mix of independent sources, with or
without Gaussian noise. It returns a list with functions for the
characteristic function, sampling and one-dimensional projections, not a
gmm. The separated coordinates are tested for pairwise independence with
the criterion of Eriksson and Koivunen (2003), reported as
diagnostics$dependence.gmm_cf_variance() and gmm_cf_tail_prob() now return the variance and the
tail probability of a per-unit counterfactual law instead of an error. They
rest on the same assumption as gmm_cf_mean(): the unit keeps its structural
residual when the treatment changes. A printed counterfactual law shows its
variance.gmm_counterfactual(), gmm_cf_mean() and
proxy_retrospective_uplift() no longer calls the per-unit counterfactual
mean identified. It assumes the unit keeps its structural residual when the
treatment changes; if the residual were redrawn, the same mixture would give
a different mean. A printed counterfactual law now labels its mean
"residual kept".gmm_independence_graph() tests each partial correlation with Fisher's z
at level alpha (default 0.05) when the mixture was fitted to data, taking
the sample size from a "sample" or "moment" fit or from the new n
argument. Giving threshold keeps the old rule |r| > threshold, which is
also used when no sample size is known.from_kde() uses is_size = 10000 + 2500 * N importance draws instead of
5000, which reduces overfitting at larger N.from_kde() holds out validation draws by default, so its fits report
validation_kld. Use validation_size = 0L for the old behaviour.fit_kld_em() and fit_proxymix(regime = "kld") draw validation points in
batches of is_size until the standard error of validation_kld is at most
a tenth of the estimate, up to 10 * is_size draws. The old default of
ceiling(is_size / 4) draws left a standard error about as large as the
estimate on a close fit. Give validation_size to fix the number of draws.proxymix_low_quality) no longer fires on
relative ESS below 0.05. With a fixed proposal, relative ESS describes the
proposal and not the fit, so the advisory fired on about half of good
regime "kld" fits and more often as is_size grew. It now fires when a fit
did not converge, is degenerate, or has a held-out KL divergence above 0.3.
gmm_fit_quality() reports that estimate as the new field heldout_kld.
When the target's normalising constant is unknown, it is estimated from the
same validation draws. A fit made with validation_size = 0 is judged on
the new field kld_approx, half the weighted variance of log f - log g.
This approximates the KL divergence only when the proxy is close to the
target, and a proxy wider than the target can keep it below 0.3.proxy_mnar_sensitivity() reports, for each slope, the observed-data
log-likelihood of the selection model (loglik) and whether the fit
converged (converged), and warns when a fit did not converge. Its new
max_iter argument defaults to 500, because fits near missing at random
could stop before converging at the gmm_impute() default of 100.
from_objective() gains scale = "loglik" for negative log-likelihoods:
the default temperature ladder then ends at 1, so the returned map is the
likelihood surface on the box.
proxy_functional_ci() gains tail_warn (default 1e-3) and warns, with
class proxymix_tail_functional, when the functional is a probability
below that level. Its interval covers refit variability only, not the
error of the mixture's tails.
glance() on a fitted proxy reports validation_kld, the held-out KL
estimate.
fit_em_samples() and fit_kld_em() scale the ridge ridge_eps by each
coordinate's variance. They scaled one ridge by the mean variance, which
inflated the fitted spread of coordinates whose variance is far below the
mean: with standard deviations from 0.026 to 5078, the smallest came out 85
times too wide. Fits whose coordinates have similar variances barely
change.
gmm_independence_graph() scales its ridge by each coordinate's variance
for the same reason. A partial correlation of 0.6 between two coordinates
with standard deviation 0.001, beside one with standard deviation 10000,
was reported as 0.
gmm_evidence() keeps proxy draws at which the target density is zero as
zero terms of the average. It dropped them before, which biased log Z
upwards for a target with bounded support and understated its standard
error. The returned n is now the number of draws requested.
A target's log_normalizer is log Z, the log of the integral of
exp(log_density). kld_final_absolute and validation_kld_absolute are
now the raw estimates minus log_normalizer. They added it before, so
passing the correct log Z doubled the offset. The help now states that the raw
estimate against an unnormalised log-density is the KL divergence plus
log Z. Code that stored minus log Z to get a small absolute KL must
change sign. This includes log_normalizer in the internal
from_fb_posterior() specifications.
kld_trace() stays finite when some fitting draws fall where the target
density is zero. Those draws are left out of each trace entry, as they
already were from kld_final.
gmm() rejects a covariance with a non-finite entry or one whose Cholesky
factorisation fails, with an error naming covariances. Such a mixture
was constructed before and failed later in rgmm() or dgmm().
fit_em_samples() and fit_moment_match() refuse a sample target with no
rows.
pgmm(lower.tail = FALSE) sums the upper-tail probabilities of the
components. It subtracted the lower-tail probability from one before,
which lost accuracy far in the tail (a relative error of 7% at eight
standard deviations).
from_objective() no longer evaluates the objective outside
[lower, upper]. Points outside the box get the same penalty as
non-finite values, so a map can no longer settle outside the box.
from_objective() raises at most one low effective-sample-size warning
per call, and only when the final cooling step is low. The effective
sample size at each step is stored in fit@metadata$from_objective$ess.
hellinger_mc() draws n_mc fresh samples from the fit, using seed, for
regime "kld" fits too. The estimate is no longer biased low and its
standard error is no longer NaN.
proxy_confounding_gap() and proxy_regime_segments() warn, with class
proxymix_constant_treatment, when each regime holds a single treatment
arm. The within-regime effect is then zero, so the gap equals the whole
estimated effect and the segment effects are all zero.
The decision family (proxy_cate(), proxy_uplift(), proxy_overlap(),
proxy_decide(), proxy_confounding_gap(), proxy_policy_value(),
proxy_retrospective_uplift(), proxy_regime_segments(),
proxy_identification_report()) previously defaulted the treatment
contrast to t1 = 1, t0 = 0 regardless of the treatment coding
fit_uplift() observed, silently returning an effect on the wrong scale
for any treatment not coded {0, 1} (e.g. a {0, 100}-coded dose). The
default now reads model@treatment_levels, the arms actually observed at
fit time; a t1/t0 matching neither observed level now aborts with a
classed proxymix_treatment_scale_error instead of returning a silently
mis-scaled effect.
John van der Hoek is added as an author.
kld_final is described as an in-sample estimate that reads low, with
validation_kld as the figure to quote.
The from_kde() help states that a Gaussian-kernel KDE is already a
mixture that can be marginalised, conditioned and sampled exactly, and
that compression reduces its cost from n components to N.
The gmm_eos_test() help states that the Andrews subsampling calibration
is asymptotic, assumes stationary ergodic innovations, and can over-reject
in short series. It also states how method = "andrews" differs from
Andrews' (2003) P-test: its p-value counts the tested block, and the block
statistics are computed at the supplied model without re-estimating it.
The mnar() and proxy_mnar_sensitivity() help no longer says the data
cannot identify the missing-not-at-random slope. Under the mixture
selection model the observed data inform the slope through the assumed
shape of the outcome distribution, which the new loglik column shows.
gmm_entropy() and gmm_conditional_entropy() are titled by the quantity
they return by default, the order-2 Renyi entropy.
The package website carries an extended version of each of the twelve
vignettes. Each adds a numerical comparison with other R packages on the
same task. The comparisons are read from stored simulation results, and
the code that produced them is in data-raw/articles/ in the source
repository.
All twelve vignettes were brought to a common quality bar: each now
opens with the question it answers, carries a Why / What / Do / Read /
Limits / What to read next / Reproduce shape, states every number in
prose from a chunk or inline code rather than by hand, draws its figures
with ggplot2 on a colourblind-safe palette with captions and alt text,
tabulates through knitr::kable(), refers to its companions by title,
and closes with its seed and sessionInfo(). One mixture, many
methods gains a section on the decision family (fit_uplift(),
proxy_cate(), proxy_decide(), proxy_identification_report()), which
no vignette previously covered and which shows the
proxymix_treatment_scale_error refusal of an unobserved treatment
contrast; How well a mixture proxies four awkward shapes shows the
proxymix_degenerate_fit refusal of a collapsed importance sample.
Claims that had never been computed are now computed, and three that did
not survive the computation were corrected: the regime (i) and regime
(iii) fits at N = 1 do not agree to within Monte Carlo error, the
donut's residual divergence needed a grid-quadrature reference to be stated
at all, and the fraction of x2 deleted in the imputation example was
near a half rather than the two fifths claimed.
Every figure chunk across the vignette set now carries a fig.cap stating
what the figure shows, and the four uncaptioned or under-interpreted plots
(calibration, from_kde, missing_data, operator_calculus) gained a
trailing prose sentence reading the result off the plot.
entropy.Rmd's order-2 Renyi entropy formula is now typeset as LaTeX
(\(H_2(\mathcal{N}(\mu, \Sigma)) = \dots\)) rather than as inline code.
The in-text short-form citation "Hoek and Elliott (2024)" is now
"van der Hoek and Elliott (2024)" in README.Rmd and the three vignettes
that used the short form (quickstart, three_regimes, calibration);
the surname was silently dropped when the co-author's Authors@R entry
was removed at 0.15.1.
README.Rmd and the two vignettes still calling the deprecated is_mvt()
alias (quickstart, three_regimes) now use the preferred proposal_mvt()
name introduced at 0.13.0.
README.Rmd's stray Unicode em/en dashes are now the LaTeX -- convention
used throughout the rest of the package's prose.
gmm_complete(), proposal_uniform()/proposal_mvn()/proposal_mvt(),
and the gmm_imputation class gained @examples.
Function names in nine vignettes printed with literal square brackets
([gmm_impute()]); they now print as plain code.
gmm_product(), gmm_convolve(), and gmm_mix()
previously carried only the first operand's certificate (and gmm_mix()
carried none); they now combine every operand's certificate conservatively --
not-converged if any operand is, degenerate if any operand is, and the
worst-case of the numeric fields -- with the individual certificates
retained under quality_sources and regime = "composite".fit_proxymix(regime = "kld") now flags when an auto-selected proposal
replaces an unbounded target-support end with a finite, data-derived
working box. Because such a box does not sample the target tail beyond it,
the fit is a truncated-target approximation rather than exact importance
sampling; this is reported with a classed warning
(proxymix_truncated_support_proposal), a support_truncated flag in
ess_summary(), and the finite box in fit@diagnostics$working_bounds.
Pass an explicit dominating proposal for exact importance sampling over
the full support.CITATION, the DESCRIPTION references, and the
documentation.gmm_fit_ensemble() quantifies the
sampling variability of the fitted mixture itself by a Bayesian
(weighted) bootstrap: each replicate re-weights the fit's own
observations with Dirichlet weights and refits by a warm-started
weighted EM. In regime "kld" the observations are the fit's cached
importance draws, so a replicate costs zero new target evaluations.
proxy_functional_ci() turns the ensemble into a percentile interval
for any functional of the proxy (moments, tail probabilities via the
operator calculus, entropies, conditional means) -- functional-space
intervals sidestep component label-switching entirely. Coverage
validated against analytic truths on a correlated Gaussian target
(validation/ensemble_coverage_study.R, R = 200 replications):
nominal 90% intervals cover the two mean coordinates and a tail
probability at 0.885 / 0.855 / 0.855.gmm_eos_test() (estimated-parameter
over-rejection of the chi-square calibration; the subsampling p-value
floor) are now stated in the manual and vignette, and the README
positions the package against its nearest CRAN neighbour (AdMit)
and the MCMC-then-fit route explicitly.select_N(): component-count selection across regimes. With
samples, the smallest BIC (agreeing with mclust::mclustBIC on shared
data); with an evaluable-only target, each candidate is scored on its
held-out validation draw (whose Monte Carlo standard error is now
reported as validation_mc_se) and chosen by the one-standard-error
rule. On an evaluable three-component target the selector recovered
N = 3 in ten out of ten seeds. The scored table is returned, so the
recommendation can be ignored by callers who prefer to choose by eye.Adaptive proposals for regime (iii): fit_kld_em(adapt = "pmc").
The importance proposal is refreshed from the current EM iterate as a
defensive mixture (the fitted components with inflated covariances plus
a defensive_gamma share of the original proposal as a heavy-tailed
anchor) -- the mixture population-Monte-Carlo scheme of Cappé et al.
(2008) with the Owen--Zhou (2000) defensive safeguard. While a batch is
degenerate the refresh fires every iteration with an escalating
inflation floored at a growing fraction of the batch's sample
covariance, so a collapsed iterate walks back out to the target; and
convergence is only accepted on an adapted batch. The per-batch ESS
trace ships in diagnostics$ess_history and the total evaluation cost
in diagnostics$n_target_evals.
Validated against analytic conjugate Gaussians and a ten-dimensional
warped target (validation/adaptive_pmc_study.R in the development
repository): with a centred but badly over-wide proposal, adaptation
restores the relative effective sample size from 0.006 to 0.26 at
p = 20 (5--40x across dimensions) with exact moment recovery; with a
proposal mislocated by three standard deviations per coordinate it
restores 0.0003 to 0.36 on the warped p = 10 target and beats the
fixed proposal by 6--65x through p = 15. The documented boundary: a
badly mislocated light anchor at p = 20 leaves too few effective
draws to carry directional signal at these budgets -- supply a located
proposal there. On the low-dimensional shipped scenarios adaptation is
uniformly at least as good as the fixed proposal (held-out validation
KLD), and adapt = "none" (the default) is byte-identical to the
historical behaviour.
gmm_product() (the normalised pointwise product of two mixtures -- the
conjugate Bayes update, with the marginal evidence returned as
metadata$log_integral and an optional reduce to cap the K1 * K2
component growth), gmm_convolve() (the exact distribution of the sum of
independent mixture variables), gmm_mix() (model averaging /
mixture-of-mixtures), and gmm_mean() / gmm_cov() (the exact first two
moments). pgmm() / qgmm() complete the d/p/q/r quartet in one
dimension. All follow the operator metadata policy: the quality
certificate travels and the provenance chain grows.gmm_evidence(): the normalising constant as a first-class output.
With the fitted proxy as the importance proposal, log Z of the target
(the log marginal likelihood, for a posterior handed over as
likelihood x prior) is estimated in the log domain with a delta-method
standard error and a heavy-tail diagnostic (classed warning
proxymix_heavy_tail when the proxy's tails look lighter than the
target's). Validated against known constants for unimodal and multimodal
targets.gmm_weights(), gmm_means(),
gmm_covariances() read the component parameters without reaching into
the @ property layout; broom-style tidy() (component table) and
glance() (one-row fit summary) register against the generics package
when installed.proposal_uniform(), proposal_mvn(), proposal_mvt() are the
preferred names of the importance-proposal constructors (the historical
is_* prefix reads as a logical predicate); the is_* names remain as
aliases.from_aggregate_likelihood(),
fit_kld_em_collider(), to_apsim_scenarios(), from_simulator()) and
the posterior-producer seam helpers (from_fb_posterior(),
fb_log_posterior_spec(), fb_producer_available(),
mock_fb_posterior()) are internal as of this release: an exported
function whose only behaviour is to error, or whose contract awaits an
unreleased counterpart, does not belong on the public surface. Their
signatures and behaviour are unchanged and remain under test; the
general-purpose S3 generic gmm_target_from_posterior() remains the
public route for external posteriors..as_gate() rejects
unknown gate types with a clear error instead of letting a third-party
gate object fail deep inside the engine.gmm_impute() on complete data now does what its warning says: the
completions are the data verbatim, with no bootstrap refits, no
imputation draws, and no random-number consumption (previously it fitted
and drew anyway).proxy_pool(method = "analytic") on a gated (MNAR / censored) mechanism
announces the downgrade to Rubin's rules instead of switching silently.gmm validator takes a Cholesky-first fast path and pays for an
eigendecomposition only when the Cholesky fails -- mixtures are
constructed inside hot loops (every filter step builds several).gmm_fit_quality(). Every closed-form operator carries the certificate
through unchanged, together with a provenance vector recording the chain
of operations, so the certificate can be read off a marginal, a filtered
belief, or any other derived mixture. Downstream verbs (gmm_entropy(),
gmm_mutual_information(), gmm_independence_graph(), gmm_intervene(),
gmm_counterfactual(), gmm_filter()) read the same certificate and raise
a one-shot advisory (class proxymix_low_quality) when the source fit is
flagged.ESS < min_ess) now flags the fit as degenerate: the warning is classed
(proxymix_low_ess), converged is forced to FALSE, and the new
on_low_ess = "abort" refuses to return the fit at all (classed error
proxymix_degenerate_fit). Previously a fit sitting on one effective draw
could report converged = TRUE with a warning as the only trace.validation_size now defaults to
ceiling(is_size / 4) rather than 0, so the overfit-vs-generalise
diagnostic exists on every regime-(iii) fit. Pass validation_size = 0L
to disable.autoplot() is now registered for plain gmm objects too, so
operator-calculus results (marginals, conditionals, filtered beliefs) plot
directly rather than only freshly fitted proxies.fit_kld_em() notes
p > 5 and warns (classed proxymix_high_dimension) beyond p = 10,
instead of only the wrapper entry points knowing the scaling story.proxymix_low_ess,
proxymix_degenerate_fit, proxymix_high_dimension, proxymix_support,
proxymix_nonmonotone, proxymix_low_quality -- pipelines can
condition-handle instead of matching message text.fit_kld_em(seed =)
previously seeded only the importance-sampling draw; the initialisation
resample, the kmeans pass, and empty-component reseeds consumed the ambient
random-number stream, so two calls with the same seed could return
different fits. All internal draws now derive from the seed. Fits under a
given seed therefore differ from 0.11.x. fit_em_samples(seed =) likewise
drives its multi-start (previously hard-coded restart seeds).from_objective() derives a distinct seed per cooling step: previously the
importance-sampling stream was re-seeded identically at every temperature,
so with the default exploration = 0.5 half of the draws were
byte-identical across the whole ladder and a basin missed by the first
exploration draw was never probed again.1e-8 and 1e+8. The floor is constant within a fit -- a ridge relative
to the component's own diagonal would shrink together with a collapsing
component and stop flooring exactly when needed. At unit data scale the
behaviour is unchanged to first order; the operators' tiny hygiene ridge
is unchanged.proxymix_nonmonotone) if their objective decreases beyond numerical
tolerance, which the documentation previously claimed and the code did
not do.p = 115 at scale 1e5), which made every
merge cost 0/Inf/NaN and the merge order arbitrary.hellinger_mc() reports the standard error of the self-normalised
weighted estimator; the previous sd(ratio)/sqrt(n) treated skewed
importance weights as uniform and could understate the Monte Carlo error
by orders of magnitude.O(K^2) pairs were
re-costed, each with its own Cholesky, after every single merge -- the
bottleneck inside long Gaussian-sum-filter runs).pmax reduce
rather than a per-row loop, on every E-step's hot path.gmm_kld() at K = 1 (the package's namesake divergence previously had
only positivity checks), closed-form Gaussian Hellinger for
hellinger_mc(), quadrature for the logit-link gated moments (the shipped
default link previously had no moment-level check), and a probit-link
end-to-end MNAR recovery.H(AX + b) = H(X) + log|det A|), marginal-vs-joint quadrature
consistency, gmm_observe with a near-noiseless selection row against
gmm_conditionalise, post-operator mass preservation, and an
rgmm/dgmm Kolmogorov-Smirnov cross-check.fit_kld_em() previously judged convergence on the relative change
of the importance-sampled KLD trace, whose magnitude carries the arbitrary
-log Z(f) offset of an unnormalised target -- so the same target, shifted
by a constant, could stop after two iterations where the unshifted run used
a hundred. Convergence is now judged on the importance-weighted EM objective
Q(theta) = sum_n W_n log g(x_n), which never touches the constant.
Iteration counts (and therefore fits that previously stopped early or late
for this reason) can change; a regression test pins the invariance.gmm_eos_test() now resolves its model specification through the same
machinery as gmm_filter(): dynamics and measurement offsets (b, d) are
honoured (they were previously ignored without notice, silently changing the
innovations), function-valued (time-varying) specifications are accepted,
and Gaussian-sum (mixture) noise is rejected with an informative error since
both calibrations assume Gaussian innovations.NaN entries from the sum
(previously they were neutralised only in the row maximum and could still
poison the result) and returns +Inf for a row containing +Inf (previously
-Inf, the exact inversion).NaN through
the M-step and abort the EM with an opaque Cholesky error.fb_producer_available()
now reads the producer package name from
getOption("proxymix.producer_package") instead of hard-coding one.README gains a GitHub installation section, a "Why not MCMC?" positioning
note with the dimensional limits of regime (iii) stated up front, a
"Which function do I need?" routing table, and the full eleven-vignette
list; the quick-start example now runs to convergence.mnar() /
censored() mechanisms available since 0.11.0 and cross-references the
companion vignette.inst/validation/contract-conformance.R is now self-contained: a minimal
serial sweep driver ships in the file, so the two-sided case set is
executable anywhere (an external driver can still be injected via
PROXYMIX_CONFORMANCE_ENGINE).withr is now a hard dependency (Imports, previously Suggests). It is
called on the core seeded-sampling paths (annealing, entropy estimation,
importance-sampled KLD-EM fitting, mixture reduction, the gmm_target
sample constructors, initialisation, and the seeded diagnostics), so a clean
installation without withr would have failed at runtime. Declaring it as the
runtime dependency it already is removes that failure mode._pkgdown.yml site url: now resolves to the
package's GitHub Pages host (max578.github.io/proxymix), and CITATION.cff
version and date-released track the released version.inst/validation/ (the
authored validation cases an external conformance harness executes to measure
coverage of the documented function surface). No change to package code; this
gives the previously local-only validation material a tracked, travelling home.gmm_independence_graph() returns the undirected second-order
conditional-independence (Gaussian graphical model) structure of a fitted
mixture -- the partial-correlation skeleton of its overall covariance, in
closed form. Composed with fit_kld_em(), it recovers the dependency
structure of a target that can only be evaluated (an unnormalised energy /
Gibbs density), where no sample exists to drive a sampling-based estimator.
It is a graphical-model diagnostic, not a causal-discovery method: it recovers
the undirected Markov skeleton, not edge directions.gmm() now rejects a covariance that is not positive-definite. A component
covariance with a negative variance (or, more generally, a negative eigenvalue
beyond numerical tolerance) was previously accepted silently, deferring the
failure to a later evaluation. The validator now checks the smallest eigenvalue
of each finite covariance against a scale-relative tolerance, so a clearly
indefinite matrix is caught at construction while a numerically near-singular
fit still constructs.gmm_impute(mechanism = censored(...)) is now numerically stable when a mixture
component sits well past the censoring bound. The fit could previously produce
NA responsibilities and abort, because it initialised the missing entries from
the observed mean -- the wrong side of a one-sided bound -- and the
truncated-conditional moments could overshoot the interval. It now initialises
inside the censored interval and clamps the truncated moments to it.gmm_impute() now takes a
missingness mechanism. mnar() specifies a selection model in which an entry of
a coordinate is missing with a probability that depends on its own unobserved
value, fitted jointly with the mixture so the imputations are not biased by an
ignorable model. The slope is a sensitivity parameter, supplied rather than
estimated, and the intercept is calibrated to the observed missingness rate.
censored() handles a known interval, such as a detection limit, by drawing
from the mixture conditional truncated to that interval in closed form. mar()
is the default and reproduces the earlier behaviour.proxy_mnar_sensitivity() sweeps the sensitivity slope over a grid and
pools the coordinate mean at each value, the input to a tipping-point analysis.gmm_eos_test() asks whether the last
m observations of a series are consistent with a linear-Gaussian state-space
model fitted on the rest, in the small-m regime (even m = 1) where ordinary
structural-break tests are undefined because the post-break parameters cannot
be estimated. It scores the last m standardised one-step filter innovations
and calibrates the score either parametrically (a chi-square reference, exact
under Gaussian innovations) or by the distribution-free subsampling test of
Andrews (2003), which stays valid under heavy-tailed observation noise. A new
vignette, Testing the last observation for instability, works through a
local-level example.Multiple imputation by Gaussian-mixture conditioning. gmm_impute()
fits a Gaussian mixture to a numeric dataset that contains missing values
and draws m completed datasets from the mixture conditional
p(x_missing | x_observed), the same Schur-complement algebra as
gmm_conditionalise(). Because the mixture can be multimodal and
heteroscedastic, the imputations follow the shape of the joint
distribution, which keeps downstream inference valid on data a
single-Gaussian or linear-Gaussian imputer mis-specifies. The mixture is
fitted to the incomplete data by an expectation-maximisation that uses each
row's observed margin and restores the conditional covariance of the filled
entries; each completion is drawn under a mixture fitted to a bootstrap
resample of the rows, so the pooled inference reflects both imputation and
parameter uncertainty. This release covers numeric data missing at random.
Pooling, diagnostics, and mice interoperability. gmm_complete()
extracts the completed datasets. proxy_pool() pools a column mean in
closed form -- the exact large-sample limit of the between-imputation
variance, with no Monte-Carlo noise and an imputation / parameter variance
split -- and proxy_fmi() reports its fraction of missing information. For a
regression or any other model estimand, as_mids() packages the completions
as a mice object so the joint mixture imputations flow into mice::pool()
unchanged: proxymix supplies the imputation model, mice the pooling. A new
vignette, Imputing missing data with a mixture, works through a multimodal
example.
A mixture proxy for the optima of an objective. from_objective()
fits a Gaussian-mixture proxy to the Gibbs measure exp(-f(x) / T) of a
user-supplied objective f over a bounded box, by cooling a short
temperature ladder through regime-(iii) KLD-EM (fit_kld_em()). The Gibbs
measure can be evaluated point-wise but not directly sampled, so this is
regime (iii) applied to an objective: the returned mixture is a closed-form
map over the low regions of f, and a multimodal objective is recovered as
a whole rather than one optimum at a time. A new vignette, Mapping the
optima of an objective, works through the bimodal and Himmelblau cases.
Modes of a Gaussian mixture. gmm_modes() returns the distinct local
modes of a mixture density by Gaussian mean-shift from each component mean
(Carreira-Perpinan 2000), with the mixture density at each mode. It
resolves the map from from_objective() into the recovered optima and
applies to any gmm or gmm_fit.
Closed-form entropy diagnostics. gmm_entropy() returns the closed-form
quadratic (order-2) Renyi entropy of a Gaussian mixture, or a Monte-Carlo
Shannon estimate bracketed by an analytic upper bound. gmm_divergence()
returns the closed-form, symmetric Cauchy-Schwarz divergence between two
mixtures (type = "cs"), with type = "kl" delegating to gmm_kld(). Both
rest on the Gaussian-product integral, so the closed-form quantities are exact
finite sums of Gaussian-density evaluations. A new vignette,
Entropy diagnostics with proxymix, demonstrates the layer.
Mutual information and conditional predictive entropy.
gmm_mutual_information() returns the closed-form, non-negative Cauchy-Schwarz
mutual information between two coordinate blocks of a fitted joint (the
divergence between the joint and the product of the marginals; zero exactly
under independence). gmm_conditional_entropy() returns the order-2 Renyi
entropy of the conditional mixture from gmm_conditionalise() -- the
predictive uncertainty of the target coordinates given the conditioned ones,
evaluated row-by-row.
Deterministic-annealing fitting and phase-transition component discovery.
fit_em_samples() and fit_kld_em() gain an opt-in anneal = TRUE argument
that locates the mixture components by deterministic annealing -- cooling a
temperature from a high value toward one -- before the unchanged cold EM loop
polishes the fit. The annealed warm-start is markedly less sensitive to local
optima than a cold multi-start. The companion diagnostic gmm_anneal_path()
tracks the number of distinct centroids as the temperature falls, a
physics-derived component-count read whose first phase transition has the
closed-form critical temperature lambda_max(Sigma^{-1} C).
Maximum-entropy targets and the ICL criterion. maxent_target()
constructs the least-committal gmm_target consistent with the supplied
constraints: the Gaussian under second-moment constraints on full support, the
uniform under a support constraint alone, and a truncated Gaussian under
second-moment constraints on a box (fit via regime (iii) under the
automatically selected support-matched proposal). bic_aic() now also returns
the integrated completed likelihood icl and the classification_entropy of
the fitted responsibilities (Biernacki, Celeux and Govaert 2000), which
penalises overlapping components and equals the BIC for a single component.
Kernel regression in the unifying-primitive vignette. The
One mixture, many methods vignette now covers Nadaraya-Watson kernel
regression as the nonparametric end of the conditioning axis: the conditional
mean of a one-component-per-datum kernel density estimate equals the
Nadaraya-Watson estimator exactly, so a single conditioning operation spans
ordinary least squares (K = 1) through fully-local kernel smoothing
(K = n).
Kalman filtering over time in the operator-calculus vignette. The
Affine-Gaussian operator calculus vignette now shows the predict
(gmm_affine) and update (gmm_observe) operators run as a filter over a
time series: at one component the recursion is exactly the classical Kalman
filter (verified against a textbook implementation on a constant-velocity
track), and at several components it is the Gaussian-sum filter.
Mixture reduction. gmm_reduce() collapses a Gaussian mixture to a budget
of at most k_max components. The default method = "merge" is a greedy,
moment-preserving pairwise merge, using either the Runnalls (2007)
Kullback-Leibler bound (cost = "kl") or a closed-form Cauchy-Schwarz cost
(cost = "cs"); every merge preserves the combined weight, mean and
covariance, so the reduced mixture has the same global mean and covariance as
the original, and reducing to one component returns the moment-matched
Gaussian. method = "anneal" refines the merge with an annealed re-fit of a
budget-sized proxy and keeps it when it improves on the merge (never worse),
which helps for smooth, over-parameterised mixtures. Reduction bounds the
component count of a Gaussian-sum filter built from gmm_affine() and
gmm_observe().
Bounded Gaussian-sum filtering. gmm_filter() runs a filter over an
observation series by alternating the predict operator (gmm_affine()), the
update operator (gmm_observe()) and an optional reduction (gmm_reduce()).
At one component it is the Kalman filter (verified against a textbook
implementation); with a Gaussian-sum process or measurement noise -- a gmm
supplied in place of a covariance matrix -- it is the Gaussian-sum filter of
Alspach and Sorenson (1972), and the k_max cap holds the component count at
budget over a long horizon. It returns the filtered mixture at each step, the
filtered means and covariances, and a tidy per-step summary including the
log marginal evidence. Constant and time-varying dynamics and measurements are
both supported. A new section of the Affine-Gaussian operator calculus
vignette demonstrates the verb.
(outcome, treatment, covariates) is read
-- in closed form, from that single fit -- as prediction, heterogeneous
treatment effects, optimal per-unit actions, off-line policy value, and an
identification audit. fit_uplift() assembles the joint fit and returns an
uplift_model; the proxy_* verbs score it without re-fitting:
proxy_predict() (the response / risk-scoring rung), proxy_cate() /
proxy_uplift() (heterogeneous effects with a delta-method or resampling
standard error), proxy_decide() (the revenue-maximising action plus an
action-flip probability), proxy_policy_value() (off-line value of a
targeting policy), proxy_confounding_gap() (the sensitivity to a latent
confounder), proxy_retrospective_uplift() (counterfactual-mean uplift for
observed units), proxy_regime_segments() (the fitted regimes as a segment
table), proxy_overlap() (per-unit positivity / mass coverage), and
proxy_identification_report() (the executive one-pager). The decision layer
rides only identified quantities -- the conditional average treatment effect
and counterfactual means.gmm_intervene() returns the interventional law p(. | do(.), .) -- it sets
the intervened coordinates inside every component without re-weighting the
regime gate (the graph surgery that distinguishes do(T = t) from T = t).
gmm_counterfactual() returns the K-atom counterfactual law of one observed
unit by abduction, action, and prediction. Only the counterfactual mean is
identified: the new gmm_counterfactual_law object exposes gmm_cf_mean(),
while gmm_cf_variance() and gmm_cf_tail_prob() deliberately error
(proxymix_not_identified) -- the individual counterfactual law depends on an
unidentified cross-world coupling.fit_uplift(outcome_type = "binary") fits on the latent
continuous scale; proxy_cate(scale = "response") and
proxy_predict(scale = "response") report effects and predictions on the
discretised predictive probability P(Y > threshold). Count outcomes are
supported on the same latent-scale reading.lm), clustering (kmeans / mclust), principal
components (prcomp), and ridge regression, with the trade-off of each
substitution stated. At N = 1 the conditional mean is exactly least squares,
the covariance eigenvectors are exactly the principal components, and a
covariance ridge is exactly the L2 penalty.data.table import for the decision-verb return tables.test-intervene.R), the K = 1 reduction of proxy_cate() to the
ordinary-least-squares treatment coefficient and the asymptotic agreement of
the delta standard error with the lm coefficient standard error
(test-uplift-cate.R), the audit verbs including recovery of a planted latent
confounder (test-uplift-audit.R), and a six-process synthetic validation
battery graded against known ground truth (test-uplift-validation.R).gmm_target() gains an optional support argument
(list(lower = , upper = ), with -Inf / Inf for unbounded
coordinates). When a regime-(iii) fit is given such a target and no
explicit proposal, fit_kld_em() now selects a support-matched
is_uniform() proposal automatically -- inset inside a compact box, or
data-derived for a one-sided coordinate -- instead of the default
multivariate-t, which placed importance mass where the target log-density
is -Inf and produced non-finite weights. The automatic choice is
announced with a one-line message, never silently. Unbounded targets are
unaffected: they keep the heavy-tailed default.epanechnikov_target(). A compact-support fixture (the
Epanechnikov kernel (3/4)(1 - u^2) on a box) joining banana_target()
/ donut_target() / mixture_target(). It declares its support, so it
fits via regime (iii) under the auto-selected uniform proposal with no NaN
weights -- the canonical case where no mixture of full-support Gaussians
can have compact support.inst/validation/regime_iii_pinned_fits.R gains a pinned Epanechnikov
bounded-support fit (ESS, support fraction, no NaN weights).test-support-aware-proposal.R), the exact K = 1
conditional-mean and conditional-variance match against lm
(test-gmr-k1-lm.R), and class / constructor / diagnostic contract
branches. Line coverage raised to >= 90%.New: a consumer seam for external Bayesian posteriors.
from_fb_posterior() takes a posterior addressed only through its
(unnormalised) log-density and returns a closed-form Gaussian-mixture
proxy via importance-sampled KLD-EM, generalising the input source from a
kernel-density estimate (from_kde()) to any Bayesian posterior. The
producer interface the seam expects is materialised and documented by
fb_log_posterior_spec(); fb_producer_available() is a capability probe
(degrades to FALSE, never errors, when no producer package is
installed); and mock_fb_posterior() is a synthetic producer (known
Gaussian or banana log-density) for testing the path with no producer
package present. proxymix never Imports: a producer package -- the seam
is a soft contract and R CMD check is clean with none installed.
New: autoplot() method for gmm_fit. Render a fitted proxy with
ggplot2::autoplot(fit) — a marginal density curve in one dimension, or a
viridis density raster with per-component ellipses in two. Any ambient
dimension is supported: the requested coordinates are reduced through the
package's own closed-form gmm_marginalise() before plotting (e.g.
autoplot(fit, dims = c(1L, 3L))). ggplot2 stays an optional dependency
— the method registers only when ggplot2 is installed, so R CMD check
remains clean with no sibling package present.
Second methodological extension. Brings a complete affine-Gaussian operator calculus to Gaussian-mixture proxies — pushforward, Bayesian update on a noisy linear observation, aggregation, missing-data conditioning — each closed-form and component-wise.
gmm_affine(g, A, b, noise_cov) — closed-form pushforward
of a Gaussian mixture through y = A x + b + epsilon,
epsilon ~ N(0, noise_cov). Returns the mixture in R^m with
mu'_k = A mu_k + b, Sigma'_k = A Sigma_k A' + noise_cov,
weights unchanged.gmm_observe(g, A, y, noise_cov) — Bayesian update on a
noisy linear observation. Applies the Kalman gain per component and
reweights component weights by per-component marginal evidence. The
finite-mixture analogue of a Kalman update.gmm_aggregate(g, A, noise_cov) — named alias for
gmm_affine() aimed at downscaling / aggregation pipelines.gmm_missing(g, observed, values) — Schur-complement
conditioning routed through an integer-index API for missing-data
pipelines.docs/design/operator_calculus_v0.3.md — pre-implementation design
note: maths, scope constraints (no non-affine fallbacks, no
approximate closed form), public API freeze, test obligations,
performance budget, pre-release gate.vignettes/operator_calculus.Rmd — educational vignette with
Kalman parity check, sequential vs stacked observations,
aggregation through a coarsening matrix, and a comparison to a
Gaussian-process latent.inst/validation/operator_calculus_pinned.R — three pinned
reference pipelines (Kalman parity, sequential vs stacked,
aggregate-then-observe) with hand-coded acceptance ranges.test-operator-calculus.R (12 tests, 46 expectations): A0–A2
(affine of moments), O0–O2 (Kalman parity, vanishing-evidence
guard, Bayes consistency), G0 (aggregate alias), M0 (missing vs
conditionalise), C0 (composition with marginalise), plus full
input-validation coverage.R/operator_calculus.R consolidates the four operators with
shared validation helpers (.validate_A, .validate_b,
.validate_noise_cov) and a single numerical-hygiene policy
(ridge after each output covariance, symmetrisation, chol-based
inverse with retry on near-singular matrices).gmm_observe() issues a proxymix_observe_no_update warning when
the marginal evidence is numerically zero at every component and
returns the prior unchanged with a metadata flag.Two methodological extensions that compose cleanly with regime (iii):
from_kde(). Compiles a kernel density
estimate over an n by p sample matrix into a closed-form
Gaussian-mixture proxy via regime-(iii) KLD-EM. Supports scalar and
diagonal bandwidths ("silverman", "scott", numeric scalar, or
per-coordinate numeric vector). Dimensional guard: p <= 5
recommended, p <= 10 allowed with warning, p > 10 rejected. The
KDE-target is normalised = TRUE by construction, so downstream
KLD and Hellinger diagnostics report absolute values. Companion
vignette: vignettes/from_kde.Rmd.gmm_target_from_posterior() (Contract A constructor).
S3 generic that compiles an (unnormalised) Bayesian posterior into a
gmm_target. The function method accepts a bare vectorised callable
with required parameter_names; the default method points users at
either a registered Bayesian-package method (brms,
Stan, ...) or the function-based path. Vectorisation contract is
enforced at construction by a probe call.DESCRIPTION now ships the canonical GitHub
namespace at github.com/max578/proxymix.test-from-kde.R (8 tests, 25 expectations): end-to-end recovery,
bandwidth selection branches, dimensional guards, normalisation
declaration, default proposal sanity, metadata pass-through.test-from-posterior.R (7 tests, 21 expectations): vectorisation
contract enforcement, log-normalizer pass-through, default-method
hinting, name validation, round-trip through fit_proxymix(regime = "kld"), attribute-based parameter-name support.inst/validation/from_kde_pinned_fits.R: pinned validation across
three reference KDE -> GMM pipelines (bimodal, banana, mixture)
with MC-SE-aware acceptance ranges.vignettes/from_kde.Rmd: educational walk-through covering scope,
bandwidth sensitivity, recovery on a known mixture, and the
contrast between KDE and proxy log-densities.gmm_target_from_posterior registers an S3 generic, paving the way
for class-specific methods registered by Bayesian fitting packages
(brms, Stan, pymc-via-reticulate) without coupling proxymix to
any specific Bayesian backend.Scientific hardening pass: regime (iii) is made harder to misuse.
gmm_target gains two new properties,
normalised (logical or NA) and log_normalizer (numeric or NA),
so that an unnormalised log_density can be supplied without making
downstream KLD or Hellinger diagnostics misleading. All three built-in
targets (banana_target(), donut_target(), mixture_target())
declare normalised = TRUE; the unnormalised case is now explicitly
documented at the target level.gmm_canonicalise() function
reorders the components of a gmm (or gmm_fit) by descending
weight, then by descending ||mu|| as a tiebreaker. fit_proxymix()
and the regime-specific fitters now canonicalise their outputs by
default (canonicalise = TRUE), making prints, snapshot tests, and
cross-run comparisons reproducible. Set canonicalise = FALSE to
retain the raw EM-order parameters.fit_kld_em() (and
therefore fit_proxymix(regime = "kld")) accepts validation_size
and validation_proposal. When validation_size > 0, a second
independent IS sample is drawn and the fit's diagnostics list
records validation_kld, validation_ess, and
validation_max_weight. This lets users tell the difference between
in-sample overfit and a fit that generalises across IS draws.fit_kld_em() now records
ess_relative (ESS / is_size), max_weight (largest
self-normalised weight), support_fraction (fraction of IS draws
with finite log-density under target and proposal), and a
Monte-Carlo standard error for the final KLD estimate
(mc_se_kld). A new ess_summary() helper returns the headline
numbers as a small list.kld_is_shifted and kld_shift_explanation whenever the target is
unnormalised or its normalisation is unknown, so users do not silently
read a shifted MC integral as an absolute divergence.hellinger_mc() now warns when the target is
not declared normalised = TRUE — the squared Hellinger distance is
not meaningful against an unnormalised target.fit_kld_em() issues a cli warning
when more than 5% of importance-sample draws fall outside the
proposal's support or carry non-finite weights. The most common
trigger is an is_uniform() proposal whose box does not cover the
target's mass.inst/validation/regime_iii_pinned_fits.R — a runnable validation
script that fits the three built-in targets with pinned seeds and
records final KLD, ESS, max weight, validation KLD, and runtime;
intended as the seed of a growing inst/validation/ corpus.test-canonicalise.R, test-normalisation.R,
test-validation-split.R, test-support-warning.R, and
test-monotone-objective.R. The last asserts monotonicity of the
fixed IS-weighted objective (\sum_n W_n \log g_\theta(x_n)) under
exact KLD-EM updates, which is a tighter check than the previous
generic "trace decreases" test.gmm_canonicalise() is the single source of truth for component
ordering — used by all three fitters and the dispatcher.Initial development release. Local-only; not yet on CRAN.
fit_proxymix() top-level dispatcher with three fitting regimes:
"moment" (closed-form moment matching), "sample" (classical EM on
i.i.d. samples), and "kld" (importance-sampled KLD-EM against an
evaluable-only target density). The "auto" regime picks the cheapest
applicable regime from the structure of the supplied gmm_target.gmm_target, gmm_fit, is_proposal, with
print() / format() methods and validators.gmm_ops.R: dgmm(), rgmm(),
gmm_marginalise(), gmm_conditionalise() (Schur complement),
gmm_kld() (Monte Carlo estimator with variational upper / lower
bounds for sanity).proposals.R: is_uniform(),
is_mvn(), is_mvt(); all wrap a is_proposal instance.kld_trace(), ess_trace(), hellinger_mc(),
bic_aic().init.R,
plus init_random(), init_kmeans(), init_moment_seed().banana_target(),
donut_target(), mixture_target(), plus the from-samples and
from-function constructors gmm_target_from_samples() and
gmm_target().quickstart, three_regimes, density_shapes
(the regime-(iii) demonstration), and roadmap (planned interfaces).The following functions ship with stable signatures, full documentation,
and signature-stability tests; the body raises a "not yet implemented"
condition with a pointer to vignettes/roadmap.Rmd.
from_kde() — KDE to GMM proxy via KLD-EM.from_aggregate_likelihood() — aggregate-likelihood downscaling
(Sejdinovic et al. kernel-downsizing framework).fit_kld_em_collider() — KLD-EM under DAG-implied conditional
independence constraints.to_apsim_scenarios() — Gaussian-mixture samples to APSIM scenario
tables.from_simulator() — wrap an expensive simulator as a
gmm_target via kernel-density or empirical-likelihood
bridges.Adaptive importance sampling, variational boosting, normalising-flow proposals, Stan / INLA inter-operation.