Parametric Bootstrap for Mixture-Order Selection

Why Bootstrap a Mixture-Order Comparison?

When comparing a mixture with \(K\) components to one with \(K+1\), the likelihood-ratio statistic measures the improvement in fit from allowing the additional component. Its familiar chi-square reference distribution is generally unreliable here: under the smaller model, the extra component can have zero weight, and its parameters are then unidentified. Mixture likelihoods can also have singularities unless the fitting problem is suitably constrained.

A parametric bootstrap offers a direct calibration strategy. Fit the smaller model, simulate data from that fitted model, and repeat the entire comparison on each simulated data set. The simulated statistics show how large a likelihood improvement can arise when data are generated under the fitted null model and both mixture orders are refitted. This is a plug-in approximation, not an automatic cure for model misspecification or unstable optimization.

Models and the Test Statistic

Let \(Y_1,\ldots,Y_n\) be independent observations in \(\mathbb{R}^p\). A \(K\)-component Gaussian mixture has density

g_K(y;\theta_K)=\sum_{j=1}^{K}\pi_j\phi_p(y;\mu_j,\Sigma_j), \qquad \sum_{j=1}^{K}\pi_j=1,\quad \pi_j\geq 0,

where \(\phi_p(\,\cdot\,;\mu_j,\Sigma_j)\) is the \(p\)-variate normal density with mean \(\mu_j\) and positive-definite covariance matrix \(\Sigma_j\). Write \(\ell_K(\theta_K)=\sum_{i=1}^n\log g_K(Y_i;\theta_K)\) for the log-likelihood, and let \(\hat\ell_K\) be its maximum over the specified \(K\)-component parameter space.

Definition: For comparing \(K\) components against \(K+1\), the likelihood-ratio statistic is \(T_{\mathrm{obs}}=2(\hat\ell_{K+1}-\hat\ell_K)\). If the smaller model is contained in the larger model and both maxima use the same likelihood and compatible constraints, then \(T_{\mathrm{obs}}\geq 0\).

For this containment, the \(K+1\)-component model must allow a component weight of zero: set the new weight to zero and retain the \(K\) original weights and component parameters. To make the maximization well-defined in practice, one may impose common restrictions, such as bounded means and covariance eigenvalues bounded away from zero and infinity. The restrictions must be compatible with that embedding. If the fitting procedure instead uses a penalty, define the statistic from the penalized objective and use that same objective throughout the bootstrap; it is then not the ordinary likelihood-ratio statistic above.

The Parametric Bootstrap Procedure

Fit the null model with \(K\) components to obtain \(\hat\theta_K\). For each bootstrap replicate \(b=1,\ldots,B\), generate \(n\) independent observations from \(g_K(\,\cdot\,;\hat\theta_K)\), refit both the \(K\)- and \(K+1\)-component models, and calculate

T_b^*=2(\hat\ell_{K+1,b}^*-\hat\ell_{K,b}^*).

The essential point is refitting both models in every replicate. Reusing the original estimates, or comparing the alternative fit on simulated data with the original null fit, does not reproduce the sampling behavior of the likelihood-ratio statistic.

Fit the null.
Maximize the \(K\)-component likelihood under the chosen constraints and save \(\hat\theta_K\) and \(\hat\ell_K\).

Simulate under the fitted null.
Generate \(n\) independent observations from \(g_K(\,\cdot\,;\hat\theta_K)\) for each replicate.

Refit both orders.
Apply the same constraints, convergence rules, and a comparably thorough initialization and restart strategy to both models on every replicate.

Compare the tail.
Count bootstrap statistics at least as large as \(T_{\mathrm{obs}}\), using ties as exceedances, and report the simulation count and fitting diagnostics.

A common finite-simulation estimate is

\hat p_{\mathrm{boot}}= \frac{1+\sum_{b=1}^{B} I(T_b^*\geq T_{\mathrm{obs}})}{B+1}.

The added one prevents a reported probability of zero when no simulated statistic exceeds the observed one. Because the null distribution is generated from an estimated parameter, this is a plug-in bootstrap calibration; the adjustment does not make it an exact test under every mixture null.

What the Bootstrap Calculation Guarantees

Result: Conditional on the observed data and fitted null parameter, if the bootstrap replicates are independent draws from the fitted null model and \(T_b^*\) is computed by the stated procedure, the empirical proportion of replicates satisfying \(T_b^*\geq T_{\mathrm{obs}}\) converges almost surely to the corresponding fitted-null tail probability as \(B\to\infty\).

Proof. Given the data, \(\hat\theta_K\) and \(T_{\mathrm{obs}}\) are fixed. Define \(I_b=I(T_b^*\geq T_{\mathrm{obs}})\). Under independent simulation and identical refitting rules, the indicators are independent and identically distributed Bernoulli random variables with success probability \(q=P_{\hat\theta_K}(T^*\geq T_{\mathrm{obs}})\), where \(P_{\hat\theta_K}\) denotes probability when the data are generated from the fitted null. They have finite expectation \(E[I_b\mid\text{data}]=q\). The strong law of large numbers therefore gives

\frac{1}{B}\sum_{b=1}^{B} I_b \xrightarrow{a.s.} q.

The plus-one estimate differs from this proportion by a quantity tending to zero: if \(S_B=\sum_b I_b\), then \(\left|(1+S_B)/(B+1)-S_B/B\right|=|B-S_B|/[B(B+1)]\leq 1/B\). Thus it has the same limit. This result concerns Monte Carlo approximation to the fitted-null tail probability; it does not establish that the fitted-null tail probability equals the true sampling tail probability. \(\square\)

The simulation error also matters. If \(S_B\) of \(B\) replicates exceed the threshold, then \(S_B/B\), conditional on the data, is an average of Bernoulli indicators. Its conditional variance is \(q(1-q)/B\). Replacing the unknown \(q\) by the estimated tail proportion gives the usual approximate Monte Carlo standard error:

\widehat{\operatorname{SE}}_{\mathrm{MC}} \approx \sqrt{\frac{\hat p(1-\hat p)}{B}}.

This measures simulation noise only. It does not account for uncertainty from estimating the null parameters, model misspecification, or local maxima in either fit.

Worked Bootstrap Comparisons

Worked Example: Computing a Bootstrap Tail Estimate

Suppose the constrained one-component fit has log-likelihood \(\hat\ell_1=-184.7\), and the constrained two-component fit has \(\hat\ell_2=-176.2\). The observed statistic is

T_{\mathrm{obs}}=2[-176.2-(-184.7)]=2(8.5)=17.0.

Generate \(B=999\) samples from the fitted one-component model and refit both orders on each. If 37 of the 999 bootstrap statistics are at least 17.0, then \(\hat p_{\mathrm{boot}}=(1+37)/(999+1)=38/1000=0.038\).

Using this estimated tail probability in the approximate Monte Carlo standard error gives \(\sqrt{0.038(0.962)/999}=\sqrt{0.0000365926}\approx 0.006049\). Thus the estimated tail probability has noticeable simulation uncertainty on the scale of several thousandths. Increasing \(B\) can reduce this simulation error, but cannot repair poor model fitting or an unsuitable null model.

Worked Example: Why Component Labels Do Not Matter

Consider a univariate two-component mixture with weights \(0.3,0.7\), means \(-1,2\), and variances \(1,4\). Its density is

g(y)=0.3\phi_1(y;-1,1)+0.7\phi_1(y;2,4).

If an optimizer returns the components in the opposite order, it reports weights \(0.7,0.3\), means \(2,-1\), and variances \(4,1\). Substitution gives \(0.7\phi_1(y;2,4)+0.3\phi_1(y;-1,1)=g(y)\) for every \(y\). Therefore each observation has exactly the same fitted density, and the log-likelihood is unchanged.

This permutation invariance means the likelihood-ratio calculation should compare maximized likelihoods, not component labels or component-by-component parameter differences. In a bootstrap replicate, whichever label order an optimizer returns, calculate the two fitted log-likelihoods and their difference. If parameters themselves must be compared, first align labels by minimizing a specified distance over all component permutations; that alignment is unnecessary for the likelihood-ratio statistic.

Worked Example: Handling an Unreliable Bootstrap Fit

Suppose \(T_{\mathrm{obs}}=5.0\), and four bootstrap data sets initially produce statistics \(4.8, 5.0, 7.1\), and an apparent failed fit. The failed fit is not a valid fourth statistic: a numerical routine that stops early may return a lower alternative likelihood and hence a misleadingly small statistic.

Inspect the convergence report and restart the two model fits for that replicate using the established initialization and restart protocol. Suppose the improved, converged statistic is \(4.2\). The four valid statistics are now \(4.8,5.0,7.1,4.2\); two satisfy the exceedance rule \(T_b^*\geq 5.0\). The plus-one estimate is \((1+2)/(4+1)=0.6\). This tiny simulation is only an arithmetic illustration, not a precise calibration.

If the fit still fails, record the failure and investigate rather than silently dropping it or counting it as a non-exceedance. Widespread failures can invalidate the calibration: improve optimization, revisit constraints, or report that the procedure did not yield a dependable bootstrap result.

Label Invariance, Model Checks, and Interpretation

Mixture components are not intrinsically ordered. Relabeling them leaves the density and likelihood unchanged, so ordering constraints introduced only for presentation must not change the set of densities being optimized. A consistent bootstrap uses the same parameter space for the original data and all replicates; otherwise, the observed and simulated statistics may not be comparable.

Before interpreting the tail estimate, check that:

The null and alternative use compatible constraints, and the null is genuinely embedded in the alternative.
Both orders are refit in every replicate with the same convergence criteria and a sufficiently robust restart strategy.
Likelihoods are finite, covariance restrictions are respected, and the alternative fit is no worse than the null fit beyond numerical tolerance.
Failures, boundary solutions, and unusually frequent near-zero component weights are investigated and reported.
The number of replicates is adequate for the tail probability being estimated; a small \(B\) gives coarse tail resolution and substantial Monte Carlo error.

The bootstrap calibrates the statistic under the fitted \(K\)-component Gaussian mixture, not under every distribution that might reasonably be called a \(K\)-component population. If that null model is a poor description of the data, the resulting tail probability may not answer the intended scientific question. The procedure also does not remove the singular likelihood problem: its fitting constraints and numerical checks remain essential.

Common Pitfall: Treating the Chi-Square Reference as Automatic

A tempting argument is that the alternative adds parameters, so the statistic must follow a chi-square distribution with degrees of freedom equal to the parameter-count difference. That regular likelihood-ratio argument requires conditions such as identifiable parameters and a nonsingular local information matrix. At a mixture-order null, the extra component may have zero weight, while its mean and covariance are then unidentified; the null is therefore singular in precisely the way that undermines the usual reference law.

The correction is not to assume that a bootstrap always fixes the problem. Instead, simulate from the fitted null and repeat the same constrained optimization to approximate the statistic's distribution under that fitted model. Then assess whether the fitting procedure is reliable and whether the null model is scientifically adequate. A bootstrap tail estimate is evidence conditional on those choices, not a universal replacement for checking assumptions.

Key takeaway: For mixture-order selection by a likelihood-ratio statistic, generate data from the fitted smaller mixture, refit both orders to every simulated data set, and compare likelihoods rather than component labels. Report the bootstrap tail estimate together with Monte Carlo uncertainty and fitting diagnostics.

Check Your Understanding

Try these questions before applying the procedure to a new mixture comparison.

Why must both the null and alternative mixture be refitted on every bootstrap data set?
What feature of a mixture-order null makes the usual chi-square likelihood-ratio argument unreliable?
If two component-label orderings produce the same density, what quantity should be used in the likelihood-ratio comparison?
With \(B=499\) replicates and 12 statistics at least as large as the observed statistic, what is the plus-one bootstrap tail estimate?
Why should a failed alternative fit not automatically be counted as a bootstrap non-exceedance?