Addressing the problem
Which fitting method applies
The package has three ways of fitting the bell curves, described by
van der Hoek and Elliott (2024). fit_proxymix() picks one
according to what you supply.
| What you have | What the package does | regime setting |
|---|---|---|
| Data, and you want one bell curve | Fits one bell curve with the same centre and spread as the data | "moment" |
| Data, and you want several bell curves | Moves and reshapes the bell curves, step by step, until together they match the data | "sample" |
| Only the formula, no data | Draws trial points, weights each one by the formula, then fits the bell curves to the weighted points | "kld" |
If you already have a sample from the target, the first two methods
apply, and established packages such as mclust,
mixtools and flexmix can also fit the mixture.
When you have only the formula, the third method applies. The package
was built for this case.
A target with no sampler
The banana target has an exact density formula and no sample, so only the third method applies.
Fit the proxy
The third method draws 2,000 trial points from a broad distribution that is easy to sample. Here it is a Student-t distribution, a relative of the normal with heavier tails, made wide enough to cover the banana. Each trial point is then weighted by how much more likely it is under the target than under the broad distribution. The weights do the same job as survey weights that correct an unrepresentative sample. The mixture is refitted to the weighted points in rounds, which stop when the fit no longer improves.
The closeness of the proxy to the target is measured by the Kullback-Leibler (KL) divergence. It is zero when the two distributions match and increases as they become more different.
The call below asks for a proxy with three components and sets a seed so the result is reproducible.
proposal <- proposal_mvt(n_dim = 2L, mean = c(0, 0),
sigma = 4 * diag(2), df = 5)
fit <- fit_proxymix(tgt, N = 3L, regime = "kld",
proposal = proposal,
is_size = 2000L,
max_iter = 60L,
seed = 1L)
fit
#> <gmm_fit>: regime = "kld", K = 3, p = 2
#> target : banana
#> iterations : 38
#> converged : TRUE
#> [1] w = 0.5345, |mu| = 0.3018, tr(Sigma) = 1.3502
#> [2] w = 0.3173, |mu| = 0.9255, tr(Sigma) = 2.0299
#> [3] w = 0.1483, |mu| = 1.5183, tr(Sigma) = 3.1396Check the fit before using it
Weighted draws have a weakness familiar from survey work. If a
handful of respondents carry very large weights, the survey estimate
rests on those few people and becomes unstable. The same happens here if
a few trial points carry most of the weight.
gmm_fit_quality() checks for this and for other signs of a
poor fit.
| Check | Value |
|---|---|
| fitting method | kld |
| rounds settled before the limit | TRUE |
| weights collapsed onto a few draws | FALSE |
| effective sample size | 820.6 |
| effective sample size as a share of all draws | 0.410 |
| smallest effective sample size of any component | 347.7 |
| largest share of the weight held by one draw | 0.00206 |
| KL divergence on the fitting draws | 0.00275 |
| half the variance of the log density ratio (a local KL approximation) | 0.0163 |
| KL divergence on fresh draws | 0.0126 |
| fresh draws minus fitting draws | 0.00983 |
The effective sample size is the number of equally weighted draws that the weighted sample is worth. The fit was tuned to the draws it was fitted on, so the KL computed on them is too low. The KL on a fresh set of draws is the one to report. The package flags a fit when this KL exceeds 0.3, when the fit did not converge, or when it is degenerate. Half the variance of the log density ratio on the fitting draws approximates the KL when the proxy is already close to the target. It is only a rough check.
Compare the proxy with the target
With two variables, the quickest check is a plot.
grid_x <- seq(-3, 3, length.out = 120L)
grid_g <- expand.grid(x1 = grid_x, x2 = grid_x)
grid_mat <- as.matrix(grid_g)
grid_g$target <- exp(tgt@log_density(grid_mat))
grid_g$proxy <- dgmm(grid_mat, fit)ggplot2::ggplot(grid_g, ggplot2::aes(x1, x2)) +
ggplot2::geom_contour_filled(ggplot2::aes(z = target), bins = 10L,
alpha = 0.85) +
ggplot2::geom_contour(ggplot2::aes(z = proxy), colour = "white",
linetype = "dashed", linewidth = 0.45, bins = 8L) +
ggplot2::scale_fill_viridis_d(option = "mako", guide = "none") +
ggplot2::coord_equal() +
ggplot2::labs(
title = "Target (filled) and fitted proxy (dashed)",
x = expression(x[1]), y = expression(x[2])
) +
ggplot2::theme_minimal(base_size = 11)
The banana target (filled contours) with the three-component proxy overlaid as dashed contours. The dashed contours follow the curve of the banana, which a single bell curve could not do. The KL divergence on fresh draws is 0.013.
Use the proxy
Questions that were hard to answer for the target have exact answers
for the proxy, because it is built from normal distributions.
gmm_marginalise(keep = 1L) gives the distribution of the
first variable on its own.
gmm_conditionalise(given = c(NA, 0.5)) gives the
distribution of the first variable when the second equals 0.5, with
NA marking the variable left free. Neither call goes back
to the target formula.
gmm_marginalise(fit, keep = 1L)
#> <marginalise(kld_em[N=3] on banana)>: K = 3 components in p = 1 dimensions
#> [1] w = 0.5345, |mu| = 0.1866, tr(Sigma) = 0.4417
#> [2] w = 0.3173, |mu| = 0.9236, tr(Sigma) = 0.4325
#> [3] w = 0.1483, |mu| = 1.3524, tr(Sigma) = 0.5179
gmm_conditionalise(fit, given = c(NA, 0.5))
#> <conditionalise(kld_em[N=3] on banana)>: K = 3 components in p = 1 dimensions
#> [1] w = 0.5597, |mu| = 0.2553, tr(Sigma) = 0.4338
#> [2] w = 0.3179, |mu| = 1.0803, tr(Sigma) = 0.2310
#> [3] w = 0.1224, |mu| = 1.2816, tr(Sigma) = 0.1535Drawing from the proxy is fast.
Comparison with the Laplace approximation, Stan and BayesianTools
In a simulation, proxymix was compared with three established methods
for a distribution that is known only by its formula. The Laplace
approximation (Tierney and Kadane, 1986) is a single normal distribution
centred on the highest point of the target, with a spread set by how
sharply the target falls away from that point. Stan (Carpenter et al.,
2017), run from R through cmdstanr, draws a sample from the
formula with the no-U-turn sampler, or NUTS (Hoffman and Gelman, 2014).
NUTS is a Markov chain Monte Carlo method: it produces random values by
a long chain of small, linked steps. A mixture was then fitted to the
NUTS draws with mclust (Scrucca et al., 2016). DEzs (ter
Braak and Vrugt, 2008), from BayesianTools (Hartig et al.,
2026), is another Markov chain Monte Carlo method. Each method was run
200 times on the banana target without further tuning. proxymix chose
its number of components automatically with select_N(), and
mclust chose its own from 1 to 6. Only this one target in
two variables was used, and the results may not carry over to targets
with more variables. The first measure is the KL divergence of each
fitted density from the target, computed by summing over a fine grid of
points. It cannot be computed for the two samplers, which return draws
but no density. The second measure is the error in the estimated
probability that \(x_1\) is greater
than 2, which is 0.0227 for the target. Smaller is better for both.
| Method | KL divergence, mean (sd) | Error in \(P(x_1 > 2)\) | Seconds per run |
|---|---|---|---|
| proxymix | 0.0081 (0.0031) | 0.00350 | 0.70 |
| NUTS draws, then mclust | 0.0139 (0.0046) | 0.00473 | 2.70 |
| Laplace approximation | 0.3749 (0.0000) | 0.00001 | < 0.01 |
| NUTS draws | – | 0.00460 | 1.04 |
| DEzs draws | – | 0.00546 | 0.53 |
proxymix gave the smallest KL divergence, 0.0081 against 0.0139 for the mixture fitted to the NUTS draws. The Laplace approximation was far behind on this measure, at 0.3749, because one bell curve cannot follow the curve of the banana. On the tail probability, however, the Laplace approximation was the most accurate, with an error of 0.00001. On its own, \(x_1\) has a standard normal distribution, and on this target the Laplace approximation reproduces it almost exactly. proxymix came second on the tail, with an error of 0.0035 against 0.0046 to 0.0055 for the two samplers and the mixture fitted to the NUTS draws. The held-out KL of 0.0126 reported for the three-component fit above is itself estimated from random draws, with a standard error of about 0.0014. On the grid used for the table, that fit has a KL divergence of 0.0154. That value lies 2.4 standard deviations above the proxymix mean in the table, using the standard deviation across runs shown there in brackets.
The Laplace approximation was also the fastest, at under 0.01 seconds
per run. DEzs took 0.53 seconds per run and was slightly faster than
proxymix at 0.70, while NUTS took 1.04 and NUTS followed by
mclust took 2.70 (medians on one computer, with other
programs running on it at the same time).
The code below runs each method once and scores it against the grid.
It needs mclust and BayesianTools from CRAN,
and cmdstanr from https://stan-dev.r-universe.dev with a CmdStan
installation. It is not run when this vignette is built.
library(proxymix)
library(cmdstanr)
library(mclust)
library(BayesianTools)
tgt <- banana_target()
# midpoint-rule quadrature; the target's mass outside the box is below 1e-5
h <- 0.05
grid <- as.matrix(expand.grid(x1 = seq(-5 + h / 2, 5, by = h),
x2 = seq(-5 + h / 2, 12, by = h)))
log_f <- tgt@log_density(grid)
f_grid <- exp(log_f)
tail_ref <- sum(f_grid[grid[, 1L] > 2]) * h^2
kl_of <- function(log_g) sum(f_grid * (log_f - log_g)) * h^2
# mass above x1 = 2 under a Gaussian mixture, from its x1 marginal
tail_of <- function(w, mean1, sd1) {
sum(w * pnorm(2, mean1, sd1, lower.tail = FALSE))
}
# the same target as a Stan program
stan_file <- file.path(tempdir(), "banana.stan")
writeLines(c(
"parameters {",
" vector[2] x;",
"}",
"model {",
" target += -0.5 * (square(x[1])",
" + square(x[2] - 0.5 * (square(x[1]) - 1)));",
"}"), stan_file)
banana_stan <- cmdstan_model(stan_file)
box <- createBayesianSetup(likelihood = function(x) tgt@log_density(x),
lower = c(-6, -6), upper = c(6, 14))
set.seed(1L)
fit <- select_N(tgt, seed = 1L)$best_fit
fit_1 <- gmm_marginalise(fit, keep = 1L)
la <- optim(c(0, 0), function(x) -tgt@log_density(x),
method = "BFGS", hessian = TRUE)
la_mean <- la$par
la_cov <- solve(la$hessian)
nuts <- banana_stan$sample(seed = 1L, refresh = 0L, show_messages = FALSE)
draws <- nuts$draws("x", format = "matrix")
mc <- Mclust(draws, G = 1:6, modelNames = "VVV", verbose = FALSE)
mc_par <- mc$parameters
de <- runMCMC(box, sampler = "DEzs", settings = list(message = FALSE))
de_draws <- getSample(de, start = 1000L)
# KL divergence from the target, for the methods that return a density
c(proxymix = kl_of(dgmm(grid, fit, log = TRUE)),
Laplace = kl_of(dmvnorm(grid, la_mean, la_cov, log = TRUE)),
"NUTS + mclust" = kl_of(dens(grid, mc$modelName, parameters = mc_par,
logarithm = TRUE)))
# error in the estimated probability that x1 > 2
c(proxymix = tail_of(gmm_weights(fit_1),
vapply(gmm_means(fit_1), `[[`, numeric(1L), 1L),
sqrt(vapply(gmm_covariances(fit_1), `[[`,
numeric(1L), 1L))),
Laplace = pnorm(2, la_mean[1L], sqrt(la_cov[1L, 1L]), lower.tail = FALSE),
NUTS = mean(draws[, 1L] > 2),
"NUTS + mclust" = tail_of(mc_par$pro, mc_par$mean[1L, ],
sqrt(mc_par$variance$sigma[1L, 1L, ])),
DEzs = mean(de_draws[, 1L] > 2)) - tail_refThe extended version of this article gives the full simulation, including how many settings each method needs the user to choose.