Generalized Multiple Importance Sampling

Introduction

Importance sampling (IS) is a well-known Monte Carlo technique that can be applied to compute integrals involving target probability density functions (pdfs) [Robert and Casella 2004; Liu 2004]. The standard IS technique draws samples from a single proposal pdf and assigns them weights based on the ratio between the target and the proposal pdfs, both evaluated at the sample value. The choice of a suitable proposal pdf is crucial for obtaining a good approximation of the target pdf using the IS method. Indeed, although the validity of this approach is guaranteed under mild assumptions, the variance of the estimator depends on the discrepancy between the shape of the proposal and the target [Robert and Casella 2004; Liu 2004].

Several advanced strategies have been proposed in the literature to design more robust IS schemes [Liu 2004, Chapter 2], [Owen 2013, Chapter 9], [Liang 2002]. A powerful approach is based on using a population of different proposal pdfs. This approach is often referred to as multiple importance sampling (MIS) and several possible implementations have been proposed depending on the specific assumptions of the problem, e.g. the knowledge of the normalizing constants, prior information of the proposals, etc [Veach and Guibas 1995; Hesterberg 1995; Owen and Zhou 2000; Tan 2004; He and Owen 2014; Elvira et al. 2015a]. In general, MIS strategies provide more robust algorithms, since they avoid entrusting the performance of the method to a single proposal. Moreover, many algorithms have been proposed in order to conveniently adapt the set of proposals in MIS [Cappé et al. 2004; Martino et al. 2017; Elvira et al. 2017].

When a set of proposal pdfs is available, the way in which the samples can be drawn and weighted is not unique, unlike the case of using a single proposal. Indeed, different MIS algorithms in the literature (both adaptive and non-adaptive) have implicitly and independently interpreted the sampling and weighting procedures in different ways [Owen and Zhou 2000; Cappé et al. 2004; Cappé et al. 2008; Elvira et al. 2015a; Martino et al. 2015; Cornuet et al. 2012; Bugallo et al. 2017]. Namely, there are several possible combinations of sampling and weighting schemes, when a set of proposal pdfs is available, which lead to valid MIS approximations of the target pdf. However, these different possibilities can largely differ in terms of performance of the corresponding estimators.

In this paper, we introduce a unified framework for MIS schemes, providing a general theoretical description of the possible sampling and weighting procedures when a set of proposal pdfs is used to produce an IS approximation. Within this unified context, it is possible to interpret that all the MIS algorithms draw samples from an equally-weighted mixture distribution obtained from the set of available proposal pdfs. Three different sampling approaches and five different functions to calculate the weights of the generated samples are proposed and discussed. Moreover, we state two basic rules for possibly devising new valid sampling and weighting strategies within the proposed framework. All the analyzed combinations of sampling/weighting provide consistent estimates of the parameters of interest.

The proposed generalized framework includes all of the existing MIS methodologies that we are aware of (applied within different algorithms, e.g. in [Elvira et al. 2015a; Cappé et al. 2004; Cornuet et al. 2012; Martino et al. 2015; Martino et al. 2017; Elvira et al. 2017]) and allows the design of novel techniques (here we propose three new schemes, but more can be introduced). An exhaustive theoretical analysis is provided by introducing general expressions for sampling and weighting in this generalized MIS context, and by proving that they yield consistent estimators. Furthermore, we compare the performance of the different MIS schemes (the proposed and existing ones) in terms of the variance of the estimators.

The rest of this paper is organized as follows. In Section 2, we describe the problem and we revisit the standard IS methodology. In Section 3, we discuss the sampling procedure in MIS, propose three new sampling strategies, and analyze some distributions of interest. In Section 4, we propose five different weighting functions, some of them completely new, and show their validity. The different combinations of sampling/weighting strategies are analyzed in Section 5, establishing the connections with existent MIS schemes, and describing three novel MIS schemes. In Section 6, we analyze the performance of the different MIS schemes in terms of the variance of the estimators. Then, Section 7 discusses some relevant aspects about the application of the proposed MIS schemes, including their use in adaptive settings. Finally, Section 8 presents some descriptive numerical examples where the different MIS schemes are simulated, and Section 9 contains some concluding remarks. A running example is introduced in Section 3 and continued in Section 4, Section 5, Section 6 and Section 8 in order to clarify the flow of the paper. In addition, we perform numerical simulations on the running example, where the proposal pdfs are intentionally well chosen, to evidence the significant effects produced by the different interpretations of the sampling and weighting schemes.

Problem Statement and Background

IS is a general Monte Carlo technique for the approximation of a pdf of interest by a random measure composed of samples and weights [Robert and Casella 2004]. In its original formulation, a set of $N$ samples, $\{{\bf x}_{n}\}_{n=1}^{N}$, is drawn from a single proposal pdf, $q({\bf x})$, with heavier tails than those of the target pdf, $\pi({\bf x})$. A particular sample, ${{\bf x}_{n}}$, is assigned an importance weight given by

which represents the ratio between the target pdf, $\pi$, and the proposal pdf, $q$, both evaluated at ${\bf x}_{n}$. The samples and weights form the random measure $\chi=\{{\bf x}_{n},w_{n}\}_{n=1}^{N}$ that approximates the measure of the target pdf as

where $\delta_{{\bf x}_{n}}({\bf x})$ is the unit delta measure concentrated at ${\bf x}_{n}$ and $\hat{Z}=\frac{1}{N}\sum_{j=1}^{N}w_{j}$ is an unbiased estimator of $Z=\int\pi({\bf x})d{\bf x}$ [Robert and Casella 2004]. Fig. 1 (a) displays an example of a target pdf and a proposal pdf, as well as the samples and weights that form a random measure approximating the posterior. Note that, unlike Markov Chain Monte Carlo (MCMC) methods, all the generated samples are used to build the estimators, e.g., there is no burn-in period.

2 Estimators in importance sampling

Otherwise, if the target distribution is only known up to the normalizing constant, $Z$, one can use the self-normalized estimator

where $Z$ is approximated by the estimate

Sampling in Multiple Importance Sampling

MIS schemes consider a set of $N$ proposal pdfs, $\{q_{{n}}({\bf x})\}_{n=1}^{N}\equiv\{q_{1}({\bf x}),\ldots,q_{N}({\bf x})\}$, and proceed by drawing $M$ samples, $\{{\bf x}_{m}\}_{m=1}^{M}$ (where $M\neq N$, in general) and properly weighting them. As a visual example, Fig. 1 (b) displays a target pdf and two proposal pdfs, as well as the samples and weights that form a random measure approximating the posterior.

Unequal weights could also be considered in the mixture. In [He and Owen 2014], the weights can be optimized to minimize the variance for a certain integrand.

Let us consider a generic mechanism for the simulation of $N$ samples from the set of $N$ proposals. Starting with $n=1$:

Choose an index $j_{n}\in\{1,\ldots,N\}$, which corresponds to the selection of the proposal pdf $q_{j_{n}}$.

Generate a sample ${\bf x}_{n}$ from the selected proposal pdf, i.e., ${\bf x}_{n}\sim q_{j_{n}}({\bf x}_{n})$.

Note that, in step 1, the probabilities associated to each possible value of $j_{n}$ are not specified yet. The graphical model corresponding to this sampling scheme is shown in Fig. 2.

Therefore, obtaining the set of samples $\{{\bf x}_{n}\}_{n=1}^{N}\equiv\{{\bf x}_{1},...,{\bf x}_{N}\}$ is in general a two step sequential procedure. First, the $n$-th index $j_{n}$ is drawn according to some conditional pdf, $P(j_{n}|j_{1:n-1})$, where $j_{1:n-1}\equiv\{j_{1},\ldots,j_{n-1}\}$ is the sequence of the previously generated indexes. We use a simplified argument-wise notation, where $p({\bf x}_{n})$ denotes the pdf of the continuous random variable (r.v.) ${\bf X}_{n}$, while $P(j_{n})$ denotes the probability mass function (pmf) of the discrete r.v. $J_{n}$. Also, $p({\bf x}_{n},j_{n})$ denotes the joint pdf and $p({\bf x}_{n}|j_{n})$ is the conditional pdf of ${\bf X}_{n}$ given $J_{n}=j_{n}$. If the argument of $p(\cdot)$ is different from ${\bf x}_{n}$, then it denotes the evaluation of the pdf as a function, e.g., $p({\bf z}|j_{n})$ denotes the pdf $p({\bf x}_{n}|j_{n})$ evaluated at ${\bf x}_{n}={\bf z}$. Then, the $n$-th sample is drawn from the selected proposal pdf as ${\bf x}_{n}\sim p({\bf x}_{n}|j_{n})$. The joint probability distribution of the current sample and all the indexes used to generate the samples from $1$ to $n$ is

where $p({\bf x}_{n}|j_{n})=q_{j_{n}}({\bf x}_{n})$ is the $n$-th selected proposal pdf, $q_{j_{n}}({\bf x}_{n})$.

2 Selection of the proposal pdfs

In the sequel, we describe three mechanisms for obtaining the sequence of indexes, $j_{1:N}$. All the mechanisms share the property that

i.e., all the indexes have the same (marginal) probability of being selected.

Random index selection with replacement: The $N$ indexes are independently drawn from the set $\{1,\ldots,N\}$ with equal probability. Thus, we have

With this type of index sampling, there may be more than one sample drawn from some proposal, and there may be proposal pdfs that are not used at all.

Random index selection without replacement: The indexes are uniformly and sequentially drawn from different sets as $j_{1}\in\mathcal{I}_{1}=\{1,\ldots,N\}$, … , $j_{n}\in\mathcal{I}_{n}=\{1,\ldots,N\}\setminus\{j_{1:n-1}\}$, i.e., removing the proposals previously used. Hence, the conditional probability mass function (pmf) of the $n$-th index given the previous ones is now

where $|\mathcal{I}_{n}|=N-n+1$. Note that the marginal pmf of the $j$-th index is still given by (3.4). There are $N!$ equiprobable configurations (permutations) of the sequence $\{j_{1},\ldots,j_{N}\}$, and in $(N-1)!$ the $k$-th index is drawn at the $n$-th position $\forall k,n=1,\ldots,N$. Therefore, $P(J_{n}=k)=\frac{(N-1)!}{N!}=\frac{1}{N}$ $\forall k,n=1,\ldots,N$. However, exactly one sample is drawn from each of the proposal pdfs by following this strategy.

Deterministic index selection without replacement: This sampling is a particular case of sampling $\mathcal{S}_{2}$, where a fixed deterministic sequence of indexes is drawn. For instance, and without loss of generality: $j_{1}=1,j_{2}=2,\ldots,j_{n}=n,\ldots,j_{N}=N$. Therefore, ${\bf x}_{n}\sim q_{j_{n}}({{\bf x}_{n}})=q_{n}({\bf x}_{n})$, and the conditional pmf of the $n$-th index given the $n-1$ previous ones becomes

The connexions of the sampling mechanisms with some resampling schemes are discussed in Appendix B.

3 Running example

Let us consider $N=3$ Gaussian proposal pdfs $q_{1}(x)=\mathcal{N}(x;\mu_{1},\sigma_{1}^{2})$, $q_{2}(x)=\mathcal{N}(x;\mu_{2},\sigma_{2}^{2})$ and $q_{3}(x)=\mathcal{N}(x;\mu_{3},\sigma_{3}^{2})$ with predefined means and variances. In $\mathcal{S}_{1}$, a possible realization of the indexes is the sequence $\{j_{1},j_{2},j_{3}\}=\{3,3,1\}$. Therefore, in this situation, ${\bf x}_{1}\sim q_{3}$, ${\bf x}_{2}\sim q_{3}$, and ${\bf x}_{3}\sim q_{1}$. In $\mathcal{S}_{2}$, the realization could result from the permutation $\{j_{1},j_{2},j_{3}\}=\{3,1,2\}$. In $\mathcal{S}_{3}$, the sequence is deterministically obtained as $\{j_{1},j_{2},j_{3}\}=\{1,2,3\}$.

4 Distributions of interest of the nn-th sample, 𝐱n{\bf x}_{n}

In the following, we discuss some important distributions related to the set of samples drawn. These distributions are of utmost importance to understand the different methods for weighting the samples discussed in the following section.

Note that the distribution of the $n$-th sample given all the knowledge of the process up to that point is $p({\bf x}_{n}|j_{1:n-1},{\bf x}_{1:n-1})=p({\bf x}_{n}|j_{1:n-1})$. In $\mathcal{S}_{1}$, this distribution corresponds to $p({\bf x}_{n}|j_{1:n-1})=\psi({{\bf x}_{n}})$. We recall that $\psi$ is the mixture of proposals defined in Eq. (3.1). In $\mathcal{S}_{2}$, we have $p({\bf x}_{n}|j_{1:n-1})=\frac{1}{|\mathcal{I}_{n}|}\sum_{\forall k\in\mathcal{I}_{n}}q_{k}({\bf x})$. Finally, under $\mathcal{S}_{3}$, $p({\bf x}_{n}|j_{1:n-1})=q_{n}({\bf x}_{n})$. Once the $n$-th index $j_{n}$ has been selected, the $n$-th sample, ${\bf x}_{n}$, is distributed as $p({\bf x}_{n}|j_{n})=q_{j_{n}}({\bf x}_{n})$ in any sampling method within the proposed framework. The marginal distribution of this $n$-th sample, ${\bf x}_{n}$, is then given by

5 Distributions of interest beyond 𝐱n{\bf x}_{n}

The traditional IS approach focuses just on the distribution of the r.v. ${\bf X}_{n}$. In MIS, we are also interested in the statistical properties of the set of samples, regardless of their index $n$, since the $N$ samples are used jointly in the estimators, regardless their order of appearance. Hence, we introduce a generic r.v.,

where $\mathcal{U}\{1,2,\ldots,N\}$ is the discrete uniform distribution on the set $\{1,2,\ldots,N\}$. The density of ${\bf X}$ is then given by

where $p_{{\bf x}_{n}}({\bf x})$ denotes the marginal pdf of ${\bf X}_{n}$, given by Eq. (3.7), evaluated at ${\bf x}$, and $\psi({\bf x})$ is the mixture pdf. For the sake of clarity, in Eq. (3.9) we have used the notation $p_{{\bf x}_{n}}({\bf x})$, instead of $p({\bf x})$ as in Eq. (3.7) and the rest of the paper, to denote the marginal pdf of ${\bf X}_{n}$ evaluated at ${\bf x}$. Moreover, one can also obtain the conditional pdf of ${\bf X}$ given the sequence of indexes as

Note that, in this case, $f({\bf x}|j_{1:N})=\psi({\bf x})$ for the schemes without replacement at the index selection ($\mathcal{S}_{2}$ and $\mathcal{S}_{3}$), but $f({\bf x}|j_{1:N})=\frac{1}{N}\sum_{n=1}^{N}q_{j_{n}}({\bf x})$ for the case with replacement ($\mathcal{S}_{1}$), i.e., some proposal pdfs may not appear while others may appear repeated.

(Sampling): In the proposed framework, we consider valid, any sequential sampling scheme for generating the set $\{{\bf X}_{1},\ldots,{\bf X}_{N}\}$ such that the pdf of the r.v. ${\bf X}$ defined in Eq. (3.8) is given by $\psi({\bf x})$. Further considerations about the r.v. ${\bf X}$ and connections with variance reduction methods [Robert and Casella 2004; Owen 2013] are given in Appendix A.

Table 5 summarizes all the distributions of interest. Note that, the pdf of the r.v. ${\bf X}$ is always the mixture $\psi({\bf x})$, but different sampling procedures yield different conditional and marginal distributions that will be exploited to justify different strategies for calculation of the importance weights in the next section. Finally, the last row of the table shows the joint distribution $p({\bf x}_{1:N})$ of the variables ${\bf X}_{1},\ldots,{\bf X}_{N}$, i.e., $p({\bf x}_{1:N})=\prod_{n=1}^{N}\psi({{\bf x}_{n}})$ and $p({\bf x}_{1:N})=\prod_{n=1}^{N}q_{n}({{\bf x}_{n}})$ for $\mathcal{S}_{1}$ and $\mathcal{S}_{3}$, respectively. For $\mathcal{S}_{2}$,

with $\mathcal{I}_{n}=\{1,\ldots,N\}\setminus\{j_{1:n-1}\}$.

Weighting in Multiple Importance Sampling

Our approach is based on analyzing which weighting functions yield proper MIS estimators. We consider that the set of weighting functions $\{w_{n}\}_{n=1}^{N}$ is proper if

This is equivalent to imposing the restriction

where we use the joint distribution of indexes and samples from Eq. (3.2).

Here we present several possible functions $\varphi_{\mathcal{P}_{n}}$, that yield an unbiased estimator of ${I}$ according to Eq. (4.3). The different choices for $\varphi_{\mathcal{P}_{n}}$, used in the denominator of the weight $w_{n}=\frac{\pi({\bf x}_{n})}{\varphi_{\mathcal{P}_{n}}({\bf x}_{n})}$, come naturally from the sampling densities discussed in Section 3. More precisely, they correspond to the five different functions in Table 5 related to the distributions of the generated samples. From now on, $p(\cdot)$ and $f(\cdot)$, which correspond to the pdfs of ${\bf X}_{n}$ and ${\bf X}$ respectively, are used as functions and the argument represents a functional evaluation.

$\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{j_{1:n-1}}({\bf x}_{n})=p({\bf x}_{n}|j_{1:n-1})$

Since the sampling process is sequential, this option is of particular interest. It interprets the proposal pdf as the conditional density of ${\bf x}_{n}$ given all the previous proposal indexes of the sampling process.

$\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{j_{n}}({\bf x}_{n})=p({\bf x}_{n}|j_{n})=q_{j_{n}}({\bf x}_{n})$

It interprets that if the index $j_{n}$ is known, $\varphi_{\mathcal{P}_{n}}$ is the proposal $q_{j_{n}}$.

$\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=p({\bf x}_{n})$

It interprets that ${\bf x}_{n}$ is a realization of the marginal $p({\bf x}_{n})$. This is probably the most “natural” option (as it does not assume any further knowledge in the generation of ${\bf x}_{n}$) and is a usual choice for the calculation of the weights in some of the existing MIS schemes (see Section 5).

$\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{j_{1:N}}({\bf x}_{n})=f({\bf x}_{n}|j_{1:N})=\frac{1}{N}\sum_{k=1}^{N}q_{j_{k}}({\bf x}_{n})$

This interpretation makes use of the distribution of the r.v. ${\bf X}$ conditioned on the whole set of indexes (defined in Section 3.5).

$\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi({\bf x}_{n})=f({\bf x}_{n})=\frac{1}{N}\sum_{k=1}^{N}q_{k}({\bf x}_{n})$

This option considers that all the ${\bf x}_{n}$ are realizations of the r.v. ${\bf X}$ defined in Section 3.5 (see Appendix A for a thorough discussion of this interpretation).

Table 1 summarizes the discussed functions $\varphi_{\mathcal{P}_{n}}$. Although some of the selected functions $\varphi_{\mathcal{P}_{n}}$ may seem more natural than others, all of them yield valid estimators. The proofs can be found in Appendix C. Other proper weighting functions are described in Section 7.2.

2 Connection with Liu-properness of single IS

We consider the definition of properness by Liu [Liu 2004, Section 2.5] and we extend (or relax) it to the MIS scenario. Namely, Liu-properness in standard IS states that a weighted sample $\{{\bf x}_{n},w_{n}\}$ drawn from a single proposal $q$ is proper if, for any square integrable function $g$,

i.e., $w$ can be in any form as long as the condition of Eq. (4.4) is fulfilled. Note that, for a deterministic weight assignment, the only proper weights are the ones considered by the standard IS approach. Note also that the MIS properness is a relaxation of the one proposed by Liu, i.e., any Liu-proper weighting scheme is also proper a according to our definition, but not vice versa.

3 Running example

Here we follow the running example of Section 3.3. For instance, let us consider the sampling method $\mathcal{S}_{1}$ and let the realization of the indexes be the sequence $\{j_{1},j_{2},j_{3}\}=\{3,3,1\}$. Under the weighting scheme $\mathcal{W}_{2}$, the weights would be computed as $w_{1}=\frac{\pi({\bf x}_{1})}{q_{3}({\bf x}_{1})}$, $w_{2}=\frac{\pi({\bf x}_{2})}{q_{3}({\bf x}_{2})}$, and $w_{3}=\frac{\pi({\bf x}_{3})}{q_{1}({\bf x}_{3})}$. However, under $\mathcal{W}_{4}$, $w_{1}=\frac{\pi({\bf x}_{1})}{\frac{1}{3}\left(q_{1}({\bf x}_{1})+2q_{3}({\bf x}_{1})\right)}$, $w_{2}=\frac{\pi({\bf x}_{2})}{\frac{1}{3}\left(q_{1}({\bf x}_{2})+2q_{3}({\bf x}_{2}\right)}$, and $w_{3}=\frac{\pi({\bf x}_{3})}{\frac{1}{3}\left(q_{1}({\bf x}_{3})+2q_{3}({\bf x}_{3})\right)}$. Note that all weighing schemes require the same number of target evaluations (which are usually more expensive) but different numbers of proposal evaluations.

Multiple Importance Sampling Schemes

In this section, we describe the different possible combinations of the three sampling strategies considered in Section 3 and the five weighting functions devised in Section 4. Once combined, the fifteen possibilities only lead to six unique MIS methods. Three of the methods are associated to the sampling scheme with replacement ($\mathcal{S}_{1}$), while the other three methods correspond to the sampling schemes without replacement ($\mathcal{S}_{2}$ and $\mathcal{S}_{3}$). Table 6 summarizes the possible combinations of sampling/weighting and indicates the resulting MIS method within brackets. The six MIS methods are labeled either by an R (sampling with replacement) or with an N (sampling with no replacement). We remark that these schemes are examples of proper MIS techniques fulfilling Remarks 3.1 and 4.1.

In all R schemes, the $n$-th sample is drawn with replacement (i.e., $\mathcal{S}_{1}$) from the whole mixture $\psi$:

Sampling with replacement, $\mathcal{S}_{1}$, and weight denominator $\mathcal{W}_{2}$: For the weight calculation of the $n$-th sample, only the proposal selected for generating the sample is evaluated in the denominator.

Sampling with replacement, $\mathcal{S}_{1}$, and weight denominator $\mathcal{W}_{4}$: With the $N$ selected indexes $j_{n}$, for $n=1,...,N$, one forms a mixture comprising all the corresponding proposal pdfs. The weight calculation of the $n$-th sample considers this a posteriori mixture evaluated at the $n$-th sample in the denominator, i.e., some proposals might be used more than once while other proposals might not be used.

Sampling with replacement, $\mathcal{S}_{1}$, and weight denominator $\mathcal{W}_{1}$, $\mathcal{W}_{3}$, or $\mathcal{W}_{5}$: For the weight calculation of the $n$-th sample, the denominator applies the value of the $n$-th sample to the whole mixture $\psi$ composed of the set of initial proposal pdfs (i.e., the function in the denominator of the weight does not depend on the sampling process). This is the approach followed by the so called mixture PMC method [Cappé et al. 2008].

2 MIS schemes without replacement

In all N schemes, exactly one sample is generated from each proposal pdf. This corresponds to having a sampling strategy without replacement.

Sampling without replacement (random or deterministic), $\mathcal{S}_{2}$ or $\mathcal{S}_{3}$, and weight denominator $\mathcal{W}_{2}$ (for $\mathcal{S}_{2}$) or $\mathcal{W}_{1}$, $\mathcal{W}_{2}$, or $\mathcal{W}_{3}$ (for $\mathcal{S}_{3}$): For calculating the denominator of the $n$-th weight, the specific proposal used for the generation of the sample is used. This is the approach frequently used in particle filtering [Gordon et al. 1993] and in the standard PMC method [Cappé et al. 2004].

Sampling without replacement (random), $\mathcal{S}_{2}$, and weight denominator $\mathcal{W}_{1}$: This MIS implementation draws one sample from each proposal, but the order matters (it must be random) since the calculation of the $n$-th weight uses for the evaluation of the denominator the mixture pdf formed by the proposal pdfs that were still available at the generation of the $n$-th sample.

Sampling without replacement (random or deterministic), $\mathcal{S}_{2}$ or $\mathcal{S}_{3}$, and weight denominator $\mathcal{W}_{3}$, $\mathcal{W}_{4}$, or $\mathcal{W}_{5}$ (for $\mathcal{S}_{2}$), or $\mathcal{W}_{4}$ or $\mathcal{W}_{5}$ (for $\mathcal{S}_{3}$): In the calculation of the $n$-th weight, one uses for the denominator the whole mixture. This is the approach, for instance, of [Martino et al. 2015; Cornuet et al. 2012]. As shown in Section 6, this scheme has several benefits over the others.

Table 2 summarizes the six resulting MIS schemes and their references in literature, indicating the sampling procedure and weighting function that are applied to obtain the $n$-th weighted sample ${\bf x}_{n}$. We consider N1 and N3 associated to $\mathcal{S}_{3}$ (they can also be obtained with $\mathcal{S}_{2}$) since it is simpler than $\mathcal{S}_{2}$. All the different algorithms in the literature (as far as we know) correspond to one of the MIS schemes described above (see Section 7.3). Moreover, several new valid schemes have also appeared naturally (R1, R2, and N2), and new ones can be proposed within this framework.

3 Running example

Let us consider the example from Section 3.3 where the realizations of the sequence of indexes for the sampling schemes $\mathcal{S}_{1}$, $\mathcal{S}_{2}$, and $\mathcal{S}_{3}$ are respectively $\{j_{1},j_{2},j_{3}\}=\{3,3,1\}$, $\{j_{1},j_{2},j_{3}\}=\{3,1,2\}$, and $\{j_{1},j_{2},j_{3}\}=\{1,2,3\}$. Figure 3(a) shows the three first schemes of Table 2 related to the sampling with replacement, $\mathcal{S}_{1}$. The figure shows a possible realization of all MIS schemes with $M=N=3$ samples and pdfs. For the $n$-th sample, we show the set of available proposals, the index $j_{n}$ of the proposal pdf that was actually selected to draw the sample, the function $\varphi_{n}$, and the importance weight. Similarly, Figs. 3(b)-(c) depict the three schemes of Table 2 related to the sampling without replacement, $\mathcal{S}_{2}$ and $\mathcal{S}_{3}$, where exactly one sample is drawn from each available proposal.

Variance Analysis of the Schemes

Although the six different MIS schemes that appear in Section 5 yield the estimator $\hat{I}$ of Eq. (2.4) unbiased (see Appendix C), the performance of each of the possible obtained estimators can be dramatically different. In this section, we provide an exhaustive variance analysis of the MIS schemes presented in the previous section. The details of the derivations are in Appendix D.2. The estimators of the three methods with replacement present the following variances:

On the other hand, the variances associated to the estimators of the three methods with no replacement are

One of the goals of this paper is to provide the practitioner with solid theoretical results about the superiority of some specific MIS schemes. In the following, we state two theorems that relate the variance of the estimator with these six methods, establishing a hierarchy among them. Note that obtaining an IS estimator with finite variance essentially amounts to having a proposal with heavier tails than the target. See [Robert and Casella 2004; Geweke 1989] for sufficient conditions that guarantee this finite variance.

For any target distribution $\pi({\bf x})$, any square integrable function $g$, and any set of proposal densities $\{q_{n}({\bf x})\}_{n=1}^{N}$ such that the variance of the corresponding MIS estimators is finite,

For any target distribution $\pi({\bf x})$, any square integrable function $g$, and any set of proposal densities $\{q_{n}({\bf x})\}_{n=1}^{2}$ such that the variance of the corresponding MIS estimators is finite,

First, let us note that the scheme N3 outperforms (in terms of the variance) any other MIS scheme in the literature that we are aware of. Moreover, for $N=2$, it also outperforms the other novel schemes R2 and N2. While the MIS schemes R2 and N2 do not appear in Theorem 6.1, we hypothesize that the conclusions of Theorem 6.2 might be extended to $N>2$. The intuitive reason is that, regardless of $N$, both methods partially reduce the variance of the estimators by placing more than one proposal at the denominator of some or all the weights. A possible interpretation of the superiority of N3 is that it uses the whole mixture at the denominator of each weight, thus providing an exchange of information between all the proposals. This exchange of information is essential in multimodal scenarios, where the whole set of proposals, seen as a mixture, should mimic the whole target, but each proposal should adapt locally to the target. Since the variance of the IS weight depends on the mismatch of the target (numerator) w.r.t. the proposal (denominator), the use of the whole mixture in the denominator reduces the variance of the weight in general, and therefore, also the variance of the estimator (see the variance analysis in Appendix D). The scheme N3 goes a step further w.r.t. R3, drawing deterministically one sample from each mixand of $\psi({\bf x})$, which can be seen as drawing $N$ samples from the mixture $\psi({\bf x})$ with a modified version of stratified sampling, a well-known variance reduction technique (see Appendix A and [Owen 2013, Section 9.12]), which is also related to the residual resampling.

Here we focus on computing the exact variances of estimators related to the running example. We simplify the case study to $N=2$ proposals, for the sake for conciseness in the proofs. The proposal pdfs are then $q_{1}(x)=\mathcal{N}(x;\mu_{1},\sigma^{2})$ and $q_{2}(x)=\mathcal{N}(x;\mu,\sigma^{2})$ with means $\mu_{1}=-3$ and $\mu_{2}=3$, and variance $\sigma^{2}=1$. We consider a normalized bimodal target pdf $\pi(x)=\frac{1}{2}\mathcal{N}(x;\nu_{1},c_{1}^{2})+\frac{1}{2}\mathcal{N}(x;\nu_{2},c_{2}^{2})$ and set $\nu_{1}=\mu_{1}$, $\nu_{2}=\mu_{2}$, and $c_{1}^{2}=c_{2}^{2}=\sigma^{2}$. Then, both proposal pdfs can be seen as a whole mixture that exactly replicates the target, i.e., $\pi(x)=\frac{1}{2}q_{1}(x)+\frac{1}{2}q_{2}(x)$. This is the desired situation pursued by an AIS algorithm: Each proposal is centered at each target mode, and the scale parameters perfectly match the scale of the modes. The goal consists in estimating the normalizing constant with the six schemes described in Section 5. We use the $\hat{Z}$ estimator of Eq. (2.6) and the estimator $\hat{I}$ of (2.4) when $g=x$, both with $N=2$ samples. The closed-form variance expressions of the six schemes are presented in the following:

The variances of the estimators of the normalizing constant (true value $Z=\int\pi(x)dx=1$) are given by

The variances of the estimators of the target mean (true value $I=\int x\pi(x)dx=0$) are given by

The derivations can be found in Appendix D.4. We observe that, for a very simple bimodal scenario where the proposals are perfectly placed in the target modes, the schemes R3 and N3 present a good performance while the other schemes do not work.

Applying the MIS Schemes

Table 3 compares the total number of target and proposal evaluations in each MIS scheme. First, note that the estimators of any MIS scheme within the proposed general framework perform $N$ target evaluations in total. However, depending on the function $\varphi_{\mathcal{P}_{n}}$ used by each specific scheme at the weight denominator, a different number of proposal evaluations is performed. We see that R3 and N3 always require the largest number of proposal evaluations. In R2, the number of proposal evaluations is variable: although each weight evaluates $N$ proposals, some proposals may be repeated, whereas others may not be used.

In many relevant scenarios, the cost of evaluating the proposal densities is negligible compared to the cost of evaluating the target function. In this scenario, the MIS scheme N3 should always be chosen, since it yields a lower variance with a negligible increase in computational cost. For instance, this is the case in the Big Data Bayesian framework, where the target function is a posterior distribution with a large amount of data in the likelihood function. However, in some other scenarios, e.g. when the number of proposals $N$ is too large and/or the target evaluations are not very expensive, limiting the number of proposal evaluations can result in a better cost-performance trade off.

Unlike most MCMC methods, several strategies of parallelization can be applied in IS-based techniques. In the adaptive context, the adaptation of all proposals usually depends on the performance of all previous proposals, and therefore the adaptivity is the bottleneck of the parallelization. The six schemes proposed in this paper can be parallelized to some extent. Once all the proposals are available, the schemes R1, N1, R3, and N3 can draw and weight the $N$ samples in parallel, which represents a large advantage w.r.t. MCMC methods. In the schemes R2 and N2, the samples can be drawn independently, but the denominator of the weight cannot be computed in a parallel way. However, since the target evaluation in the numerator of the weights is fully parallelizable, the drawback of these schemes can be considered negligible for a small/medium number of proposals.

2 A priori partition approach

The extra computational cost of some MIS schemes occurs because each sample must be evaluated in more than one proposal $q_{n}$, or even in all of the available proposals (e.g. the MIS scheme N3). In order to limit the number of proposal evaluations, let us first define a partition of the set of the indexes of all proposals, $\{1,\ldots,N\}$, into $P$ disjoint subsets of $L$ elements (indexes), $\mathcal{J}_{p}$ with $p=1,\ldots,P$, s.t.

After this a priori partition, one could apply any MIS scheme in each (partial) subset of proposals, and then perform a suitable convex combination of the partial estimators. This general strategy is inspired by a specific scheme, partial deterministic mixture MIS (p-DM-MIS), which was recently proposed in [Elvira et al. 2015a]. That work applies the idea of the partitions just for the MIS scheme N3, denoted there as full deterministic mixture MIS (f-DM-MIS). The sampling procedure is then $\mathcal{S}_{3}$, i.e., exactly one sample is drawn from each proposal. The weight of each sample in p-DM-MIS, instead of evaluating the whole set of proposals (as in N3), evaluates only the proposals within the subset that the generating proposal belongs to. Mathematically, the weights of the samples corresponding to the $p$-th mixture are computed as

Note that the number of proposal evaluations is $N\leq\frac{N^{2}}{P}\leq N^{2}$. Specifically, we have the particular cases $P=1$ and $P=N$ corresponding to the MIS schemes N3 (best performance) and N1 (worst performance), respectively. In [Elvira et al. 2015a], it is proved that for a specific partition with $P$ subsets of proposals, merging any pair of subsets decreases the variance of the estimator $\hat{I}$ of Eq. (2.4).

The previous idea can be applied to the other MIS schemes presented in Section 5 (not only N3). In particular, one can make an a priori partition of the proposals as in Eq. (7.1), and apply independently any different MIS scheme in each set. For instance, and based on some knowledge about the performance of the different proposals, one could make two disjoint sets of proposals, applying the MIS scheme N1 in the first set, and the MIS scheme N3 in the second set. Recently, a novel partition approach has been proposed in [Elvira et al. 2016]. In this case, the sets of proposals are performed a posteriori, once the samples have been drawn. The variance of the estimators is reduced at the price of introducing a bias.

3 Generalized Adaptive Multiple Importance Sampling

The MIS schemes considered in Section 5 can be directly applied to the adaptive context. Moreover, the a priori partition approach of Section 7.2 can be very useful to limit the computational cost of the different MIS schemes when the number of iterations grows (and therefore also the total number of proposals).

Let us assume that, at the $t$-th iteration, one sample is drawn from each proposal $q_{j,t}$ (sampling $\mathcal{S}_{3}$), i.e.,

$j=1,\ldots,J$ and $t=1,\ldots,T$. Then, an importance weight $w_{j,t}$ is assigned to each sample ${\bf x}_{j,t}$. As described exhaustively in Section 4, several strategies can be applied to build $w_{j,t}$ considering the different MIS approaches. Figure 4 provides a graphical representation of this scenario, by showing both the spatial and temporal evolution of the $J=NT$ proposal pdfs. In a generic AIS algorithm, one weight

is associated to each sample ${\bf x}_{j,t}$. In the MIS scheme N1, the function employed in the denominator is

In the following, we focus on the MIS scheme N3 in the adaptive framework, considering several choices of the partitioning of the set of proposals, since this scheme attains the best performance, as shown in Section 6. In the full N3 scheme, the function $\varphi_{j,t}$ is

where $\psi({\bf x})$ is now the mixture of all the spatial and temporal proposal pdfs. This case corresponds to the blue rectangle in Fig. 4. However, note that the computational complexity can become prohibitive as the product JT increases. Furthermore, two natural alternatives of partial N3 schemes appear in this scenario. The first one uses the following partial mixture

with $j=1,\ldots,J$, as mixture-proposal pdf in the IS weight denominator, i.e. using the temporal evolution of the $j$-th single proposal $q_{j,t}$ at the weight denominator. In this case, there are $P=J$ mixtures, each one formed by $L=T$ components (red rectangle in Fig. 4). Another possibility is considering the mixture of all the $q_{j,t}$’s at the $t$-th iteration, i.e.,

with $t=1,\ldots,T$, so that we have $P=T$ mixtures, each one formed by $L=J$ components (green rectangle in Fig. 4). The function $\varphi_{j,t}$ in Eq. (7.4) is used in the standard PMC scheme [Cappé et al. 2004]; Eq. (7.6), in the particular case of $J=1$, has been considered in the adaptive multiple importance sampling (AMIS) algorithm [Cornuet et al. 2012]. Note that the schemes that consider at the denominator of the weight the temporal sequence of adapted proposals can introduce a bias in the IS estimators (see [Cornuet et al. 2012, Section 5] for more details). The choice in Eq. (7.7) has been applied in the adaptive population importance sampling (APIS) [Martino et al. 2015], the layered adaptive importance sampling (LAIS) [Martino et al. 2017], and the deterministic mixture population Monte Carlo (DM-PMC) [Elvira et al. 2017] algorithms. In other techniques, such as mixture PMC (M-PMC) [Douc et al. 2007a; Douc et al. 2007b; Cappé et al. 2008], a similar strategy is employed, but using sampling $\mathcal{S}_{1}$ in the mixture $\phi_{t}({\bf x})$, i.e., with the MIS scheme R3.

Note that using $\psi({\bf x})$ and $\xi_{j}({\bf x})$ the computational cost per iteration increases as the total number of iterations $T$ grows. Indeed, at the $t$-th iteration all the previous proposals $q_{j,1},\ldots,q_{j,t-1}$ (for all $j$) must be evaluated at all the new samples ${\bf x}_{j,t}$. Hence, algorithms based on these proposals quickly become unfeasible as the number of iterations grows. On the other hand, using $\phi_{t}({\bf x})$ the computational cost per iteration is controlled by $J$, remaining constant regardless of the number of adaptive steps performed.

4 Guidelines for applying MIS

The superiority of N3 is theoretically proved for the unnormalized estimator in Theorems 6.1 and 6.2, and practically shown by means of several numerical simulations for the self-normalized estimator (see next section). However, the associated computational complexity is also increased w.r.t. the other MIS schemes in terms of proposal evaluations. If $N$ is small or the target evaluations are expensive (w.r.t. the cost of the proposal evaluations), N3 should be used. However, when the target evaluation is cheap and/or the number of proposals is large, the use of N3 increases notably the computational complexity. In this case, the novel schemes R2 or N2 seem to provide very good results, and their theoretical properties are superior to those of N1 and R1. However, future studies will be required to characterize these novel schemes and investigate efficient parallelization techniques. We also recommend to combine adaptive schemes with the partition approach proposed in [Elvira et al. 2015a] and [Elvira et al. 2016], and summarized in Section 7.2. Note that further investigation is also needed for efficiently constructing the partitions of the proposals that allow to reduce the computational complexity while retaining most of the variance reduction associated to the N3 scheme.

In the adaptive context, there is a big potential for the MIS schemes where all spatial and temporal proposals are used at the denominator of all weights (blue square in Fig. 4). However, the computational complexity for large number of proposals is prohibitive, and further theoretical analysis about the bias of the estimators is needed (see [Cornuet et al. 2012, Section 5]). The adaptivity of MIS algorithms is essential in challenging high-dimensional setups. The N3 scheme has exhibited a very good performance when used within adaptive MIS algorithms due to two main reasons. First, the variance of the estimators at each iteration is reduced as proved in Theorems 6.1 and 6.2, which explains part of the variance reduction attained in AMIS [Cornuet et al. 2012], LAIS [Martino et al. 2017], or GAPIS [Elvira et al. 2015b]. Second, when the IS weights are used for adaptive purposes (e.g. in APIS [Martino et al. 2015] or DM-PMC [Elvira et al. 2017]), the use of the whole mixture of proposals in the denominator of the weights can be seen as a cooperative adaptive procedure (see [Elvira et al. 2017] for further details).

Finally, one of the strengths of the N3 scheme is its performance in multimodal scenarios, where N1 should always be avoided. If $N$ is comparable to the number of modes, an adaptive N3 scheme should be employed; the aforementioned cooperation in the proposals adaptation has an implicit repulsive behavior that promotes the adaptation to different modes. However, if $N$ is much larger, the adaptive algorithm may use R2 or N2 with potentially similar performance but less computational complexity.

Numerical Examples

In the previous sections, we have provided several theoretical results for comparing different MIS schemes according to different quality measures, e.g., ranking them in terms of the variance of the corresponding estimators. In this section, we provide different numerical results in order to quantify numerically the gap among these methods. In the following, we show that even in the case where the different proposals are well tuned (in the sense of a small or no mismatch with a multimodal target), the choices of the sampling and weighting procedures dramatically affect the performance of the MIS estimator.

Let us consider again the target pdf of the running example

with means $\nu_{1}=-3$, $\nu_{2}=0$, and $\nu_{3}=3$, and variances $c_{1}^{2}=c_{2}^{2}=c_{3}^{2}=1$. As proposal functions we use $q_{i}(x)=\mathcal{N}(x;\mu_{i},\sigma)$, with $\mu_{i}=\nu_{i}$ and $i=1,2,3$ and $\sigma^{2}=1$, i.e., the proposal pdfs can be seen as a whole mixture that exactly replicates the target, i.e., $\pi(x)=\psi(x)=\frac{1}{3}q_{1}(x)+\frac{1}{3}q_{2}(x)+\frac{1}{3}q_{3}(x)$.

The goal is to estimate the mean of the target pdf with the six MIS schemes. Fig. 5(a) shows the MSE of the estimator $\hat{I}$ for all the methods w.r.t. the number of total samples (note that some schemes require that the total number of samples is multiple of $M=3$). The results have been averaged over $5\cdot 10^{6}$ runs. The solid black line shows the variance of the natural estimator, i.e. sampling directly from the target pdf (since this is possible in this easy example). Note that the method $\hat{I}_{\texttt{R3}}$ exactly replicates the performance of $\bar{I}$: this method samples from the mixture of Gaussians in the traditional way and the weights, due to the perfect match, are always $w=1$, i.e., $\hat{I}_{\texttt{R3}}$ and $\bar{I}$ are equivalent. We can see that $\hat{I}_{\texttt{N3}}$ is the best estimator in terms of variance, while $\hat{I}_{\texttt{R1}}$ and $\hat{I}_{\texttt{N1}}$ present a high variance. Note that, surprisingly, $\hat{I}_{\texttt{N3}}$ has better performance than sampling from the target, i.e., estimator $\bar{I}$. This is because the sampling $\mathcal{S}_{3}$ can be seen as a sampling from the mixture of proposals $\psi({\bf x})$ (which coincides with the target in this example) with a variance reduction technique, as we discuss in Appendix A. Note also that the inequality proved in Theorem 6.1 holds since all methods are unbiased and therefore the MSE is due only to the variance. We can see that $\hat{I}_{\texttt{R2}}$ and $\hat{I}_{\texttt{N2}}$ also behave badly in terms of variance.

2 Applying the MIS schemes in adaptive IS (AIS)

We apply the different MIS schemes within an AIS context. In particular, we focus on the LAIS algorithm, recently proposed in [Martino et al. 2017]. The method consists of an upper layer with a MCMC that draws samples from the target, while a lower layer uses those samples as location parameters (means) of some proposal pdfs for applying IS. In its basic version, $J$ Metropolis-Hastings chains independently run at the upper layer, and hence MIS is applied in the lower layer with $J$ proposals at each iteration. In the following, we implement the six adaptive MIS schemes in a spatial manner for two different target pdfs. For instance, the N3 scheme is implemented by sampling exactly one sample from each of the $J$ proposals at the $t$-th iteration, and applying at the denominator of the IS weight the whole mixture of $J$ proposals as in Eq. (7.7) (see the green square of Fig. 4).

Let us first consider a mixture of five bivariate Gaussians,

2.2 Multidimensional banana-shaped distribution.

We consider the banana shape target example used in [Haario et al. 1999; Haario et al. 2001] which “can be be calibrated to become extremely challenging” [Cornuet et al. 2012]. The target is based on a $d_{x}$-dimensional multivariate Gaussian ${\bf x}\sim\mathcal{N}({\bf x};\textbf{0}_{d_{x}},{\bf\Sigma})$ with $\Sigma=\text{diag}(\sigma^{2},1,...,1)$, where the second variable is nonlinearly transformed from $x_{2}$ to $x_{2}-b(x_{1}^{2}-\sigma^{2})$. This transformation leads to a banana-shaped distribution with zero mean and uncorrelated components (note that the target dimension $d_{x}\geq 2$).

3 Discussion on the experimental results

The numerical experiments confirm that N3 provides the best performance. The scheme R3 also presents a good performance in most cases. The performance of R1 and N1 is, in general, much worse than the performance of the other schemes. Both schemes account at the weight denominator only for the proposal from which the sample is drawn, which in a multimodal scenario can be problematic. While R1 is a novel scheme that has naturally arisen in this work, and it probably has little interest from a practical point of view, N1 has been applied in different adaptive MIS algorithms, such as the original version of PMC [Cappé et al. 2004].

The novel schemes R2 and N2 have appeared in this new framework and deserve a further analysis. The hierarchy theoretically proved for $N=2$ proposals in Theorem 6.2 still holds in the numerical examples for $N>2$, e.g. in Figs. 5(a) and 5(b). In some scenarios, for instance where there is a big number of proposals compared to the modes of the target, these schemes can attain most of the variance reduction of N1 and N3 while reducing the number of proposal evaluations w.r.t. N3. In the example with AIS methods, both R2 and N2 present a very competitive performance w.r.t. to N3.

Finally, observe that in Fig. 5, when a small number of samples $M$ is employed, the schemes N1, N2 and N3, i.e., those with index selection without replacement ($\mathcal{S}_{2}$ and $\mathcal{S}_{3}$), behave better. This occurs because the variance associated to the index selection is reduced by guaranteeing that all proposal pdfs are always used.

Conclusions

In this work, we have introduced a unified framework for sampling and weighting in the context of multiple importance sampling (MIS). This approach extends the concept of a proper weighted sample, enabling the design of a wide range of sampling/weighting combinations. In particular, we have considered three specific sampling procedures and we have proposed five types of generic weighting functions (related to different conditional and marginal distributions which depend on the sampling scheme). As a result of the combinations of sampling and weighting procedures, we have analyzed the six unique resulting schemes (three of them are not present in the literature to the best of our knowledge). We have provided a theoretical comparison of these schemes in terms of variance, establishing a ranking of the different methods in terms of performance and computational complexity. Moreover, we have discussed the application of the MIS schemes within adaptive procedures. In addition, we have provided the practitioner with several useful and easy-to-follow guidelines for applying the MIS schemes in different scenarios. We have analyzed the behavior of the MIS schemes in three different numerical examples which corroborate the previous theoretical analysis.

ACKNOWLEDGMENTS

We thank the Editor, Associate Editor, and referees for their constructive comments that helped to improve the paper. V.E. acknowledges support from the Agence Nationale de la Recherche of France under PISCES project (ANR-17-CE40-0031-01). M.F.B. thanks the support of the National Science Foundation (NSF) under Award CCF-1617986.

References

A Further observations about the sampling 𝒮3\mathcal{S}_{3}

In the sampling procedure $\mathcal{S}_{3}$, ${\bf X}_{n}\sim q_{n}({\bf x})$ for $n=1,\dots,N$, i.e., the selection of the index is deterministic. Note that the set of samples $\{{\bf x}_{n}\}_{n=1}^{N}$ is used in the IS estimators regardless of the order they are drawn. It can be interpreted that the $N$ samples are drawn from the mixture $\psi({\bf x})=\frac{1}{N}\sum_{n=1}^{N}q_{n}({\bf x})$ via Rao-Blackwellization (see [Owen 2013, Section 9.12] for more details). More formally, if we define the r.v. ${\bf X}={\bf X}_{n}\quad\mbox{ with }\quad n\sim\mathcal{U}\{1,2,\ldots,N\}$, then ${\bf X}\sim\psi({\bf x})$. The procedure $\mathcal{S}_{3}$ follows a similar principle as a well-known variance reduction method, known as the stratified sampling [Robert and Casella 2004; Liu 2004], where the domain of ${\bf X}$ is divided into different regions that, in the case of sampling $\mathcal{S}_{3}$, are unbounded and overlapped [Owen 2013, Section 9.12]. Finally, note that the approach $\mathcal{S}_{3}$ can also be seen as the application of a quasi-Monte Carlo technique [Niederreiter 1992] for generating the deterministic sequence of indexes $j_{1}=1,j_{2}=2,\ldots,j_{N}=N$ (uniform, in the sense of low-discrepancy sequence) and then drawing ${\bf x}_{n}\sim q_{j_{n}}({\bf x})=q_{n}({\bf x})$ for $n=1,\ldots,N$. Note also, that $\mathcal{S}_{3}$ can be seen as a residual resampling step of the indexes of the proposals. Since all weights of the proposals are the same, the resampling is fully deterministic, which explains part the variance reduction of the MIS schemes with sampling $\mathcal{S}_{3}$.

B Connections with resampling methods

Resampling methods are used in PFs to replace a set of weighted particles with another set of equally weighted particles. The way we address the sampling process in MIS has clear connections with the resampling step in PFs (e.g., see [Douc and Cappé 2005]). An important difference of the proposed framework is that the MIS proposals are equally weighted in the mixture. The sampling method $\mathcal{S}_{1}$ is then equivalent to the multinomial resampling, whereas the sampling methods $\mathcal{S}_{2}$ and $\mathcal{S}_{3}$ correspond to residual resampling (note that, since $M=N$ and all the proposals are equally weighted, exactly one sample per proposal is drawn). In future works, it would be interesting to analyze sampling schemes related to residual, stratified and systematic resamplings, which can be incorporated quite naturally in MIS schemes, when the weights of the proposals are different (see for instance [He and Owen 2014]).

C Proofs of unbiasedness of the MIS estimators

In this appendix we prove the unbiasedness of the estimator $\hat{I}$ of Eq. (2.4) for the five weighting options described in Section 4. We recall that the general expression for the expectation of $\hat{I}$ within the proposed framework is

Option 1 ($\mathcal{W}_{1}$): $\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{j_{1:n-1}}({\bf x}_{n})=p({\bf x}_{n}|j_{1:n-1})$. We first marginalize in Eq. (C.1) over all indexes that do not affect the $n$-th weight ($j_{n:N}$):

Then, substituting $\varphi_{j_{1:n-1}}({\bf x}_{n})=p({\bf x}_{n}|j_{1:n-1})$ into Eq. (C.2), canceling terms and marginalizing $j_{1:n-1}$, we have:

Option 2 ($\mathcal{W}_{2}$): $\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{j_{n}}({\bf x}_{n})=p({\bf x}_{n}|j_{n})$. We substitute $\varphi_{j_{n}}({\bf x}_{n})=p({\bf x}_{n}|j_{n})$ into Eq. (C.1), which cancels the denominator:

Option 3 ($\mathcal{W}_{3}$): $\varphi_{\mathcal{P}_{n}}({\bf x}_{n})=\varphi_{n}({\bf x}_{n})=p({\bf x}_{n})$. Since $\varphi_{n}$ does not depend on any index, we can first marginalize over the whole set of indexes $j_{1:N}$ in Eq. (C.1):

Then, substituting $\varphi_{n}=p({\bf x}_{n})$ in Eq. (C.3):

Option 4 ($\mathcal{W}_{4}$): $\varphi_{\mathcal{P}_{n}}({\bf x})=\varphi_{j_{1:N}}({\bf x})=f({\bf x}|j_{1:N})=\frac{1}{N}\sum_{n=1}^{N}q_{j_{n}}({\bf x})$. In this case, the expectation of $\hat{I}$ can be expressed as:

Substituting $\varphi_{j_{1:N}}({\bf x})=f({\bf x}|j_{1:N})=\frac{1}{N}\sum_{n=1}^{N}q_{j_{n}}({\bf x})$ in Eq. (C.4), and cancelling the denominator:

Option 5 ($\mathcal{W}_{5}$): $\varphi_{\mathcal{P}_{n}}({\bf x})=\varphi({\bf x})=f({\bf x})=\frac{1}{N}\sum_{n=1}^{N}q_{n}({\bf x})=\psi({\bf x})$. Now, the expectation of $\hat{I}$ becomes

where, in the last step, we have used the identity

for any valid sampling procedure within this framework (see Remark 3.1 and Section 3.5 for more details). Substituting $\varphi({\bf x})=\psi({\bf x})$ in Eq. (C.5)

D Variance analysis of the MIS estimators

that approximates $I$. The variance of $\hat{I}$ can be expressed in the general form as

In the general case of Eq. (D.2), the $N$ terms of the sum of the estimator in ${\hat{I}}$ are dependent. However, in the specific cases where they are independent, the variance of a sum of r.v.’s can be simplified as the sum of the variances, i.e.,

In some MIS schemes, the $N$ terms are dependent (due to a sampling without replacement or because the $n$-th weight depends on several indexes $j_{k}$, with at least one $k\neq n$). However, conditioned to the whole set of indexes $j_{1:N}$, the terms of the sum in Eq. (D.1) are always conditionally independent, so we can apply

In the following, we analyze the variance of the six MIS schemes discussed through this paper under the assumptions described in Theorem 6.1 (see Section 6 for more details). Since some schemes arise under more than one sampling/weighting combination (see Table 6), here we always use the combination that facilitates the analysis.

1. [R1] Sampling 1 / Weighting 2: In this scheme, all the terms of the sum in Eq. (D.1) are independent, so we can use Eq. () for computing the variance of ${\hat{I}}$. Substituting $\varphi_{j_{n}}({\bf x}_{n})=p({\bf x}_{n}|j_{n})=q_{j_{n}}({\bf x}_{n})$ in ,

were we have used that $P(j_{n})=\frac{1}{N}$, $\forall j_{n}\in\{1,...,N\}$.

2. [R2] Sampling 1 / Weighting 4: The expression for the conditional independence of Eq. () is used substituting $\varphi_{j_{1:N}}({\bf x}_{n})=f({\bf x}_{n}|j_{1:N})=\frac{1}{N}\sum_{k=1}^{N}q_{j_{k}}({\bf x}_{n})$ and averaging it over the $N^{N}$ equiprobable sequences of indexes $j_{1:N}$:

where we have used the identity $f({\bf x}|j_{1:N})=\frac{1}{N}\sum_{n=1}^{N}q_{j_{n}}({\bf x}_{n})$. This expression for the variance resembles that of scheme [N3], averaged over the $N^{N}$ possible mixtures (combinations) that can arise with sampling $\mathcal{S}_{1}$.

3. [R3] Sampling 1 / Weighting 3: All the elements are independent in the sum, and the weights do not depend on any index of the set $j_{1:N}$. Therefore, we can start with Eq. (), marginalize over the indexes, and substitute $\varphi_{n}({\bf x}_{n})=p({\bf x}_{n})=\psi({\bf x}_{n})$,

4. [N1] Sampling 3 / Weighting 3: The methods that use sampling without replacement introduce correlation at the selection of the proposals. However, under the perspective of the deterministic sampling ($\mathcal{S}_{3}$), the $n$-th sample ${\bf x}_{n}$ is a realization of the r.v. $X_{n}\sim q_{n}$ and is independent of the other samples. Marginalizing first Eq. () over the indexes, and substituting $\varphi_{n}({\bf x}_{n})=p({\bf x}_{n})=q_{n}({\bf x}_{n})$:

5. [N2] Sampling 2 / Weighting 1: In this scheme, we use again the expression for conditional independence of Eq. (). Substituting $\varphi_{j_{1:n-1}}=p({\bf x}_{n}|j_{1:n-1})$,

Since the the integrals only depend on the set of indexes $j_{1:n}$, each term of the sum has been first marginalized over $j_{n+1:N}$. The first term in the sum can then be further marginalized over $j_{n}$ to obtain the final expression. Note that the variance is the average of the variance of all the $N!$ possible sequences of indexes in the sampling without replacement.

6. [N3] Sampling 3 / Weighting 5: We have followed the same arguments of scheme N1. Marginalizing Eq. () over all the set of indexes $j_{1:N}$, and substituting $\varphi_{n}({\bf x}_{n})=f({\bf x}_{n})=\psi({\bf x}_{n})$:

where we have used the identity $\psi({\bf x})=\frac{1}{N}\sum_{n=1}^{N}q_{n}({\bf x})d{\bf x}$.

D.2 Proof of Theorem 6.1

The proof of Theorem 6.1 is split in the next three propositions.

Proof: See that Eqs. (D.5) and (D.8) are equivalent. ∎

Proof: Subtracting Eqs. (D.7) and (D.8), we get

Now, let us note that the left-hand side of Eq. (D.11) is the inverse of the arithmetic mean of $q_{1}({\bf x}),\ \ldots,\ q_{N}({\bf x})$,

whereas the right hand side of Eq. (D.11) is the inverse of the harmonic mean of $q_{1}({\bf x}),\ \ldots,\ q_{N}({\bf x})$,

Therefore, the inequality in Eq. (D.11) is equivalent to stating that $\frac{1}{A_{N}}\leq\frac{1}{H_{N}}$, or equivalently $A_{N}\geq H_{N}$, which is the well-known arithmetic mean–harmonic mean inequality for positive real numbers [Hardy et al. 1952; Abramowitz and Stegun 1972; Gwanyama 2004]. ∎

Proof: Subtracting (D.7) and (D.10), we get

with $a_{n}=\int\frac{{\pi}({\bf x})g({\bf x})}{\psi({\bf x})}q_{n}({\bf x})d{\bf x}$. The inequality of Eq. (D.12) holds, since it is the definition of the Cauchy-Schwarz inequality [Hardy et al. 1952],

with $b_{n}=1$ for $n=1,...,N$. $\hfill\Box$

Proof of Theorem 6.1. The proof is obtained by applying Propositions D.1, D.2, and D.3. ∎

D.3 Proof of Theorem 6.2

Let us first particularize the variance expression for $N=2$. From Eq. (D.8),

Proof: See that Eqs. (D.17) and (D.18) are equivalent. ∎

Proof: Analyzing Eqs. () and (), we see that Eq. (D.17) can be rewritten as

Proof of Theorem 6.2. The proof is obtained by applying Propositions D.4 and D.5. ∎

We hypothesize that Theorem 6.2 might also hold for $N>2$. The MIS schemes R2 and N2 seem to average estimators with variance reduction (related to N3) with estimators with worse variance (related to N1).

Note that the scheme R3 does not appear in Theorem 6.2. Eq. (D.17) can be rewritten as

D.4 Example with closed-form variances

Let us derive the expressions of the example of Section 6.1 by considering the targeted distribution

We consider $N=2$ proposal densities, $q_{1}({\bf x})=\mathcal{N}({\bf x}|-\mu,\sigma^{2})$ and $q_{2}({\bf x})=\mathcal{N}({\bf x}|\mu,\sigma^{2})$. Note that the mixture of proposals is exactly the targeted distribution, i.e. $\psi({\bf x})=\pi({\bf x})$. We address the case where we want to estimate a specific moment $g$ of $\pi$ with the $M=2$ samples. In the following, we provide explicit variances of the unnormalized estimator of Eq. (2.4) for the six MIS schemes. From Eq. (D.5),

Moreover, from Proposition D.4, $\hat{I}_{\texttt{N2}}=\hat{I}_{\texttt{R2}}$.

E Multidimensional mixture of generalized Gaussian distributions

Let us consider a mixture of multivariate generalized Gaussian distributions (GGD) as a target pdf. In particular

where ${\bm{\mu}}_{k}=[\mu_{k,1},...,\mu_{k,d_{x}}]^{\top}$, ${\bm{\alpha}}_{k}=[\alpha_{k,1},...,\alpha_{k,d_{x}}]^{\top}$, and ${\bm{\beta}}_{k}=[\beta_{k,1},...,\beta_{k,d_{x}}]^{\top}$ are respectively the mean, scale, and shape parameters of each component of the mixture. Each component of the mixture factorizes in all dimensions, i.e., the multivariate GGD pdf is the product of $N$ unidimensional GGD pdfs. Namely,

where $\kappa_{k,d}=\frac{\beta_{k,d}}{2\alpha_{k,d}\Gamma\left(\frac{1}{\beta_{k,d}}\right)}$, $\Gamma(\cdot)$ is the gamma function, and $x_{d}$ is the $d$-th dimension of ${\bf x}$. This family of distributions includes both Gaussian and Laplace distributions with $\beta=2$ and $\beta=1$, respectively. In this example, $\mu_{1,d}=-3$, $\mu_{2,d}=1$, $\mu_{3,d}=5$, $\beta_{1,d}=1.1$, $\beta_{2,d}=1.8$, $\beta_{3,d}=5$, $\alpha_{1,d}=\alpha_{2,d}=\alpha_{3,d}=1$ for all $d=1,...,d_{x}$. The expected value of the target ${\pi}({\bf x})$ is then $E_{\pi}[{X_{d}}]=1$ for $d=1,...,d_{x}$. In order to study the performance of the different MIS schemes, we vary the dimension of the state space in Eq. (E.1) testing different values of $d_{x}$ (with $2\leq d_{x}\leq 10$). We consider the problem of approximating via Monte Carlo the expected value of the target density, and we compare the performance of all MIS schemes. In this example, we use $N=500$ non-standardized t-student densities as proposal functions, where each location parameter has been selected uniformly within the $^{d_{x}}$ square, and the scale parameters and the degree of freedom parameters have been selected as $\sigma_{n,d}=5$ and $\nu_{n,d}=5$, respectively, for $n=1,...,N$ and $d=1,...,d_{x}$. For each method, we draw $M=kN$ samples, with $k=32$, and we average all the results over $200$ runs.

Fig. 7 shows the MSE in the estimation of the mean of the target (averaged over all dimensions) when we increase the dimension $d_{x}$. Note that the hierarchy established in Section 6 also holds in this example regardless the dimension. In this case, methods R1 and N1 behave poorly even at lower dimensions, while the other MIS schemes have a similar behavior. When we increase the dimension, all the methods degrade, and, at certain point ($d_{x}\geq 6$), the performance of all of them is similar. Note that the proposal pdfs are fixed in random locations of the space, which is well covered at low dimensions (since we are using $N=500$ pdfs), but this coverage becomes worse as the dimension increases. This can probably explain the similar performance of all the methods in higher dimensions.