Information criterion for approximation of unnormalized densities

TL;DR

This paper introduces a new method for approximating unknown densities using a novel information criterion called CIC.

Contribution

The novel contribution is the development of the cross-entropy information criterion (CIC) for model selection in density approximation.

Findings

The CIC is proven to be an asymptotically unbiased estimator of cross-entropy under regularity conditions.

The proposed method effectively selects parametric models that closely approximate target densities in numerical studies.

Abstract

This paper considers the problem of approximating an unknown density when it can be evaluated up to a normalizing constant at a finite number of points. This density approximation problem is ubiquitous in statistics, such as approximating a posterior density for Bayesian inference and estimating an optimal density for importance sampling. We consider a parametric approximation approach and cast it as a model selection problem to find the best model in pre-specified distribution families (e.g., select the best number of Gaussian mixture components and their parameters). This problem cannot be addressed with traditional approaches that maximize the (marginal) likelihood of a model, for example, using the Akaike information criterion (AIC) or Bayesian information criterion (BIC). We instead aim to minimize the cross-entropy that gauges the deviation of a parametric model from the target density. We propose a novel information criterion called the cross-entropy information criterion (CIC) and prove that the CIC is an asymptotically unbiased estimator of the cross-entropy (up to a multiplicative constant) under some regularity conditions. We propose an iterative method to approximate the target density by minimizing the CIC. We demonstrate how the proposed method selects a parametric model that well approximates the target density through multiple numerical studies in the Supporting Information.

Figures

Illustration of our EM algorithm (for k = 3 at t = 1) that updates a randomly initialized density in (a) through Iterations 1–6 in (b)–(g) to approximate the unknown target density in (h).Yellow dots in (a)–(g) are the 2-dimensional pilot data Xi(0),i=1,…, 1000 sampled from the initial distribution Qθ^(0)=Qη. Gray-scale filled contour plots represent the Gaussian mixture density with k = 3 components in (a)–(g) and the target density in (h). White contour line plots in (a)–(g) represent the three component densities of the Gaussian mixture density. Red dashed line is the reference line that marks the shape of the target density in (h). The target density is the optimal importance sampling density in the structural safety example (with b = 1 . 5) in Appendix C.

Plot of the CIC (cumulative version), CIC¯(t)(d), in (17) versus the number of components k, which determines the model dimension d, of the Gaussian mixture model.As the iteration counter t increases from 1 (red solid line) to 4 (blue dash-dot line) to 7 (green dashed line), the CIC is calculated using a larger sample. The CIC is minimized at k = 7 for t = 1, k = 6 for t = 4, and k = 8 for t = 7 in this example. The circles in the plot correspond to the approximate densities shown in Fig 3(a)–3(e).

Gaussian mixture models (GMMs) with k components in (a)–(e) correspond to the circles in Fig 2 and are compared with the target density in (f), which is the optimal importance sampling density in the structural safety example (with b = 1 . 5) in Appendix C.

Equations

Full text

1 Introduction

We consider the problem of finding a parametric approximation of an unknown target density $[eqn]$ from which we cannot sample directly. Instead, we assume the ability to evaluate the unknown non-negative function r satisfying $[eqn]$ where ρ is an unknown normalizing constant. The function evaluation is expensive enough to worth minimizing it. To refer to this problem in shorthand, we coin a term, the Boltzmann Approximation problem (BA problem), owing to its origin in physics. The parametric approximation of $[eqn]$ takes two steps, namely, 1) choosing a parametric family and 2) minimizing the Kullback-Leibler (KL) divergence from a member density in the chosen family to the target density $[eqn]$ . The latter is a well-studied optimization problem and the former is an under-explored model selection problem [1]. A criterion for the model selection is the primary contribution of this paper.

The BA problem is ubiquitous in statistics, including three motivating examples:

Example 1: Simulation-based inference. Simulation models are widely used to estimate a mean $[eqn]$ . We assume the input X follows a known nominal density p ( x ) and run the simulation model to observe a non-negative output, v ( X ) .

When the simulator is computationally expensive, the importance sampling estimator $[eqn]$ is widely used by sampling Xi from q instead of p. This unbiased estimator has theoretically zero variance if Xi, i = 1 , … , n , follows a density $[eqn]$ [3]. This optimal density is the approximation target. By defining r ( x ) = v ( x ) p ( x ) , we encounter the BA problem.

Example 2: Causal inference. Consider a simple causal inference problem where we define a response *Y *(A) under a binary treatment A ∈ { 0 , 1 } . Under the potential outcome model [4],

to investigate the causal effect of the treatment on the response, we need to know the densities p(Y(0)=y0) and p(Y(1)=y1). However, our data only allows us to estimate $[eqn]$ and $[eqn]$ .

Luckily, Tukey’s factorization allows for approximating the counterfactual density

$$[eqn]$$

and similarly for $[eqn]$ . With a model on λ(y0), we aim to sample from the density implied by the factorization. Defining $[eqn]$ and r ( ⋅ ) = λ ( ⋅ ) p ( Y ( 0 ) = ⋅ | A = 0 ) yields the BA problem.

Example 3: Bayesian inference. The most prominent BA problem is approximating a posterior density for Bayesian inference. Let r ( ⋅ ) = λ ( ⋅ ) p ( Y ( 0 ) = ⋅ | A = 0 ) be the prior density of the parameter r ( ⋅ ) = λ ( ⋅ ) p ( Y ( 0 ) = ⋅ | A = 0 ) and $[eqn]$ be the likelihood function. The posterior (density) is $[eqn]$ . In this case, $[eqn]$ is the posterior density and $[eqn]$ is the prior density times the likelihood.

When r is computationally light to evaluate, Markov chain Monte Carlo methods have been used extensively, such as the Metropolis-Hastings algorithm that probabilistically accepts or rejects samples. But there is growing interest in the BA problems where evaluating r is computationally heavy. In Bayesian inference, r may involve computing the likelihood for a large dataset or a complex Bayesian model.

The former case often resorts to the variational inference [5] and the latter case gave rise to the approximate Bayesian computation to bypass evaluating the likelihood function and r [6]. For causal inference, if the model involves covariates, we need to sample counterfactual variables several times per each possible covariate value. For simulation-based inference, running a simulator is often computationally expensive (e.g., minutes or hours). In some of these scenarios where we can only afford a relatively small number (e.g., hundreds or thousands) of evaluations of r , existing rejection-sampling methods are untenable. This paper assumes we want to use every evaluation of r to approximate $[eqn]$ . As a result, the effort (including computational cost) to better use each evaluation is considered worthwhile and likely negligible compared with the cost of evaluating r. Across our simulation studies in the Supporting Information, the average computational time for our algorithm is less than a second per each evaluation of r, while evaluating r can take minutes or hours in the Examples above.

We are interested in both approximating $[eqn]$ and estimating the normalizing constant ρ because ρ is, in some applications, even more important to estimate than the density $[eqn]$ itself. ρ is the model evidence for Bayesian inference and the estimand for importance sampling. The importance sampling context particularly motivates this study and its emphasis on unbiased estimation of ρ (e.g., a rare event probability), which is typically of lesser interest in Bayesian inference.

The existing literature on approximating $[eqn]$ in the BA problem can be generally grouped into parametric and nonparametric approaches [7,8]. This paper considers a class of parametric approaches where we posit a parametric family of densities and find its member closest to $[eqn]$ in terms of the closeness measured by the KL divergence. This parametric framework is very general and includes maximum likelihood estimation. For importance sampling, this framework is used for the so-called cross-entropy method [9].

For this parametric framework, we propose a novel information criterion called the cross-entropy information criterion (CIC). We prove that the criterion is an asymptotically unbiased estimator of the KL divergence (up to a multiplicative constant and an additive constant). We justify that the minimization of the CIC allows for selecting a good family of densities to approximate $[eqn]$ using a limited number of evaluations of r.

We also show that the CIC reduces to the Akaike information criterion (AIC)[10] if we have the ability to sample from $[eqn]$ instead of evaluating r. Rigorous theoretical analysis of the AIC is a long-standing problem in the literature. Theoretical analyses of information criteria akin to the AIC often impose uniform integrability conditions directly on estimators to express the model complexity penalty in terms of the free parameter dimension d . For details, see Conditions A7–A8 in [11], Theorem 1 in [12], a discussion under Eq (7.28) in [13], an approximation from Eqs (2.16) to (2.17) in [14], and more references cited on p. 1157 in[15] and on p. 416 in [17]. For autoregressive models, sufficient conditions for such uniform integrability conditions were established much later than the seminal paper that introduced AIC [10], where the uniform integrability was not established rigorously [17,18]. For general parametric models, establishing general versions of the sufficient conditions remains an open problem. This paper establishes such sufficient conditions for CIC, which is theoretically more general than AIC. Practically, the CIC and AIC are useful for two different classes of model selection problems and hence not interchangeable.

The remainder of this paper is organized as follows. Sect 2 briefly reviews the relevant background. Sect 3 proposes the CIC. Sect 4 explains how the CIC can be used in practice. Proofs, implementation details, and numerical experiments are included in the Supporting Information.

2 Background

This section briefly explains the origin of the BA problem and reviews KL divergence, the maximum likelihood estimator (MLE), and the AIC to a) introduce the minimum cross-entropy estimator (MCE), which is a generalized version of MLE, and b) pave the way for generalizing the AIC to the CIC.

The name of the BA problem originates from statistical physics where the nonnegative function r ( x ) is what we can evaluate and is expressed as $[eqn]$ where ϕ ( x ) is often called the energy of a state x. The target density $[eqn]$ is called the Boltzmann distribution. Of particular significance in physics is evaluating ρ = ∫ r ( x ) d x , called the partition function. The BA problem often has a practically significant consideration: evaluating r is (computationally) expensive, so we want to make use of all evaluations.

KL divergence is commonly used to gauge a difference between two distributions in statistical inference. Consider two probability measures, Q^∗^ and Q, on a common measurable space such that Q^∗^ is absolutely continuous with respect toQ( written Q^∗^≪Q).

The MLE is a prominent example of using the KL divergence. When the data X1,…,X**n are drawn from an unknown Q^∗^ (note that such direct sampling is not possible in the BA problem), we can approximate Q^∗^ by a member in a parametric family $[eqn]$ by minimizing the KL divergence from Qθ to Q^∗^ over $[eqn]$ . Suppose $[eqn]$ for all $[eqn]$ so that the KL divergence is well defined over $[eqn]$ . Also, for a dominating measure μ (e.g., counting or Lebesgue measure), suppose $[eqn]$ for all $[eqn]$ and $[eqn]$ so that densities $[eqn]$ and $[eqn]$ exist. Then, the KL divergence is

$$[eqn]$$

Note that only the second term in (1), called cross-entropy, depends on μ. Therefore, minimizing the KL divergence over $[eqn]$ is equivalent to minimizing the cross– entropy over $[eqn]$ . Because $[eqn]$ is unknown, the cross-entropy should be estimated based on $[eqn]$ . An unbiased, consistent estimator of the cross-entropy is

$$[eqn]$$

The MLE of μ, denoted by $[eqn]$ , is the minimizer of the cross-entropy estimator in (2).

Another example of using the KL divergence or cross-entropy is the Akaike information criterion (AIC) [10]. (We continue to use the same notations as above.) As the dimension d of the free parameter space $[eqn]$ (or equivalently, the model degrees of freedom) increases, the bias of $[eqn]$ approximating Q^∗^ will reduce. To compare the different approximating distributions (or models), we could use a plug-in estimator

$$[eqn]$$

of the cross-entropy. But, this causes an overfitting problem because of the downward bias created from using the data twice (once for $[eqn]$ and another for estimating the cross-entropy).

The AIC remedies this issue by correcting the asymptotic bias of the estimator in (3). The AIC is defined (up to a multiplicative constant) as

$$[eqn]$$

where the bias correction term d ∕ n penalizes the model complexity, balancing it with the goodness-of-fit represented by the first term. Minimizing the AIC is minimizing an asymptotically unbiased estimator of the cross-entropy. Therefore, both the MLE and AIC aim at minimizing the cross-entropy from an approximate distribution Qθ to Q^∗^.

Hereafter, we continue to use the notations Qθ and Q^∗^ for the approximate and target distribution, respectively. But, note that in the BA problem we cannot generate/sample data directly from Q^∗^. Instead, we can evaluate r, which is equal to $[eqn]$ up to an unknown normalizing constant ρ, at a limited number of points.

3 Cross-entropy information criterion for the BA problem

We first introduce an iterative approximation of the cross-entropy in the BA problem, which is suitable for parameter estimation for a fixed-dimension parametric family, in Sects 3.1 and 3.2. In Sect 3.3, we then present the CIC, which can be used for model selection between families with different parameter dimensions. These concepts are combined to yield a sequential algorithm for parameter estimation and model selection. As part of the development, we first define the algorithm in a form which estimates parameters using data from only one stage at a time. In Sect 3.4, we then provide a practical version of the algorithm which uses the cumulative data from all stages.

3.1 Approximate Cross-Entropy (ACE) for the BA problem

Analogous to the MLE and AIC, our approximation task in this paper considers minimizing the KL divergence (or cross-entropy) from a parametric distribution Qθ to the target distribution Q^∗^ over $[eqn]$ (to well-define the KL divergence, hereafter assume $[eqn]$ for all $[eqn]$ ) when the target density $[eqn]$ is proportional to a nonnegative function r, i.e., $[eqn]$ for a positive unknown constant ρ = ∫ r d μ. Recall that in the BA problem, we cannot directly sample from Q^∗^ but can evaluate r at a limited number of points. Minimizing the KL divergence in (1) over r is equivalent to minimizing the cross– entropy - $[eqn]$ and equivalent to minimizing

$$[eqn]$$

which is unknown in practice because r can be evaluated only at observed data points. Using importance sampling, we approximate r in (5) by an unbiased and consistent estimator

$$[eqn]$$

where $[eqn]$ with any parameter $[eqn]$ as long as the support of qη covers the support of $[eqn]$ regardless of any model misspecification (i.e., $[eqn]$ ) in practice. Because the estimator in (6) is to approximate the cross-entropy in the BA problem, we call it the approximate cross-entropy (ACE).

Therefore, by minimizing $[eqn]$ in (6) over $[eqn]$ , we can approximately minimize the KL divergence (or cross-entropy) from Qθ to Q^∗^. Thus, we call

$$[eqn]$$

the minimum cross-entropy estimator (MCE) because it minimizes the ACE. Note that if the random sample is directly drawn from the target distribution (i.e., $[eqn]$ ), then the MCE reduces to the MLE because minimizing (6) is equivalent to minimizing (2) due to $[eqn]$ . More properties of the MCE (e.g., consistency, asymptotic normality, limiting behavior) are characterized in Appendix A.

In the importance sampling literature, the density qη in (6) is called the importance sampling (or proposal) density and $[eqn]$ the optimal importance sampling density. The density $[eqn]$ is optimal because the variance of the following importance sampling estimator of ρ is reduced to zero if $[eqn]$ are drawn from $[eqn]$ :

$$[eqn]$$

In practice, $[eqn]$ is unknown and thus approximated by qη. Therefore, finding qη closest to $[eqn]$ is of primary interest.

3.2 Iterative procedure for minimizing the ACE

A bad choice of sampling density qη may lead to a large variance of the ACE $[eqn]$ in (6). Inspired by the cross-entropy method [9], we propose to update η iteratively by the closest estimate of θ^∗^ as described in Box 1. We use the same sample size n for each iteration for notational simplicity without loss of generality. Later, Sect 3.4 will discuss a more efficient algorithm that uses data in a cumulative fashion.

Box 1: Iterative procedure for approximating Q ^∗^ by minimizing the estimator of cross– entropy from a parametric distribution *Q_θ_

• to Q ^∗^. Iterative procedure for approximatingQ^∗^ Inputs: iteration counter t = 1, the number of iterations τ, the sample size n, and the initial parameter $[eqn]$ .

1. Sample $[eqn]$ from $[eqn]$ .
2. Find the MCE $[eqn]$ , where

$$[eqn]$$

3. If t = τ, output the approximate distribution $[eqn]$ . Otherwise, increment t by 1 and go to Step 1.

An important property of the procedure in Box 1 is that the iterative update of η by the best estimate of θ^∗^ leads to an information criterion that mimics the AIC, as detailed in the next section. If η is fixed, we will not obtain such a nice property.

The iterative approach can be used to estimate ρ, the normalizing constant, by modifying the estimator in (8) as follows:

$$[eqn]$$

Furthermore, if $[eqn]$ , the estimator in (10) is the optimal importance sampling estimator having zero variance. Because the iterative procedure refines $[eqn]$ to be closer to $[eqn]$ , $[eqn]$ will generally have a smaller variance as t gets larger.

Suppose the time complexity for sampling a single observation r ( X ) based on X $[eqn]$ is O ( a ) and the time complexity of optimizing the MCE is O ( b ) . Then the total time complexity for running the algorithm in Box 1 for τ iterations is O ( τ ( na + b ) ) . So when the sample size n is large and/or evaluating r ( ⋅ ) is expensive, most of the cost is from sampling rather than the optimization. So we are allowed to use a more time-consuming optimization program to find the MCE without increasing too much of the total computational cost. In particular, when we use the EM algorithm to find the MCE (see Appendix B for the algorithm details), the time-complexity O ( b ) = O ( K ) , where K is the number of iterations in the EM algorithm. In this case, setting K = O ( n ) will maintain the total time complexity to be O ( τn ) .

As for the space complexity, Step 1 of the algorithm in Box 1 requires O ( n ) to store sampled values. Step 2’s space complexity depends on the optimization algorithm choice. When we use the EM algorithm, its expectation and maximization steps only need to store function values evaluated for n observations at the last EM iteration’s parameters, not the entire EM iteration history. As a result, the required space for the t^th^ iteration isO(dn), whered is the dimension of n. Over τ iterations in Box 1, only the previous iteration’s information needs to be kept, resulting in the total space complexity for Step 2 to remain at O ( dn ) . In practice, because r ( ⋅ ) is expensive to evaluate, storing all its evaluations across τ iterations is sensible, although not required. Considering this addition of O ( τn ) , the total space complexity is O ( dn + τn ) .

3.3 Cross-entropy information criterion (CIC)

To simplify the model complexity penalty term in the AIC, Akaike [10] assumes that the true data-generating distribution belongs to the parametric distribution family being considered. We make a similar assumption $[eqn]$ (i.e., assumption (A1) in Appendix A) to simplify the asymptotic bias of $[eqn]$ in estimating $[eqn]$ .

In what follows we describe our main result as Theorem 1 that quantifies the asymptotic bias. The assumptions and proof are deferred to Appendix A.

Theorem 1 (Asymptotic bias of $[eqn]$ in estimating $[eqn]$ ). Suppose that assumptions (A1-6) and (B1-5) hold. Then

$$[eqn]$$

for each t = 2 , … , τ.

In the above theorem, there are two sources of randomness in the expectation. The first one is $[eqn]$ , and the second is $[eqn]$ . The bias occurs because these two are dependent, and this is the key insight leading to the derivation of an information criterion akin to AIC.

The asymptotic bias, − ρd ∕ n, is proportional to the free parameter dimension d of the parameter space $[eqn]$ , similar to the penalty term of the AIC in (4). In practice, $[eqn]$ is unknown, but we can use a consistent estimator of ρ to estimate the asymptotic bias, such as the estimator in (10).

As a bias-corrected estimator of the cross-entropy $[eqn]$ (up to a multiplicative constant), we define the cross-entropy information criterion (CIC) as follows:

Definition 1 (Cross-entropy information criterion (CIC)).

$$[eqn]$$

for t = 1 , … , τ , where $[eqn]$ with $[eqn]$ in (9). $[eqn]$ is a consistent estimator of ρ, such as the estimator in (10).

We note that the CIC reduces to the AIC up to an additive $[eqn]$ if the samples are all drawn from the target distribution, that is, $[eqn]$ for t = 1 , … , τ in Box 1. If so, the first term of the CIC in (11) becomes

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

because $[eqn]$ in (12) and $[eqn]$ in (13). Plugging the expression in (14) into the CIC in (11) shows that the CIC is equal to ρ times the AIC in (4) up to an additive $[eqn]$ . Note that unless the exact sampling from the target distribution is possible, the CIC remains different from the AIC. Thus, it is generally indefensible to use the AIC in lieu of the CIC for the model selection under consideration in this paper.

The asymptotic bias expression in Theorem 1 holds only for t ≥ 2, because when t = 1, the initial sample is drawn from $[eqn]$ , which is not a distribution converging to Qθ^∗^. Regardless, in practice, one may still use the CIC to select a reasonable parameter dimension d at the first iteration (t = 1).

3.4 The CIC based on cumulative data

If we use the equal sample size n for each iteration, the model dimension d for later iterations may vary only a little from the earlier iterations. Alternatively, we can aggregate the samples gathered through iterations to obtain a cumulative version of the CIC as discussed in this subsection.

In the t^th^ iteration, the cumulative version uses all the observed data up to the current iteration to estimate d, instead of using only the current iteration’s data $[eqn]$ (recall Box 1). The benefit of the cumulative version is the tendency of the aggregated estimator of d to have a smaller variance than the non– aggregated estimator $[eqn]$ in (9). This approach, in turn, can reduce the variance of the MCE as well, which minimizes the aggregated estimator of d.

For more flexibility, we can allocate a different sample size for each iteration, that is, n_t_ for the t^th^ iteration, t = 0 , 1 , … , τ (for example, a large n0 for the initial sample to broadly cover the support of Qη and equal sample sizes $[eqn]$ for the later iterations). Then, we can find the MCE

$$[eqn]$$

where the aggregated estimator of t = 0 , 1 , … , τ is denoted as

$$[eqn]$$

for t = 1 , … , τ with $[eqn]$ . Note that $[eqn]$ in (16) is an unbiased estimator of t = 1 , … , τ.

By using all data gathered up to the t^th^ iteration, we can determine the model parameter dimension d at the t^th^ iteration with the following CIC:

Definition 2 (CIC: Cumulative version).

$$[eqn]$$

for t = 1 , … , τ , where $[eqn]$ in (15). $[eqn]$ is a consistent estimator of ρ, such as the estimators in (18) and (19).

As t increases, the accumulated sample size $[eqn]$ increases so that the free parameter dimension d can increase. Thus, the cumulative version of the CIC allows the use of a highly complex model if it can better approximate Q^∗^.

Note that the cumulative version of the CIC in (17) is motivated by the non-cumulative version in (11), but it remains unclear if a result similar to Theorem 1 will appear in this case. The major challenge is that the latest iteration depends on information from all previous iterations, so the analysis on the bias becomes more complicated.

As a consistent and unbiased estimator of ρ (akin to the estimators in (6) and (8)), we can use

$$[eqn]$$

at the 1*^st^* iteration. At the t^th^ iteration for t = 2 , … , τ, we can use

$$[eqn]$$

where we do not use the data $[eqn]$ from the initial distribution $[eqn]$ because they could potentially increase the variance of the resulting estimator if Qη is too different from Q^∗^. The estimator in (19) is an importance sampling estimator of $[eqn]$ . The potential for the increased variance has been well studied in the importance sampling literature, including defensive techniques [19,20]. Owen and Zhou’s method (SEIS) is implemented and tested in S2 Additional experiments (Interested readers can find cemSEIS.py in our Python package linked in S3 Code). Note also that one can choose the initial distribution (denoted by $[eqn]$ thus far) to be any distribution whose support covers the support of Q^∗^ to keep the estimator in (18) unbiased, although a judicious choice can reduce the estimator’s variance. Future research may investigate how to assign greater weights on newer observations to improve the estimation depending on the asymptotic behavior of $[eqn]$ approaching $[eqn]$ .

4 Application of the cross-entropy information criterion

This section details how the cumulative version of CIC can be useful in practice. Henceforth, CIC refers to the cumulative version (Definition 2), not the non-cumulative version (Definition 1), unless specified otherwise. We first present how the CIC can help choose the number of components, k, for a mixture model in conjunction with an expectation-maximization (EM) algorithm that finds the MCE for a given k. Then, we present the summary of how to use the CIC to iteratively approximate a target distribution. In Appendix C, we provide numerical examples to illustrate the use of the CIC for approximating the optimal importance sampling density and the posterior density in Bayesian inference. Additional numerical experiments are included in S2 Additional experiments, which investigate when the CIC-based importance sampling using a Gaussian mixture model works well or not. Interested readers are also referred to the work of [21], which applies the CIC-based importance sampling to stochastic simulation models, where the evaluation of r is stochastic (i.e., r ( x ) is not a function, but it follows a distribution that depends on x).

4.1 Simulation: Mixture model and an EM algorithm

To approximate a target distribution, we consider a parametric mixture model with a parameter dimension d. Parametric mixture models are often used to approximate a posterior density for Bayesian inference [5] and an optimal importance sampling density [3,22,23]. The density approximation quality hinges on the number of mixture components, k (or equivalently, the parameter dimension d). Prior studies either assume that k is given [3,22] or use a rule of thumb to choose k based on “some understanding of the structure of the problem at hand” [23].

We can use the CIC to select k for any parametric mixture model, considering various parametric component families. For example, exponential families are especially convenient because the MCE can be found by using an expectation-maximization (EM) algorithm. In this paper, we use the Gaussian mixture model (GMM) for illustration. Appendix B details our version of the EM algorithm to find the MCE. Future research may investigate computationally more efficient methods based on the recent advancement of distributed computing and gradient-based methods [2].

Fig 1 illustrates our EM algorithm in action for the first outer iteration (t = 1), where a GMM (gray-scale filled countour plot) with three component densities (white countour lines) is updated over EM iterations to approximate an unknown target density. We can see that in contrast to the conventional EM algorithm that maximizes the likelihood of a model (i.e., goodness-of-fit) to approximate the distribution of observed data, our algorithm uses the data (yellow dots), $[eqn]$ , to estimate and minimize the cross-entropy from the approximate density to the target density. Out of the 1000 observations (yellow dots) in Fig 1(a) (note that the same data are plotted in (b)–(g) as well), only the small portion of them that fall above the red dashed line contribute to the cross-entropy estimate (in (16) because r ( x ) is zero below the red dashed line for the structural safety example in Appendix C). See Appendix D for more details of the EM algorithm (including its parameter initialization strategy based on [16]) implemented for the structural safety example in Appendix C.

Illustration of our EM algorithm (for k = 3 at t = 1) that updates a randomly initialized density in (a) through Iterations 1–6 in (b)–(g) to approximate the unknown target density in (h).Yellow dots in (a)–(g) are the 2-dimensional pilot data Xi(0),i=1,…, 1000 sampled from the initial distribution Qθ^(0)=Qη. Gray-scale filled contour plots represent the Gaussian mixture density with k = 3 components in (a)–(g) and the target density in (h). White contour line plots in (a)–(g) represent the three component densities of the Gaussian mixture density. Red dashed line is the reference line that marks the shape of the target density in (h). The target density is the optimal importance sampling density in the structural safety example (with b = 1 . 5) in Appendix C.

4.2 Summary of the CIC-based distribution approximation procedure

This subsection summarizes how we can use the CIC to approximate a target distribution in practice. Using the EM algorithm in Sect 4.1 for different k’s (or d’s) in the t^th^ iteration for t = 1 , … , τ, we can find the MCE in (15) and calculate the CIC in (17). At the minimum of the CIC, we can then find the best number of components, $[eqn]$ ( or the best model dimension $[eqn]$ to use in the t^th^ iteration.

The CIC tends to decrease and then slowly increase as d increases, subject to the randomness of the data. Fig 2 shows such a pattern, where k is the number of mixture components in the GMM with unconstrained means and covariances. Note that k is proportional to the free parameter dimension d = ( k − 1 ) + k ( p + p ( p + 1 ) ∕ 2 ) , with p denoting the dimension of the GMM density support.

Within the t^th^ iteration, we calculate the CIC over k as shown in Fig 2 for t = 1, 4, 7 for the structural safety example in Appendix C. As the iteration counter t increases, $[eqn]$ in (17) uses a larger sample that accumulated data over iterations to more accurately estimate the cross-entropy.

Plot of the CIC (cumulative version), CIC¯(t)(d), in (17) versus the number of components k, which determines the model dimension d, of the Gaussian mixture model.As the iteration counter t increases from 1 (red solid line) to 4 (blue dash-dot line) to 7 (green dashed line), the CIC is calculated using a larger sample. The CIC is minimized at k = 7 for t = 1, k = 6 for t = 4, and k = 8 for t = 7 in this example. The circles in the plot correspond to the approximate densities shown in Fig 3(a)–3(e).

At t = 1, the GMM with k = 7 in Fig 3(b) achieves the minimum CIC (as shown in Fig 2), while k = 1 and k = 10 result in seemingly over-simplified and over-complicated densities in Figs 3(a) and 3(c), respectively, for the given sample size, 1000 (note that the effective sample size is much smaller because only a small portion of the data fall above the red dashed line as explained earlier with Fig 1). A good choice of k (neither too small nor too large) for the given sample size at the current iteration helps subsequent iterations by preventing sampling from an overly simplified/complicated distribution that could misguide the later iterations. Thus, it is beneficial to refine the approximate distribution proportionally (neither too much nor too little) for the given data size as guided by the CIC. Over iterations, sampled data (yellow dots in Fig 1) should increasingly cover the entire support of the target density. The CIC-minimizing densities at t = 4 in Fig 3(d) and t = 7 in Fig 3(e) capture the overall shape of the target density in Fig 3(f).

Box 2 summarizes the CIC-based distribution approximation procedure. Note that in addition to approximating the target distribution Q^∗^, if we want to estimate a quantity of interest such as ρ, we can sample $[eqn]$ and use the estimator such as $[eqn]$ in (19). The numerical examples in the Supporting Information use this additional step.

Gaussian mixture models (GMMs) with k components in (a)–(e) correspond to the circles in Fig 2 and are compared with the target density in (f), which is the optimal importance sampling density in the structural safety example (with b = 1 . 5) in Appendix C.

Box 2: CIC-based approximation of the target distributionQ ^∗^. CIC-Based Approximation of the Target Distribution Q^∗^ Inputs: iteration counter t = 1, the number of iterations τ, the sample size per iteration $[eqn]$ , the initial parameter dimension d^(0)^, and the initial parameter $[eqn]$ .

1. Sample $[eqn]$ .
2. Find the best model dimension $[eqn]$ to use, where $[eqn]$ is in (17) with the MCE $[eqn]$ in (15). $[eqn]$ is a consistent estimator of ρ, such as the estimators in (18) and (19).
3. If t = τ, output the approximate distribution $[eqn]$ . Otherwise, increment t by 1 and go to Step 1.

5 Conclusion

This paper proposed the cross-entropy information criterion (CIC) to find a parametric density that has the asymptotically minimum cross-entropy to a target density to approximate. The CIC is the sum of two terms: an estimator of the cross-entropy (up to a multiplicative constant) from the parametric density to the target density, and a model complexity penalty term. Under certain regularity conditions, we proved that the penalty term corrects the asymptotic bias of the first term in estimating the true cross-entropy. Empirically, we demonstrated that minimizing the CIC leads to a density that well approximates a target density.

The CIC allowed us to develop a principled algorithm to near-automatically approximate an unknown density that can be evaluated up to a normalizing constant at a limited number of points. The necessary manual “tuning” in practice is minimal as it pertains primarily to the selection of the initial sampling distribution, $[eqn]$ . It can be chosen judiciously based on a priori knowledge about where r is expected to be large (or at least, non-zero) once one determines the mixture component distribution (e.g., Gaussian vs. something else) and the algorithm for minimizing the CIC (e.g., EM algorithm). This paper’s CIC-based algorithm is made publicly available as the first off-the-shelf software package for solving the BA problem.

Supporting information

S1 AppendicesIt includes Appendix A (Assumptions and Proofs), Appendix B (EM Algorithm for Minimizing the Cross-Entropy), Appendix C (Numerical Experiments: Importance Sampling and Bayesian Inference), and Appendix D (Implementation Details of Numerical Experiments).(PDF)

S2 Additional experimentsIt includes additional numerical experiments that investigate the empirical performance of the CIC-based importance sampling through two numerical examples and one case study(PDF)S3 Code. The code for reproducing all experimental results in the paper is publicly available as a Python package on https://pypi.org/project/cicriterion/. Its archived version’s DOI is doi: 10.5281/zenodo.13901261