Addressing the problem
Why some quantities are exact
The Shannon entropy of a density \(f\) is the average of \(-\log f(x)\) over draws from \(f\). The Rényi-2 entropy is \(-\log \int f(x)^2\, dx\). Both are measured in nats, the unit of natural logarithms. Write \(\phi(x; m, S)\) for the normal density with mean \(m\) and covariance matrix \(S\). The integral of a product of two normal densities has an exact value: \(\int \phi(x; a, A)\, \phi(x; b, B)\, dx = \phi(a; b, A + B)\). The square of a mixture, or the product of two mixtures, is a sum of such products, and its integral is therefore exact too. The Rényi-2 entropy, the Cauchy-Schwarz divergence and the Cauchy-Schwarz mutual information are built from these integrals. The Shannon entropy needs the logarithm of a sum of densities, which has no exact formula.
Entropy of a mixture
The code below builds a mixture of two normal distributions with centres 4 apart, and a single normal distribution with correlation 0.3. For a single normal distribution in \(p\) variables with covariance matrix \(\Sigma\), the Rényi-2 entropy is \(\tfrac{p}{2}\log(4\pi) + \tfrac{1}{2}\log\det\Sigma\), an independent check on the package.
g <- gmm(
weights = c(0.5, 0.5),
means = list(c(-2, 0), c(2, 0)),
covariances = list(diag(2), diag(2))
)
h2_mixture <- gmm_entropy(g)
sigma_one <- matrix(c(1, 0.3, 0.3, 1), 2L, 2L)
one <- gmm(weights = 1, means = list(c(0, 0)), covariances = list(sigma_one))
h2_closed <- gmm_entropy(one)
h2_analytic <- 0.5 * (2 * log(4 * pi) +
as.numeric(determinant(sigma_one, logarithm = TRUE)$modulus))
h2_gap <- abs(h2_closed - h2_analytic)
sh <- gmm_entropy(g, order = "shannon", n_mc = 5000L, seed = 1L)
sh_slack <- sh$upper_bound - sh$mc
sh_slack_se <- sh_slack / sh$mc_se| Quantity | Value (nats) |
|---|---|
| Rényi-2, two-component mixture | 3.2060 |
| Rényi-2, single normal, from the package | 2.4839 |
| Rényi-2, single normal, from the formula | 2.4839 |
| Shannon, two-component mixture, simulation estimate | 3.4779 |
| Shannon, standard error of the estimate | 0.0134 |
| Shannon, exact upper bound | 3.5310 |
Distance between two mixtures
The Cauchy-Schwarz divergence is symmetric, and zero only when the
two mixtures are equal. The Kullback-Leibler divergence is not symmetric
and has no exact formula for mixtures. With type = "kl",
gmm_divergence() returns a simulation estimate and the
approximation of Hershey and Olsen (2007), which needs no
simulation.
q <- gmm(
weights = 1, means = list(c(0, 0)), covariances = list(diag(2) * 2)
)
d_cs <- gmm_divergence(g, q)
d_self <- gmm_divergence(g, g)
d_kl <- gmm_divergence(g, q, type = "kl", n_mc = 2000L)| Quantity | Value (nats) |
|---|---|
| Cauchy-Schwarz, g against q | 0.3881 |
| Cauchy-Schwarz, g against itself | 0.0000 |
| Kullback-Leibler, simulation estimate | 0.5703 |
| Kullback-Leibler, standard error of the estimate | 0.0218 |
| Kullback-Leibler, Hershey-Olsen approximation | 0.5003 |
Dependence between variables
The Cauchy-Schwarz mutual information compares the joint distribution
with the distribution the variables would have if they were independent.
Here the conditional entropy of \(x_1\)
at a value of \(x_2\) is the Rényi-2
entropy of \(x_1\) when \(x_2\) is held at that value. It is computed
at each value separately, not averaged over \(x_2\).
gmm_conditional_entropy() takes one row per value, with
NA marking the free variable.
sigma_joint <- matrix(c(1, 0.7, 0.7, 1), 2L, 2L)
joint <- gmm(
weights = 1, means = list(c(0, 0)), covariances = list(sigma_joint)
)
mi <- gmm_mutual_information(joint, 1L, 2L)
independent <- gmm(
weights = 1, means = list(c(0, 0)), covariances = list(diag(2))
)
mi_independent <- gmm_mutual_information(independent, 1L, 2L)
given_grid <- rbind(c(NA, 0), c(NA, 1), c(NA, 2))
h_cond <- gmm_conditional_entropy(joint, given = given_grid)| Quantity | Value (nats) |
|---|---|
| mutual information, correlation 0.7 | 0.1030 |
| mutual information, independent variables | 0.0000 |
| conditional entropy of \(x_1\) at \(x_2 = 0\) | 0.9288 |
| conditional entropy of \(x_1\) at \(x_2 = 1\) | 0.9288 |
| conditional entropy of \(x_1\) at \(x_2 = 2\) | 0.9288 |
Cooling the fit
The EM algorithm, the usual fitting method for a mixture, repeats two
steps: it shares each data point among the components, then refits each
component to its share. Deterministic annealing (Rose, 1998) softens the
shares with a temperature \(T\): point
\(i\) goes to component \(k\) in proportion to \(\pi_k\, \phi(x_i;\ \mu_k,
\Sigma_k)^{1/T}\), where \(\pi_k\) is the weight of the component. At
a high temperature every component sits at the mean of the data. As
\(T\) falls towards one, the components
split apart at critical temperatures. The quantity minimised is the free
energy \(F = \langle E \rangle - T H\),
where \(\langle E \rangle\) is the
average of \(-\log \phi(x_i;\ \mu_k,
\Sigma_k)\) over the shares. \(H\) is the entropy of the shares relative
to the component weights, \(H = -n^{-1}
\sum_{i,k} \gamma_{ik} \log(\gamma_{ik} / \pi_k)\), where \(\gamma_{ik}\) is the share of point \(i\) given to component \(k\). It is zero when every point is shared
in proportion to the weights. With anneal = TRUE,
fit_em_samples() and fit_kld_em() cool first
and start the usual fit from the result. The data below are three
far-apart clusters of 100 points each.
x_three <- rbind(
matrix(rnorm(200L), ncol = 2L) +
matrix(rep(c(-7, -7), each = 100L), ncol = 2L),
matrix(rnorm(200L), ncol = 2L) +
matrix(rep(c(7, -7), each = 100L), ncol = 2L),
matrix(rnorm(200L), ncol = 2L) +
matrix(rep(c(0, 8), each = 100L), ncol = 2L)
)
tgt_three <- gmm_target_from_samples(x_three)
fit_annealed <- fit_em_samples(tgt_three, N = 3L, anneal = TRUE, seed = 1L)
fit_annealed@diagnostics$annealed
#> [1] TRUEgmm_anneal_path() counts the distinct component centres
at each temperature. It returns the count that held for the largest
number of cooling steps, which are evenly spaced in log temperature. The
temperature of the first split can be computed from the covariance
matrix of the data, \(C\). It is \(T_c = \lambda_{\max}(\Sigma^{-1} C)\), the
largest eigenvalue of \(\Sigma^{-1}
C\), where \(\Sigma = \sigma^2
I\) is the covariance matrix of every component during cooling.
With the default \(\sigma = 1\), \(T_c\) is the largest eigenvalue of \(C\), which is the variance of the data
along the direction in which they are most spread out.
path <- gmm_anneal_path(x_three, k_max = 6L, n_steps = 60L, seed = 1L)
k_found <- path$k_selected
t_empirical <- path$first_critical_temperature
t_analytic <- path$t_critical_analytic| Quantity | Value |
|---|---|
| components found | 3 |
| first critical temperature, recorded | 47.60 |
| first critical temperature, exact | 51.62 |
| cooling steps | 60 |
anneal_df <- rbind(
data.frame(
temperature = path$path$temperature,
value = path$path$n_effective,
panel = "distinct centres"
),
data.frame(
temperature = path$path$temperature,
value = path$path$free_energy,
panel = "free energy"
)
)
ggplot2::ggplot(anneal_df, ggplot2::aes(temperature, value)) +
ggplot2::geom_vline(
xintercept = t_analytic, linetype = "dotted",
colour = "#0072B2", linewidth = 0.7
) +
ggplot2::geom_vline(
xintercept = t_empirical, linetype = "dashed",
colour = "#D55E00", linewidth = 0.7
) +
# a count recorded at one temperature holds until the next, cooler step
ggplot2::geom_step(direction = "vh", linewidth = 0.8,
colour = "#000000") +
ggplot2::facet_wrap(~ panel, ncol = 1L, scales = "free_y") +
ggplot2::scale_x_log10() +
ggplot2::labs(
x = "temperature (log scale, cooling from right to left)",
y = NULL,
title = "Where the components split as the fit cools"
) +
ggplot2::theme_minimal(base_size = 11)
Cooling on three well-separated clusters. Upper panel: the number of distinct component centres at each temperature. Lower panel: the free energy. The dashed line marks the temperature at which the first split was recorded, the dotted line the exact critical temperature.
Maximum-entropy targets
Among all densities that meet a set of constraints, the one with the
largest entropy assumes nothing beyond them (Jaynes, 1957).
maxent_target() builds it. For a given mean and covariance
matrix it is the normal distribution, for a given range alone the
uniform, and for a mean and covariance matrix on a bounded box a
truncated normal. A bounded target records its range. When such a target
is fitted from its formula alone (regime = "kld", van der
Hoek and Elliott, 2024), the trial points that the fit weights are drawn
from within that range.
me_gauss <- maxent_target(moments = list(mean = c(0, 0), cov = diag(2)))
me_unif <- maxent_target(support = list(lower = c(0, 0), upper = c(1, 1)))
unif_density <- exp(me_unif@log_density(matrix(c(0.5, 0.5), nrow = 1L)))
gauss_density <- exp(me_gauss@log_density(matrix(c(0, 0), nrow = 1L)))| Constraint | Family returned | Density at the centre |
|---|---|---|
| mean and covariance matrix | maxent_gaussian | 0.159 |
| the unit square as range | maxent_uniform | 1.000 |
Choosing the number of components
The BIC (Bayesian information criterion) and the AIC (Akaike
information criterion, with a lighter penalty) balance how well a
mixture fits against its number of parameters. In
bic_aic(), smaller is better. The ICL (Biernacki et al.,
2000) adds a penalty for overlap, \(\mathrm{ICL} = \mathrm{BIC} + 2 E_N\) with
\(E_N = -\sum_{i,k} \gamma_{ik} \log
\gamma_{ik}\), where \(\gamma_{ik}\) is the share of point \(i\) given to component \(k\). \(E_N\) is zero when every point belongs
clearly to one component and grows with overlap. The ICL equals the BIC
for one component, is never smaller than the BIC, and favours
well-separated components.
x_two <- rbind(
matrix(rnorm(200L, -4), ncol = 2L),
matrix(rnorm(200L, 4), ncol = 2L)
)
fit_two <- fit_em_samples(gmm_target_from_samples(x_two), N = 2L, seed = 1L)
crit <- bic_aic(fit_two)| Criterion | Value |
|---|---|
| BIC | 1492.24 |
| AIC | 1455.95 |
| ICL | 1492.24 |
| \(E_N\) | \(1.4 \times 10^{-13}\) |
| free parameters | 11 |
Comparison with FNN, mclust and pcalg
The FNN package estimates the Shannon entropy from the distances
between each data point and its nearest neighbours (Kozachenko and
Leonenko, 1987). Its mutinfo() estimates the mutual
information in the same way (Kraskov et al., 2004). Neither estimator
fits a model. In a simulation, 200 datasets of 500 rows were drawn from
each of two mixtures with known entropy. One has three overlapping
components in two variables. The other has two components in four
variables, and its mutual information is between the first two variables
and the last two. proxymix was given the true number of components,
which FNN does not need, and used 5000 simulation draws. FNN used 10
neighbours. The true values came from numerical integration of the true
density. The measure is the mean absolute error, and smaller is
better.
| Design | Shannon quantity | True value | Error, proxymix | Error, FNN | Difference (standard error) |
|---|---|---|---|---|---|
| three components, 2 variables | entropy | 3.533 | 0.039 | 0.056 | -0.017 (0.002) |
| two components, 4 variables | entropy | 6.230 | 0.066 | 0.074 | -0.009 (0.002) |
| three components, 2 variables | mutual information | 0.381 | 0.026 | 0.035 | -0.008 (0.001) |
| two components, 4 variables | mutual information | 0.396 | 0.028 | 0.031 | -0.003 (0.002) |
proxymix had the smaller error for the entropy in both designs and for the mutual information with two variables. The difference was 8.2 standard errors for the entropy with two variables, 4.0 for the entropy with four variables, and 5.6 for the mutual information with two variables. With four variables, the difference in mutual information is 1.3 standard errors, too small to separate the two methods.
Where the same simulation chose the number of components or recovered
a graph, other packages did better. With three overlapping components,
the BIC of the mclust package (Scrucca et al., 2016) chose the true
count in 100 per cent of datasets and proxymix’s BIC in 97 per cent.
Both ICLs mostly chose two components, and found the true count in 32
per cent of datasets for mclust and 2 per cent for proxymix. With two
components in four variables, mclust’s ICL found the true count in 100
per cent of datasets and proxymix’s in 83 per cent. On a six-variable
chain with two components, the PC algorithm of the pcalg package
(Kalisch et al., 2012), which removes edges by repeated tests of partial
correlations, was run at level 0.01 and recovered the graph exactly in
93.5 per cent of datasets. gmm_independence_graph(), at
level 0.05 per pair, did so in 63 per cent. The chain leaves 10 pairs
without an edge. Were their tests independent, all 10 would be left out
at level 0.05 with probability 0.60, close to proxymix’s 63 per cent, so
most of the gap comes from the looser level.
proxymix was slower than FNN and pcalg. On one dataset of 500 rows, the FNN estimates took 0.002 seconds and the proxymix fit with its estimates 0.029 seconds. The PC algorithm took 0.002 seconds and proxymix 0.165. proxymix’s five fits for the component count took 0.171 seconds and mclust’s ICL 0.318 (median of five runs on one computer).
The code below runs the three competitors and the proxymix estimates on one dataset from each design.
library(proxymix)
library(FNN)
library(mclust)
library(pcalg)
# one dataset of 500 rows from the design with three overlapping components
w <- c(0.4, 0.35, 0.25)
mu <- list(c(-2, 0), c(2, 1), c(0, 3))
sig <- list(matrix(c(1, 0.5, 0.5, 1), 2L),
matrix(c(1.5, -0.6, -0.6, 0.8), 2L),
diag(c(0.5, 1.2)))
set.seed(1L)
z <- sample.int(3L, 500L, replace = TRUE, prob = w)
x <- matrix(rnorm(500L * 2L), 500L, 2L)
for (k in 1:3) {
x[z == k, ] <- x[z == k, , drop = FALSE] %*% chol(sig[[k]]) +
matrix(mu[[k]], sum(z == k), 2L, byrow = TRUE)
}
# Shannon entropy and mutual information by nearest neighbours
entropy(x, k = 10L)[10L]
mutinfo(x[, 1L, drop = FALSE], x[, 2L, drop = FALSE], k = 10L)
# the same two quantities from a fitted three-component mixture
fit <- fit_em_samples(gmm_target_from_samples(x), N = 3L, seed = 1L)
h_mc <- function(g) {
gmm_entropy(g, order = "shannon", n_mc = 5000L, seed = 1L)$mc
}
h_mc(fit)
h_mc(gmm_marginalise(fit, 1L)) + h_mc(gmm_marginalise(fit, 2L)) - h_mc(fit)
# number of components by mclust's BIC and ICL, one to five offered
Mclust(x, G = 1:5, verbose = FALSE)$G
icl_m <- mclustICL(x, G = 1:5, verbose = FALSE)
which(icl_m == max(icl_m, na.rm = TRUE), arr.ind = TRUE)[1L]
# one dataset from the graph design: a six-variable chain, two components
omega <- diag(6L)
omega[abs(row(omega) - col(omega)) == 1L] <- -0.4
set.seed(1L)
z_g <- sample.int(2L, 500L, replace = TRUE, prob = c(0.5, 0.5))
x_g <- matrix(rnorm(500L * 6L), 500L, 6L) %*% chol(solve(omega))
x_g[z_g == 2L, 1L] <- x_g[z_g == 2L, 1L] + 4
# the graph by the PC algorithm, and by proxymix
pc_fit <- pc(suffStat = list(C = cor(x_g), n = nrow(x_g)),
indepTest = gaussCItest, alpha = 0.01, p = ncol(x_g))
adj <- as(pc_fit@graph, "matrix")
((adj + t(adj)) > 0) * 1L
sel <- select_N(gmm_target_from_samples(x_g), candidates = 1:4, seed = 1L)
gmm_independence_graph(sel$best_fit)It needs FNN and mclust from CRAN, and pcalg, which needs graph and RBGL from Bioconductor. It is not run when this vignette is built.
install.packages(c("FNN", "mclust", "BiocManager"))
BiocManager::install(c("graph", "RBGL"))
install.packages("pcalg")The extended version of this article gives the full code and results of this comparison, adds the graphical lasso and the huge package to the graph comparison, and applies each method to the Palmer penguins data.