{'fraction_non_alphanumeric': 0.035517131682853344, 'fraction_numerical': 0.008533093807893112, 'mean_word_length': 4.7741761210156675, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A layered multiple importance sampling scheme for focused optimal Bayesian experimental design

We develop a new computational approach for"focused"optimal Bayesian experimental design with nonlinear models, with the goal of maximizing expected information gain in targeted subsets of model parameters. Our approach considers uncertainty in the full set of model parameters, but employs a design objective that can exploit learning trade-offs among different parameter subsets. We introduce a new layered multiple importance sampling scheme that provides consistent estimates of expected information gain in this focused setting. This sampling scheme yields significant reductions in estimator bias and variance for a given computational effort, making optimal design more tractable for a wide range of computationally intensive problems.

1. Introduction. Data acquired from experiments are essential for developing and refining models. Acquiring data can be a costly process, so it is important to carefully choose experimental conditions that maximize the value of the acquired data, as defined by some optimal experimental design criterion [3]. In general, the design criterion should depend on what one intends to do with the data acquired from an experiment.
In this paper, we consider the optimal design of experiments where the observed data are used to infer uncertain model parameters. This setting includes models that may contain a large number of parameters, but where the user may care only about inference of a targeted subset. For example, in the context of chemical kinetics, experimental outcomes might be influenced by multiple rate constants and thermodynamic parameters, but only the rates of certain subsets of reactions may be of interest. Another example is the Kennedy & O'Hagan model calibration framework [21], where users may favor experimental designs that focus on reducing uncertainties in 'physical' parameters over reducing uncertainty in 'tuning' parameters. To this end, we employ a design criterion that favors experiments maximizing the expected information gain in user-selected parameters of interest; following [24], we refer to this approach as focused optimal experimental design. Our key contribution is a new multiple importance sampling framework that makes focused experimental design tractable.
We formulate our experimental design problems in a Bayesian setting, where uncertain parameters are endowed with probability distributions and parameter inference can naturally be performed when the data are noisy and limited in number [5,37]. We adopt a decisiontheoretic approach, following [27], wherein an optimal design maximizes the expected utility of an experiment. As suggested by [26], our utility function is chosen to be the information gain due to an experiment, which is equivalent to the Kullback-Leibler divergence from the prior to the posterior-or from selected prior marginals to the corresponding posterior marginals. This nonlinear Bayesian optimal design framework is a generalization of Bayesian D-optimal design, making no assumptions on the linearity of the forward model or the Gaussianity of the prior and noise distributions [9].
For models where the observables depend nonlinearly on the parameters of interest, exact evaluation of this design criterion is in general intractable, as it requires integration over the possible outcomes of the experiment and the space of model parameters. A nested Monte Carlo estimator [34] of the expected information gain (jointly over all model parameters) is proposed in [36], but has significant limitations due to computational expense. Alternatives to this nested Monte Carlo estimator include kernel estimators of mutual information as proposed in [23] and applied in [38], but these approaches involve density estimation and thus require careful parameter choices to control bias and variance, along with large sample populations. Laplace approximations to the posterior distribution [29,7] or computationally efficient surrogate models [19] can result in more tractable design criteria, but are unsuitable when these approximations break down, e.g., as the posterior distribution departs from normality or when a sufficiently accurate surrogate cannot be constructed. In general, these approximations induce asymptotic bias and do not yield consistent estimators of the expected information gain. Of course, it is possible to use Laplace approximations to construct consistent importance sampling estimators, as recently proposed in [4] and discussed below; these techniques are complementary to the Monte Carlo methods developed here. Further reviews of modern computational methods and approximations for decision-theoretic Bayesian design can be found in [35,41].
Other recent experimental design formulations avoid computing the expected information gain criterion altogether, and instead begin with the Bayesian A-optimality criterion [1] for linear problems, extending it to nonlinear settings by (i) using Laplace approximations of the posterior for tractability; and (ii) averaging over many possible realizations of the data [2]. These approaches are similar to nonlinear design criteria that average functionals of the expected or observed Fisher information matrix, as described in [9]. Though these alternatives provide computational scalability and can be rationalized in weakly nonlinear problems, they remain only approximations of a decision-theoretic expected utility. Also, while D-optimal design and its generalization to expected information gain are invariant under reparameterization of the underlying model, A-optimal design and its extensions discussed above are not [6,12]. This paper proposes a consistent, yet tractable Monte Carlo estimator of the expected information gain, not limited by linear approximations. Estimating the expected information gain requires estimating posterior normalizing constants for many realizations of the data. Yet estimating the posterior normalizing constant-even for a single realization of the data-is a challenging task. In addition, our criterion for focused experimental design introduces marginalization over nuisance parameters [25], further reducing the efficiency and accuracy of existing Monte Carlo approaches. While many schemes [10,17] have been developed to estimate posterior normalizing constants, they only target a single posterior distribution. Estimating posterior normalizing constants over many realizations of the data, however, allows for the sharing of information between estimators. We use this idea to develop a layered multiple importance sampling (LMIS) estimator of the expected information gain, which iteratively improves the biasing distributions used to estimate posterior normalizing constants. This is accomplished by optimally combining samples and cached model evaluations from earlier iterations using multiple importance sampling [39], yielding consistent and more accurate estimators.
Our approach has roots in the simpler importance sampling scheme of [15]. A useful generalization of the latter, recently proposed by [4,14], involves computing a Laplace approximation of the posterior for each realization of the data, and using this approximation to construct a biasing distribution to estimate each posterior normalizing constant. The resulting expected information gain estimators are also consistent, but require solving an optimization problem and computing a Hessian for each realization. In contrast, the LMIS scheme developed here emphasizes recycling information from one realization of the data to the next. Within this framework, many different biasing schemes can be used. The particular biasing schemes employed here avoid repeated optimization and do not require gradients or Hessians of the forward model or log-likelihood. LMIS is thus entirely non-intrusive, in the sense of requiring only pointwise evaluations of the likelihood function.
We demonstrate using both linear and nonlinear numerical examples examples that our LMIS estimator offers multiple order-of-magnitude reductions in mean squared error over comparable Monte Carlo methods using the same number of model evaluations, while remaining asymptotically exact. This new estimator enables efficient estimation of expected information gain for both joint and the more challenging case of focused experimental design. This paper is organized as follows. In section 2, we introduce the Bayesian optimal experimental design problem and formulate the focused information gain criterion. We illustrate this formulation using a linear-Gaussian toy problem-which has closed-form expressions of the expected information gain, and a clear trade-off between inferring all the model parameters and inferring a subset thereof. In section 3 we derive a nested Monte Carlo estimator of our focused design criterion, allowing for a generic choice of importance sampling scheme in its inner loops; we then analyze the asymptotic bias and variance of this estimator. Subsection 3.4 describes our layered multiple importance sampling scheme. Section 4 evaluates the numerical performance of the LMIS estimator on linear problems of varying dimension, and in subsection 4.2 we apply the LMIS algorithm to the optimal design of a Mössbauer spectroscopy experiment, which is nonlinear in its parameters.
2. Optimal experimental design. The Bayesian paradigm treats unknown parameters as random variables. Let (Ω, F, P) be a probability space, where Ω is a sample space, F is a σ-field, and P is a probability measure on (Ω, F). Let the real-valued random vector z : Ω → R n θ +n η represent the uncertain parameters of a statistical model, and let it be partitioned coordinate-wise, z := (θ, η), into θ : Ω → R n θ and η : Ω → R n η , which denote the parameters of interest and nuisance parameters, respectively. Both θ and η are parameters to be conditioned on experimental data. The random variable z is associated with a prior measure µ on R n θ +n η such that µ(A) = P(z −1 (A)) for A ∈ R n θ +n η . We can then define p(z) = dµ/dz = p(θ, η) as the joint density of θ and η with respect to the Lebesgue measure dz. For the present purposes, we will assume that this density always exists. Similarly, the data y can be treated a real-valued random variable y : Ω → R n y . Finally, let d ∈ D ⊆ R n d denote the design variables, or experimental conditions. Hence, n θ is the number of uncertain parameters of interest, n η is the number of nuisance parameters, and n d is the number of design variables. In the Bayesian setting, upon observing a realization of the data y obtained by performing an experiment under conditions d, we can update our knowledge about the model parameters using Bayes' rule: Here p(θ, η) is the prior density, p(y|θ, η, d) is the likelihood function associated with our statistical model, p(θ, η|y, d) is the posterior density, and p(y|d) is the evidence or posterior normalizing constant. It is reasonable to assume that prior knowledge on θ and η (i.e., the prior measure µ) does not depend on the experimental design, and thus we employ the simplification p(θ, η|d) = p(θ, η).
The likelihood function often combines a deterministic forward model with a statistical model for measurement and/or model errors. For example, the observations might be modeled as additive Gaussian perturbations of some deterministic model predictions G(θ, η, d): Here, G is the forward model; we will assume that it is (in general) nonlinear and that its evaluations are the dominant computational cost of evaluating or simulating from the likelihood. Thus, while evaluations of G(θ, η, d) for new values of the parameters or design variables are expensive, it is cheap to evaluate p(y|θ, η, d) for different values of y, given fixed values of θ, η, and d. This assumption is representative of many science and engineering problems, where G might represent the solution of a set of differential equations and where the noise model is relatively simpler.
2.1. Expected utility framework. Following the decision-theoretic approach of Lindley [27], an objective for optimal experimental design should have the following general form: where u(y, d, θ, η) is a utility function representing the value of observing a particular outcome y of the experiment under design d, when the parameters are θ and η. We do not know the outcome of the experiment a priori, however; nor do we know the true values of the parameters θ and η. So we instead consider the expected utility U (d), where the expectation is taken over the joint distribution of θ, η, and y, and where Θ ⊆ R n θ , H ⊆ R n η and Y ⊆ R n y are the ranges of θ, η, and y, respectively.
To focus attention on inference of the parameters θ, we choose u(y, d, θ, η) to be the Kullback-Leibler (KL) divergence (i.e., relative entropy) from the marginal prior on θ to the marginal posterior of θ, where the marginalization is over the nuisance parameters η. In other words, the divergence is from the marginal prior on the parameters of interest, to the marginal posterior on the parameters of interest, and hence our utility function is Above,θ is simply a dummy variable of integration. The final equality reflects the fact that the KL divergence does not actually depend on the values of θ and η. Thus we can write the expected utility (2.1) of a design as an expectation over the data only: Now, substituting the utility (2.2) into the objective function (2.3), we obtain which means that the expected utility U (d) is the expected information gain in the parameters of interest θ, irrespective of the information gain in the nuisance parameters η. Intuitively, the expected information gain describes the average reduction in entropy from the prior marginal of θ to the posterior marginal of θ, where the average is taken over all possible experimental outcomes y, weighed according to their prior predictive distribution. U (d) is equivalent to the mutual information between the parameters of interest θ and the data y, given a particular design d [11].

Linear-Gaussian example.
To illustrate our formulation for focused optimal experimental design, we first apply it to a toy problem that can be solved analytically. Consider the following experiment with a linear model, additive independent noise 1 , 2 ∼ N (0, σ 2 ), and a design space d ∈ [0, 1]. The prior distributions on θ and η are chosen to be independent standard normals, θ ∼ N (0, 1) and η ∼ N (0, 1). The posterior θ, η|y, d ∼ N (µ post , Σ post ) is a multivariate normal with mean and covariance available in closed form. The marginal posterior p(θ|y, d) is also Gaussian and follows immediately.
The structure of this model encodes a simple trade-off between learning about θ and learning about η. As d increases, we obtain a larger signal-to-noise ratio in the measurement of θ while the signal-to-noise ratio associated with measuring η decreases. This trade-off can be visualized by observing the anisotropy of the posterior distribution when varying d for a fixed value of y, as illustrated in Figure 1. The marginals at the bottom of the figure suggest that the design minimizing uncertainty in θ is simply d = 1. Since we are working with a linear-Gaussian model, we can obtain analytical expressions for the expected information gain in (θ, η) jointly, and in the parameter of interest θ alone, The expected information gains for both cases are plotted on the right of Figure 1.
To find the optimal experimental design, the expected utility must be maximized over the design space D = [0, 1]. When the objective is to maximize information gain in the parameters of interest, we have U = U θ and the optimal design is d * θ = 1; this is consistent with the θ marginals shown in Figure 1. On the other hand, when the objective is to maximize information gain in all uncertain parameters, we have U = U θ,η and the optimal design becomes d * θ,η = 1 2 . Notice that maximizing the expected information gain in (2.5) is equivalent to minimizing the determinant of the posterior covariance matrix, i.e., Bayesian D-optimal design [9]. We will revisit more complex versions of this example in subsection 4.1, where we benchmark the convergence of various Monte Carlo estimators of (2.4).

Monte
Carlo estimation of the expected information gain. In general, the expected information gain (EIG)-whether in the parameters θ alone or jointly in (θ, η)-does not have a closed form and must be approximated numerically, e.g., by a Monte Carlo estimate. In this section, we derive a nested Monte Carlo estimator of the focused EIG (2.4), analyze its properties, and then propose a layered multiple importance sampling (LMIS) estimator that has better properties. To make the subsequent developments clear, we note that the nested estimators employed here comprise an outer loop that involves taking an expectation over the prior distribution and an inner loop that involves computing certain posterior normalizing constants. Subsection 3.1 presents a general format for nested estimators of (2.4) that employ importance sampling in their inner loops, but deliberately leaves the choice of biasing distributions open. Subsection 3.2 analyzes the bias and variance of these estimators; this analysis also suggests what good biasing distributions ought to be. Subsection 3.3 discusses the convergence rate of the nested estimators and the optimal asymptotic scaling of inner loop and outer loop sample sizes. Subsection 3.4 then presents a way of adaptively identifying biasing distributions via multiple importance sampling, which leads to our new LMIS estimator. LMIS can offer significant increases in performance over simpler and non-adaptive choices of biasing, while preserving asympotically exact estimation of the focused EIG.
3.1. A nested importance sampling estimator. The EIG as written in (2.4) contains an inner loop of integration over the posterior distribution, and also requires evaluating the normalized posterior density of θ; both are computationally undesirable. By applying Bayes' rule to the terms involving the posterior density, we can rewrite (2.4) so that integration is instead performed over the prior distribution of θ, η, and y: This form of the focused EIG can be approximated by the following Monte Carlo estimate: from the prior as follows: for each i, draw (θ (i) , η (i) ) ∼ p(θ, η) from the joint prior on the parameters; then draw y (i) ∼ p(y |θ (i) , η (i) , d).
The sample η (i) is necessary to realize a value of y (i) but does not appear explicitly in (3.1).
The summand in (3.1) contains the posterior normalizing constant p(y (i) |d) for a given realization of the data y (i) and design parameters d, known as the marginal likelihood, along with a partial posterior normalizing constant p(y (i) |θ (i) , d), which we call the conditional likelihood. The marginal likelihood typically does not have an analytical form, but can be approximated by an importance sampling estimate for each y (i) : with corresponding importance weights: .

(3.6)
Substituting the above importance sampling estimates of the marginal and conditional likelihoods into (3.1), we arrive at a nested Monte Carlo estimator for the EIG in θ: So far, we have not explicitly specified the biasing distributions In the next sections, we will identify optimal biasing distributions and describe methods to approximate them.

Bias and variance of the nested estimator.
Choosing an appropriate biasing distribution requires understanding the bias and variance of the focused EIG estimator (3.7). Using the second-order delta method (see Appendix A), we find that the bias of U (d) is, to leading order: and that its variance is, to leading order: The higher-order terms in the bias (3. ); see Appendix A for details. Also, σ 2 marg (y, d) and σ 2 cond (y, θ, d) are variances of the summands in the marginal and conditional likelihood estimators (3.2) and (3.5), respectively: From (3.8), we see that the bias is controlled by the variance of the inner marginal and conditional likelihood estimators, which can be reduced by choosing better biasing distributions q (y) marg and q (y,θ) cond or by increasing the number of inner Monte Carlo samples M 1 and M 2 . 1 The bias can take both positive and negative values since it is the difference of two positive quantities. This property contrasts with the nested Monte Carlo estimator of unfocused EIG in [36], where the bias only takes positive values. The latter corresponds to the case of no nuisance parameters (n η = 0); then the conditional likelihood p(y|θ, d) (3.4) can be evaluated exactly, so σ 2 cond (y, θ, d) is effectively zero and (3.8) is always positive. The variance of U (d) decreases with the number of outer samples as 1/N . The first term, D 3 (d), depends only on the variability of the conditional and marginal likelihoods given (y, θ) drawn from the prior. The second and third terms, D 1 (d) and D 2 (d), depend explicitly on the variance of the inner estimators. Thus, reducing the variance of the inner estimators will also reduce the variance of U (d). 1 Here and in (3.10a)-(3.10b), we write the inner biasing distributions as q (y) marg and q (y,θ) cond to indicate that they can depend arbitrarily on the outer values of y and possibly θ. To recover the simpler notation of subsection 3.1, where i indexes the outer loop, simply put q cond and q In order to minimize the bias and variance of U (d), we therefore need to find biasing distributions q (y) marg and q (y,θ) cond that minimize the variance of the importance sampling estimates of the marginal and conditional likelihood. From (3.10a) and (3.10b), we see that the using the joint and conditional posteriors as biasing distributions (with normalized densities) would result in zero-variance estimators. Current approaches for evaluating the marginal likelihood in this setting, such as [19,36], instead use the prior distribution p(θ, η) as a biasing distribution, i.e., so that the importance weights w marg and w cond in (3.3) and (3.6) are identically one. However, when the joint or conditional posterior distributions are more concentrated than their corresponding prior distributions, e.g., when the data y (i) are informative, using the prior distribution as a biasing distribution can result in very high-variance estimates of the marginal and conditional likelihoods, which ultimately lead to large bias (and variance) when estimating the expected information gain. (We will demonstrate this issue in section 4.) To avoid this problem, a more sophisticated approach is needed to find biasing distributions that more closely resemble the joint and conditional posteriors.
3.3. Sample sizes and scaling. The nested estimator (3.7) is consistent (see [34]) and in particular, its bias and variance converge to zero as M 1 , M 2 , N → ∞. Since all three sample sizes must increase, it is natural to ask whether there is an optimal allocation of inner (M 1 , M 2 ) and outer (N ) samples for a given computational budget For simplicity, let M 1 = M 2 = M so that W = 2M N and define C := (C 1 − C 2 ) 2 > 0, where we drop dependence on d to reduce clutter. Let α 2 = 2M/N be the ratio between the number of inner and outer Monte Carlo samples. The value of α and the corresponding ratio α 2 that minimize the leading-order MSE are This result suggests that the ratio between the number of inner and outer samples should decrease slowly as the computational budget increases, in particular scaling as W −1/3 . The ratio C/D 3 reflects the relative magnitudes of leading-order terms in the bias and variance; if the bias dominates, one should choose a relatively larger number of inner loop samples. A detailed derivation of the optimal scaling and a numerical study of the asymptotic behavior are presented in Appendix B.
It is important to emphasize, however, that all these results are asymptotic in M, N → ∞. For finite M and N , higher-order terms in the bias and variance may come into play. Also, the constants in (3.11) are in practice difficult to estimate, and an incorrect value-one that does not capture the relative magnitudes of the bias and the variance-may cause the α ∼ W −1/6 scaling to perform worse than a fixed M/N ratio for finite sample sizes. An ideal approach might involve directly sampling from the posterior corresponding to each y (i) , but doing so is computationally taxing-and, moreover, should only be done when evaluating the corresponding normalized posterior density (to avoid instabilities of the harmonic mean estimator [31]), which requires knowledge of the quantity we actually wish to estimate. Instead, we will restrict our choice of biasing distributions to a parametric family with known density and use approximations of the joint and conditional posteriors, within this family, as biasing distributions. We describe our approach below using multivariate tdistributions; since these distributions are heavy-tailed, they are broadly useful in importance sampling [33]. As we shall emphasize later, however, our approach can easily be used with other parametric families.
To identify a multivariate t-distribution that is good for biasing relative to a given posterior, it is useful to estimate the posterior mean and covariance. Ideally, we would like to estimate these moments without drawing additional samples that require costly forward model evaluations. Our key idea is as follows: if we proceed sequentially through the outer loop of (3.7), indexed by i, then we can use the samples drawn from each inner-loop biasing distribution in the previous steps {1, . . . , i − 1} to help characterize the posterior corresponding to the current y (i) . To reuse these existing samples, we consider a multiple importance sampling approach [39], which combines the previously-drawn samples into a single mixture distribution that is used to estimate the mean µ To this end, we first introduce self-normalized importance sampling estimators of the posterior mean and covariance: are drawn from a biasing distribution q (i) post (θ, η) and the unnormalized importance weights are given by , which can be rescaled to obtain the self-normalized importance weights post be a mixture distribution whose components are distributions from which we have already drawn samples-e.g., the prior distribution (since we can easily draw all N outer loop samples in advance) and the inner-loop biasing distributions from earlier iterations: Instead of drawing new samples from this mixture distribution, we simply realize that the existing samples are a draw from q (i) post : Moreover, because the forward model G has already been evaluated each sample of (θ, η) in this empirical mixture, evaluating the likelihood at the current data point y (i) will be very cheap. This process results in the incremental enrichment of the samples used to estimate the posterior mean µ Since it is helpful for biasing distributions to have heavier tails, we typically choose small ν, e.g., ν = 2.5. One possible choice of q (i) cond is another multivariate t-distribution derived from conditionals of the posterior Gaussian approximation: where we used the block-matrix notation Alternatively, one could use the conditional of q (i) marg directly; note that this is just another multivariate t distribution, very similar to (3.16), but with a scaling of the scale matrix Σ η|θ := Σ η − Σ T θη Σ −1 θ Σ θη and a modification of the degrees of freedom ν [13]. Another possible choice is simply the η-marginal of q (i) marg . This last choice is more diffuse, since it does not depend on θ (i) , but may be useful in situations where the preceding approximations to the true posterior conditional are poor.
Note, finally, that the importance sampling estimates for the marginal and conditional likelihood involve drawing new samples {(θ (i,j) , η (i,j) )} We call this complete strategy layered multiple importance sampling, with two layers: 1. A biased layer (with bias due to self-normalization) to estimate posterior moments, via multiple importance sampling with sample reuse; 2. An unbiased layer to estimate the marginal and conditional likelihoods, using biasing distributions selected from a parametric family according to the estimated posterior moments.
3.4.2. Pruning the mixture biasing distribution. Although layered multiple importance sampling does not require additional model evaluations beyond a standard nested scheme, if we include all the biasing distributions from earlier iterations in (3.13), the number of likelihood evaluations-using cached model evaluations-grows as O(N 2 M 2 1 ), which can result in a nontrivial computational overhead. Instead, we can include only the most "useful" samples in the empirical mixture distribution. Since we do not know the posterior beforehand, we can use the (θ (i) , η (i) ) which gave rise to the data y (i) as our best guess for where the posterior p(θ, η|y (i) , d) is most concentrated. This selection criterion allows us to define an index set where we only include distributions that have a higher probability density on the sample (θ (i) , η (i) ) than the prior. Then, the mixture density (3.13) can be modified as 18) and the samples drawn from this mixture distribution q (i) post can be written as a disjoint union of samples drawn previously from the prior and biasing distributions: where the total number of samples is The specification of an index set J (i) allows some generality: instead of the strategy in Our overall strategy-layered multiple importance sampling with mixture pruning-is illustrated via the flowchart in Figure 2. 3.4.3. Ordering samples. Since the layered multiple importance sampling approach is sequential by construction, a natural question is whether the ordering of the prior samples {(θ (i) , η (i) , y (i) )} N i=1 will affect the behavior of the overall estimator. Since we draw these samples before evaluating the inner estimators, we have the freedom to specify a particular ordering. Intuitively, a good ordering is one that maximizes incremental collection of existing samples  enrichment-such that samples drawn from earlier biasing distributions will be reused in many subsequent mixture distributions. This goal suggests an ordering that prioritizes the posterior distributions p(θ, η|y (i) , d) that overlap the most with all other posterior distributions. However, since the N different posterior distributions are not characterized before the ordering is chosen, we suggest using the prior density of p(θ (i) , η (i) ) as a surrogate.
With sufficient data and under appropriate consistency conditions, the posterior distributions p(θ, η|y (i) , d) will be concentrated near the θ (i) , η (i) that gave rise to the y (i) ; more generally, the centers of the posteriors will (collectively) be distributed similarly to the prior. Thus, ordering the samples by decreasing prior density should achieve a similar effect as ordering the posteriors by decreasing overlap. Our complete algorithm, LMIS with mixture pruning and an initial sorting of samples, is described in Algorithm 1.

Numerical results.
In this section, we examine two experimental design problems. The first is a linear Gaussian problem with a known, exact solution; this example is used to assess the performance and convergence properties of the layered multiple importance sampling algorithm, for varying dimensions. The second example is a more challenging nonlinear experimental design problem that occurs in Mössbauer spectroscopy, where we again illustrate the computational efficieny of the LMIS approach as well as the impact of the focused design objective.

4.1.
Higher-dimensional linear Gaussian example. We revisit the linear Gaussian problem that was first introduced in subsection 2.2, but increase the dimension of the inference problem to n = 4 and n = 8 dimensions. The dimension of the parameter of interest remains the same, i.e., n θ = 1, so the dimension of the nuisance parameters is n η = n − n θ . In general, the observation model for our n-dimensional inference problem is:  The design matrix G is parameterized by the parameter d ∈ [0, 1] and a gain parameter k > 0. The diagonal terms describe a trade-off between learning about the parameter of interest and learning about the nuisance parameters where d = 1 maximizes the information gain in θ at the expense of information gain in the nuisance parameters (η 1 , η 2 , . . . , η n−1 ). Furthermore, the nonzero off-diagonal terms in the design matrix introduce coupling between the parameter of interest and the nuisance parameters so that the resulting optimal design is not trivially d = 1.

Algorithm 1: Layered multiple importance sampling estimate of the expected information gain
Function ExpectedInformationGain(d) for i ← 1 . . . N do Sample (θ (i) , η (i) ) from the prior p(θ, η) Sample y (i) from the likelihood p(y|θ (i) , η (i) , d) The EIG for all of the parameters is where Γ post is the posterior covariance matrix, given by // Compute importance weights and normalize them // Estimate posterior mean and covariance − µ  Figure 4. We see that the optimal designs for the two objectives are very different: the optimal design for learning about all the parameters is d * joint = 0 and the optimal design for learning only about the parameter of interest is d * marg ≈ 0.93. The determinant of the posterior covariance matrix is plotted in Figure 4 as a function of the design parameter d, which illustrates the relationship between the posterior covariance and the EIG.
From the analysis of the EIG estimator in (3.8) and (3.9), we expect that the bias and variance will be largest when the biasing distributions and the posterior distributions have the smallest overlap. Thus, if we use the the prior distribution as the biasing distribution, the bias and variance of the EIG estimator should be large when there is significant posterior concentration. Accordingly, in Figure 4 we observe that the volume (computed as the determinant of the posterior covariance matrix) of the posterior distribution varies over several orders of magnitude across the design space, and in Figure 5a and Figure 5b we observe large positive bias for designs that correspond to significant posterior concentration.
The correlation of estimator bias with posterior concentration is particlularly undesirable behavior-particularly in the context of focused experimental design, since bias that is correlated with the determinant of the joint posterior can mask the true shape of the objective function (which depends only on the marginal posterior). We observe this behavior in Figure 5b where the location of local and global maxima correspond to sub-optimal designs that maximize posterior concentration instead of those that maximize information gain in the parameter of interest. Furthermore, when using the prior as a biasing distribution, the large estimator variance may pose a challenge for stochastic optimization algorithms. On the other hand, LMIS yields estimates of the EIG with much smaller bias and variance while using the same number of model evaluations. More importantly, the LMIS estimates of the EIG in both the focused and unfocused cases capture the correct locations of the global maxima.

Sample allocations, bias, and variance.
The allocation of samples between the outer and inner estimators can play an important role in the bias and variance of the overall estimator. In Figure 6, we investigate the effect of the ratio between the number of outer and inner samples on the mean squared error of the EIG estimator, both using the prior as a biasing distribution and using layered multiple importance sampling. When using the prior as a biasing distribution, the bias dominates the mean squared error, so allocating more inner samples-at the expense of affording fewer outer samples-will reduce the overall mean squared error. On the other hand, with LMIS, the biasing distributions closely approximate the posterior distributions, so fewer inner samples are required, and more samples can be allocated for the outer estimator, which not only decreases the overall variance, but also provides additional opportunities for incremental enrichment, further improving the approximation accuracy of the inner estimators.
Using favorable but fixed ratios between the number of outer and inner samples, 2 we can examine the mean squared error as a function of number of model evaluations. In Figure 7, we observe a dramatic reduction in the mean squared error when using LMIS. With 10 6 model evaluations, the resulting LMIS estimate has mean squared error comparable to using the exact posteriors as biasing distributions. 3 To better understand the observed convergence rates in Figure 7, we plot the mean squared error versus the number of inner and outer samples in Figure 8. When using the prior as the biasing distribution, we observe the convergence rates predicted by (3.8) and (3.9). The estimator bias decreases as 1/M and is unaffected by the number of outer samples, N . The variance also decreases with  increasing number of inner and outer samples, which suggests that the 1/N M terms in (3.9) are not negligible.
On the other hand, LMIS allows the estimator bias to decrease with the number of outer samples. With larger N , more samples can initially be used to estimate the posterior moments used to find optimal biasing distributions. Furthermore, the LMIS estimator applies incremental enrichment a total of N times, so as N becomes larger, more samples are made available for use when estimating the posterior moments. When using LMIS, the variances of the inner estimators are significantly reduced, so the 1/N M terms in the variance for the LMIS estimator are negligible in comparison to prior-based sampling. Thus, instead of the usual bias-variance trade-off between maximizing the number of inner and outer samples, the earlier analysis and the numerical results suggest that maximizing the number of outer samples is the best way to reduce both the bias and variance! In fact, if the biasing distributions converge to the posterior distributions as N → ∞ (e.g., if they are of the same family of distributions), it is conceivable that one could set M = 1 in this limit, since importance samping using the normalized posterior results in a zero-variance estimator. However, for finite sample sizes and non-Gaussian posterior distributions, we still require M → ∞ to ensure asymptotic consistency, since the inner estimators correct for the mismatch between the biasing distributions and the posterior distributions.

Importance sampling diagnostics and incremental enrichment.
The exact value of the EIG is generally not available, so in practice, we must instead rely on importance sampling diagnostics to measure the performance of the estimator. A useful "rule of thumb" for measuring importance sampling performance is to use the effective sample size (ESS) to measure how different the biasing distribution is from the target distribution [28]. However, in our case, the target distribution is the prior distribution, so maximizing the ESS diagnostic would suggest choosing the prior as the biasing distribution, which we showed to be sub-optimal. Instead, we consider the effective sample size customized for the target function-the marginal or conditional likelihood-as suggested in [33]. For a diagnostic that is specific to the importance sampling estimator in (3.2), we define w (i,j) .
Then, the effective sample sizes customized to the marginal likelihood and the conditional likelihood are In Figure 9, we can see that LMIS greatly increases the effective sample size customized to the marginal and conditional likelihood for both small N = 320 and large N = 1000 values of N and M .
In Figure 10, we explore the mechanism by which incremental enrichment can improve the quality of the biasing distributions for the inner estimators. Since the samples of (θ (i) , η (i) ) and the corresponding y (i) are sorted in decreasing order of prior density p(θ (i) , η (i) ), we should expect that the effective sample size of the inner estimators decrease with increasing outer loop index i. This is particularly evident when using the prior as a biasing distribution, since as i increases, the ith posterior p(θ, η |y (i) , d) will be further out in the tail of the prior distribution p(θ, η), and is reflected in the monotonic decrease of cESS with i seen in Figure 10. On the other hand, when using LMIS, the cESS initially increases for small i, as a direct result of incremental enrichment. Then, the cESS decays more slowly for the bulk of the inner estimates than in the prior case. Finally, for large i, the quality of the biasing distributions deteriorates rapidly, since there are fewer available samples for reuse when estimating the moments of posterior distributions concentrated in the tail of the prior. However, as N becomes larger, the impact of the deterioration in the tails diminishes (cf. the decrease in probability mass in the lower-right corner of each subplot in Figure 10). The deterioration in the tails is also reflected in Figure 8 where the bias of the LMIS does not decrease as 1/M for fixed N . This can be attributed to the observation that the biasing distributions q (i) CL for estimating the conditional likelihoods p(y (i) |θ, d) are obtained by calculating conditionals of an estimated posterior covariance. The conditional mean and variance are very sensitive to the correlation in the estimated posterior covariance, which may not be a robust estimator in the small sample size limit. The result is that the variance of the inner estimator of the conditional likelihood may decrease more slowly than 1/M .

4.2.
Nonlinear example: Mössbauer spectroscopy. The Mössbauer effect refers to recoil-free nuclear resonance fluorescence, which involves the resonant and recoil-free emission and absorption of gamma radiation by atomic nuclei bound in a solid [30]. Mössbauer spectroscopy is responsible for the discovery of the Mössbauer isomeric shift [22], which carries valuable information regarding the electronic and chemical structure of the lattice in which emitting and absorbing nuclei are embedded. These shifts are incredibly small, usually on the order of 10 −8 eV, which is roughly 12 orders of magnitude smaller than the energy of the emitted gamma rays themselves. Spectroscopic experiments involve moving the gamma ray source relative to the absorber. Change in the energy of the emitted gamma ray from Doppler shifting can be precisely controlled by the relative velocity of the source and absorber. In practice, shifts of 10 −8 eV correspond to velocities on the order of mm/s. The general experimental setup for Mössbauer spectroscopy involves a radioisotope source, a drive system to move the source (either with constant acceleration or constant velocity), an absorber, and a gamma ray detector. The Mössbauer drive is often operated at constant velocity for precision measurements, so that integration time may be spent only at the velocities of interest, which, in general, are not uniformly or symmetrically distributed about zero [16]. If one uses long counting times and small velocity steps to obtain high precision, the total time of the run can take days, so careful selection of velocities at which to acquire data is critical for a successful experiment.

Problem formulation.
An experimental design consists of a set of n d velocities for high precision measurement of the isomeric shift-the center of the absorption peak-which we denote by θ. We use the standard parameterization of the absorption peak as having a Lorentzian profile [18], such that the the number of detector counts y i at a velocity d i can be described by the following relationship, also illustrated in Figure 11: The model parameters are assigned the following prior distributions, where log-normal priors are assigned to parameters which can only have positive values, and the mean and variance of the prior distributions are chosen to "standardize" the signal, which would be standard experimental practice before collecting data. The additive observation error is assumed to be Gaussian where i ∼ N (0, 0.1 2 ), which represents a signal-to-noise ratio of 10:1.
In Figure 12, we plot two "maps" of EIG estimated using LMIS, using sample sizes M 1 = M 2 = 100 and N = 1000. One map uses θ = 'center' and lets η comprise the remaining parameters; the other map puts θ = 'offset' and lets η comprise the rest. The optimal designs differ for these two scenarios: d * center = (−1.3, 0, 1.3) and d * offset = (−2, 0, 2). For visualization purposes, we restrict the design problem to three velocities, with one of the three velocities fixed to zero. To better understand the impact of these design choices, Figure 13 plots the posterior distribution obtained using a data realization from each optimal design. Comparing the two focused designs, we observe a distinct tradeoff between learning about the 'center' parameter and the 'offset' parameter, especially in the lower-left pairwise marginal. Furthermore, the two-dimensional marginals clearly illustrate the non-Gaussianity of the posterior distribution.
Next we evaluate the performance of the LMIS estimator, focusing on the estimation of EIG in the marginal distribution of θ = 'center.' Since we do not have a closed-form expression for this EIG, we resort to calculating a "reference" value by using a large computational budget: 10 8 model evaluations with the LMIS estimator. In Figure 14, we plot the mean squared error with respect to this reference value for both the LMIS estimator and the prior estimator. We choose to evaluate the focused EIG for the design (a) Expected utility for learning 'center' with optimal designs marked by stars and illustrated above. (b) Expected utility for learning 'offset' with optimal designs marked by stars and illustrated above.  3) since it presents the most challenging case-with the greatest amount of posterior concentration. Compared to the prior estimator, we observe a reduction of up to two orders of magnitude in the mean squared error when using LMIS. This result is consistent with the earlier linear-Gaussian example and suggests that the algorithm is able to identify high-quality biasing distributions even in the nonlinear case.
We can also use the cESS diagnostic (4.2) to evaluate the performance of the LMIS scheme. The distribution of cESS is plotted in Figure 15, and displays a marked increase over that obtained when using the prior as a biasing distribution, comparable to the improvement for the linear Gaussian example of the same dimension. Again, this result suggests that LMIS is able to find high-quality biasing distributions for estimating both the marginal and conditional likelihoods, even in the case of a nonlinear observational model and non-Gaussian posteriors.

5.
Conclusion. This paper introduces a new computational approach for focused optimal Bayesian experimental design, where the design objective is to maximize expected information gain in the marginal distribution of parameters of interest. In particular, we develop a consistent layered multiple importance sampling (LMIS) estimator of this objective, which uses incremental enrichment to construct a sequence of increasingly tailored biasing distributions. We demonstrate the effectiveness of our approach on design problems of moderate dimension. Compared to previous algorithms that draw samples from the prior, the LMIS estimator reduces both bias and variance by several orders of magnitude, for the same number of model evaluations. The reduction in variance may improve the efficiency of optimization over the design space, but more importantly, the reduction in bias diminishes the confounding effect caused by link between estimator bias and posterior concentration. The latter is particularly important in the context of focused design, as posterior concentration in general has a different dependence on the design parameters than the information gain in selected marginals of interest.
The enrichment process of LMIS suggests a different balance of inner-loop and outer-loop samples than previous nested Monte Carlo methods. We also emphasize that LMIS does not require gradients of the likelihood or deterministic models embedded therein, and thus is well suited to problems where gradient information may not be available. While here we implemented LMIS with a particular family of inner-loop biasing distributions (multivariate t), the overall algorithm can be used with any parametric family or mixture, as suited to Possible extensions of this work include using surrogate models to accelerate estimation; while previous work did so in a biased way [19], the LMIS framework could incorporate surrogate models only in its outer layer, to help construct biasing distributions more cheaply. When gradients of the posterior density are readily available, one also could use Laplace approximations of the posterior distribution to construct suitable biasing distributions for estimation of the marginal [4] and conditional likelihoods. Moreover, surrogate models and approximations of the marginal likelihood could be combined with the full-model LMIS estimator using a control variate or multilevel Monte Carlo approach, to further reduce the number of full model evaluations while preserving consistency. Another clear extension is to pair LMIS estimators with optimization methods designed for stochastic objectives (e.g., [20,40,32,8]) to facilitate search over the design space.
An interesting application of focused experimental design is for models that contain discrepancy terms [21] intended to capture model misspecification. In this setting, one could treat parameters of the discrepancy model as nuisance parameters (η), where the goal is to learn about other parameters of interest (θ) in the presence of model discrepancy. Alternatively, one could focus on learning about the discrepancy term to better capture the nature of the mismatch between the model and data.
where y (i) ∼ p(y (i) |θ (i) , η (i) , d), (θ (i) , η (i) ) ∼ p(θ, η), and S  and from the law of total variance, the variance of U (d) can be written as Using the delta method, we can obtain the second-order approximation of the expected value and variance of u(y, θ, S A , S B ; d):   .7) is consistent (see [34]) and in particular, its bias and variance converge to zero as M 1 , M 2 , N → ∞. Since all three sample sizes must increase, it is natural to ask whether there is an optimal allocation of inner (M 1 , M 2 ) and outer (N ) samples for a given computational budget W = N (M 1 + M 2 ). Let α 2 = 2M/N be the ratio between the number of inner and outer Monte Carlo samples. Then, by substituting M = N α 2 /2 and N = √ W /α, we can rewrite the MSE in terms of the total computational work W and α: The value of α and the corresponding ratio α 2 that minimize the leading-order MSE are This result suggests that the ratio between the number of inner and outer samples should decrease slowly as the computational budget increases, in particular scaling as W −1/3 . The ratio C/D 3 reflects the relative magnitudes of leading-order terms in the bias and variance; if the bias dominates, one should choose a relatively larger number of inner loop samples. The minimum MSE, realized by this value of α * , is Of course, if one does not know the constant in (B.2) and instead just uses α = α 0 /W 1/6 , the work-scaling of the MSE remains the same: Thus we can compare two scenarios. For a fixed ratio between M and N , i.e., if α 2 were not allowed to depend on W , the leading-order term of (B.1) would yield MSE = O(W −1/2 ). If, instead, we adjust the ratio according to the optimal scaling α 2 ∝ W −1/3 , we obtain a faster decay rate of MSE = O(W −2/3 ). We use the 4-dimensional linear Gaussian problem to verify the results of the optimal scaling analysis in subsection 3.2. In Figure 16, we plot the estimator MSE against the number of model evaluations using two scaling strategies: a fixed ratio (dashed line) and a scaled ratio (solid line). In practice, the constant term ( C D 3 ) 2/3 in (B.2) for the optimal ratio is not known. In Figure 16, we anchor the ratio at W = 10 4 to be M = 10N (left) and M = 1000N (right) to investigate the influence of the constant term on the convergence rate for finite sample sizes. We observe that when the MSE is dominated by the estimator bias (e.g., M = 10N ), the trade-off explored by the optimal scaling actually results in sub-optimal convergence in MSE. This is because the benefits of a faster convergence rate in the estimator variance are overshadowed by the slower convergence rate in the bias. On the other hand, when the bias and variance are of similar magnitude (e.g., M = 1000N at W = 10 4 ), we can observe benefits of a faster convergence rate. Between W = 10 6 and W = 10 7 , the MSE converges at the predicted rates: W −1/2 when using a fixed ratio and W −2/3 when using the optimal scaling. In practice, without knowing the optimal value of the constant from (B.2) or the magnitude of the estimator bias, a strategy that uses a fixed ratio between the number of samples in the inner and outer estimators may be appropriate for realistic computational budgets.