Estimation of cosmological parameters using adaptive importance sampling
We present a Bayesian sampling algorithm called adaptive importance sampling or Population Monte Carlo (PMC), whose computational workload is easily parallelizable and thus has the potential to considerably reduce the wall-clock time required for sampling, along with providing other benefits. To assess the performance of the approach for cosmological problems, we use simulated and actual data consisting of CMB anisotropies, supernovae of type Ia, and weak cosmological lensing, and provide a comparison of results to those obtained using state-of-the-art Markov Chain Monte Carlo (MCMC). For both types of data sets, we find comparable parameter estimates for PMC and MCMC, with the advantage of a significantly lower computational time for PMC. In the case of WMAP5 data, for example, the wall-clock time reduces from several days for MCMC to a few hours using PMC on a cluster of processors. Other benefits of the PMC approach, along with potential difficulties in using the approach, are analysed and discussed.
pacs:98.80.es, 02.50.-r, 02.50.Sk
In recent years we have seen spectacular advances in observational cosmology, with the availability of more and more high quality data allowing for the testing of models with higher complexity. Some of these tests have been made possible thanks to the use of Bayesian sampling techniques, and in particular Markov Chain Monte Carlo (MCMC) – an (iterative) algorithm that produces a Markov chain whose distribution converges to the target posterior . After a “burn-in” period, samples from such a chain can be regarded as samples approximately from . Proposed values for the chain or the updating scheme of MCMC can be designed to ensure that moves towards regions of higher mass under are favored, and regions with null probability (under ) are never visited. This way, most of the computational effort can be spent in the region of importance to the posterior distribution, and an MCMC approach is usually much more efficient than traditional grid sampling of the model parameter space.
The MCMC technique is now well known in cosmology, and in particular in its most simple form, the Metropolis-Hastings algorithm, thanks to the user-friendly and freely available package COSMOMC (Lewis and Bridle, 2002). Other forms of the MCMC algorithm, like Gibbs-sampling and Hybrid Monte Carlo (better known in cosmology as Hamiltonian sampling), have also been proposed and have found some interesting usage in the estimation of the posterior distribution for the Cosmic Microwave Background anisotropy power spectrum at low resolution (see (Rudjord et al., 2008) and references therein, also (Taylor et al., 2008) and (Dunkley et al., 2009)).
For all its advantages over grid sampling, the MCMC approach also
suffers from problems. One difficulty is to assess the correct
convergence of the chain. Another lies in the presence of correlations
within the chain which can greatly reduce the efficiency of the
sample (Robert and Casella, 2004). A third issue which is
particularly relevant for the usage of MCMC in cosmology is the
computational time involved. Indeed, whatever the sampling technique,
we often need to compute at least one estimate of the posterior for
each sampled point. This computation can be slow in cosmology. With
the current processing speed of computers, a point of the posterior
of, for example, the WMAP5 data set, using CAMB
(Lewis et al., 2000)
On the hardware side, however, there is a route to speed improvement that does not lie in quicker CPUs, but on the availability of cheap multi-CPU computers and the standardization of clusters of computers. This opportunity, however, is only partly opened to MCMC. Indeed, there are two ways of parallelizing the parameter exploration. First, by distributing the computation of the likelihood, which is not always possible and does not always lead to speed improvement. Second, by running multiple chains in parallel. This last option is the simplest, but is ‘forbidden’ by the iterative nature of the MCMC algorithm. More precisely, running parallel chains and mixing them in the end to build a bigger chain sample is of course possible (and can be advantageous in fully exploring the support of ), but at the condition that each of the individual chains has converged. In the absence of such a condition, significant biases in the sample can be introduced. Determining convergence for each chain is inherently difficult in practice and has largely prevented more widespread use of the approach (Rosenthal, 2000). Thus, for MCMC any speed improvements through parallelization are difficult to achieve.
In this paper, we propose another sampling algorithm suitable for cosmological applications, that is not based upon MCMC, and can be parallelised. This novel algorithm, called Population Monte Carlo (PMC) is an adaptive importance sampling technique, that has been studied recently in the statistics literature (Cappé et al., 2007). While this algorithm solves some of the issues of MCMC in cosmology, the approach of course has a different set of potential problems that we will analyse and discuss, along with its advantages.
The paper is outlined as follows. In the next section, we provide a brief introduction to the Bayesian approach, which we hope will give the casual reader some important keys for further readings, and we also discuss the challenges and issues involved with using either an MCMC or an importance sampling algorithm for estimation. We then describe details of the PMC approach. In Sect. III, we assess the performance of the PMC approach using a simulated target density with features similar to cosmological parameter posteriors, and provide a comparison to results obtained using an MCMC approach. In Sect. IV, we illustrate the results from the PMC approach using actual data, consisting of CMB anisotropies, supernovae of type Ia and weak cosmological lensing. We conclude in Sect. V with a discussion and an outlook for further work.
ii.1 Bayesian inference via simulation
A key feature of Bayesian inference is to provide a probabilistic expression for the uncertainty regarding a parameter of interest by combining prior information along with information brought by the data. Prior information, for example, could take the form of information obtained from previous experiments which cannot readily be incorporated into the current experiment or simply consist of a feasible range. The absence of prior information, however, is not restriction for the use of Bayesian inference and estimation can still be regarded as valid (Robert, 2001). Information brought by the data and prior information are entirely subsummed in the posterior probability density function obtained, up to a normalization constant, by
It is however generally difficult to handle the posterior distribution, due to (a) the dimension of the parameter vector , and (b) the use of non-analytical likelihood functions. For both of these reasons, the normalizing constant missing from the right-hand side of (1) is usually not explicitly available. A practical solution to this difficulty is to replace the analytical study of the posterior distribution with a simulation from this distribution, since producing a sample from allows for a straightforward approximation of all integrals related with , due to the Monte Carlo principle (Robert and Casella, 2004). In short, if is a sample drawn from the distribution and denotes a function (with finite expectation under ), the empirical average
is a convergent estimator of the integral
in the sense that the empirical mean (2) converges to as grows to infinity. Quantities of interest in a Bayesian analysis typically include the posterior mean, for which ; the posterior covariance matrix corresponding to ; and probability intervals, with , where is a domain of interest, and denotes the indicator function which is equal to one if and zero otherwise.
ii.2 Markov chain Monte Carlo methods
For most problems in practice, direct simulation from is not an option and more sophisticated approximation techniques are necessary. One of the standard approaches (Robert and Casella, 2004) to the simulation of complex distributions is the class of Markov chain Monte Carlo (MCMC) methods that rely on the production of a Markov chain having the target posterior distribution as limiting distribution.
MCMC can be implemented with many Markovian proposal distributions but the standard approach is the random walk Metropolis-Hastings algorithm: given the current value of the chain, a new value is drawn from , where the so-called proposal denotes a symmetric probability density function. The point is then accepted as with probability (also called acceptance rate in this context)
and otherwise, . The algorithm is implemented as follows:
Random walk Metropolis-Hastings algorithm
Choose an arbitrary value of .
Generate and .
While this algorithm is universal in that it applies to any choice of posterior distribution and proposal , its performance highly depends on the choice of the proposal that has to be properly tuned to match some characteristics of . If the scale of the proposal is too small, that is, if it takes many steps of the random walk to explore the support of , the algorithm will require many iterations to converge and, in the most extreme cases, will fail to converge in the sense that it will miss some relevant part of the support of (Marin et al., 2005). If, on the other hand, the scale of is too large, the algorithm may also fail to adequately sample from . This time, the chain may exhibit low acceptance rates and fail to generate a sufficiently diverse sample, even with longer runs. There exist monitors that assess the convergence of such algorithms but they usually are conservative – i.e., require a multiple of the number of necessary iterations – and partial – i.e., only focus on a particular aspect of convergence or on a special class of targets (Robert and Casella, 2004). MCMC algorithms are also notoriously delicate to calibrate on-line, both from a theoretical point of view and from a practical perspective (Haario et al., 2001). For these approaches, often called adaptive MCMC, some recommendations for the optimal scaling and calibration schedule for various proposals in high dimensions have been proposed (Roberts and Rosenthal, 2001), but this is still at an experimental stage.
ii.3 Population Monte Carlo
Population Monte Carlo (PMC) (Cappé et al., 2004, 2007) is an adaptive version of importance sampling (Von Neumann, 1951; Rubinstein, 1981) that produces a sequence of samples (or populations) that are used in a sequential manner to construct improved importance functions and improved estimations of the quantities of interest.
We recall that importance sampling is based on the fundamental identity (Robert and Casella, 2004)
which holds for any probability density function with support including the support of and any function for which the expectation is finite. Hence, this approach to approximating integrals linked with complex distributions is also universal in that the above identity always holds. If are drawn independently from ,
provides a converging approximation to . In this context, is called the importance function and are commonly referred to as importance weights. For Bayesian inference, one cannot directly use (6) as only the unnormalised version of (i.e., the right-hand side of eq. 1) is available. Conveniently, the self-normalised importance ratio
where the normalised importance weights are defined as
is also a converging approximation to , independent of the normalization of . For an importance function that is closely matched to the target density, significant reductions in the variance of the Monte Carlo estimates are possible in comparison to estimates obtained using MCMC (Robert and Casella, 2004). However, the importance sampling approach is equally prone to poor performances as MCMC, in that the resulting converging approximation may suffer from a large or even infinite variance if is not selected in accordance with . There is no universal importance function and most of the research in this field aims at fitting the most efficient importance functions for the problem at hand.
Population Monte Carlo offers a possible solution to this difficulty through adaptivity: given the target posterior density up to a constant, PMC produces a sequence of importance functions aimed at approximating this very target. The first sample is produced by a regular importance sampling scheme, , associated with importance weights
and their normalised counterparts (eq. 8), providing a first approximation to a sample from . Moments of can then be approximated to construct an updated importance function , etc.
The approximation can be measured in terms of the Kullback divergence (also called Kullback-Leibler divergence or relative entropy) from the target,
and the density can be adjusted incrementally such that is smaller than . The importance function should be selected from a family of functions which is sufficiently large to allow for a close match with but for which the minimization of (10) is computationally feasible. In Cappé et al. (2007) the authors propose to use mixture densities of the form
where is a vector of adaptable weights for the mixture components (with and ), and ) is a vector of parameters which specify the components; is a parameterised probability density function, usually taken to be multivariate Gaussian or Student-t (where the latter is to be preferred in cases where it is suspected that the tails of the posterior are indeed heavier than Gaussian tails). Given the vast array of densities that can be approximated by mixtures, such an importance function provides considerable flexibility to efficiently estimate a wide range of posteriors, including in this case those found in cosmological settings.
The generic PMC algorithm then consists of the following:
Population Monte Carlo algorithm
Choose an importance function .
Generate an independent sample .
Compute the importance weights .
Update the importance function to , based on the previous weighted sample .
Generate independently .
Compute the importance weights .
Unlike for MCMC, in a PMC approach, the process can be interrupted at any time as the sample produced at each iteration can be validly used to approximate expectations under using self-normalised importance sampling following (7). Further, sampling outputs from previous iterations can be combined (Veach and Guibas, 1995; Cornuet et al., 2009), and the sample size at each iteration does not necessarily need to be fixed. Both of these properties of PMC can be exploited to improve parameter estimates, either by increasing the coverage of the importance function to the target density or increasing the precision of the approximation for the integral of interest.
Also note that an approximate sample from the target density can be obtained by sampling () with replacement, using the normalised importance weights . Although this process induces extra Monte Carlo variation, there are a number of methods available which considerably reduce the variation involved (e.g. residual sampling (Liu and Chen, 1995a) or systematic sampling (Whitley, 1994)).
Updating the importance function in the Gaussian case
In this section, we particularise the generic PMC algorithm to the case where the importance function consists of a mixture of -dimensional Gaussian densities with mean and covariance ,
Using this importance function for the mixture model (11), we start the PMC algorithm by arbitrarily fixing the mixture parameters , and then sample independently from the resulting importance function to obtain our initial sample . From this stage, updates of the parameters proceed recursively.
At iteration , the importance weights associated with the sample are given by
with normalised counterparts given by eq. (8). The parameters ( and ) of the importance function are then updated according to
Appendix A provides derivations of these expressions and further details on the general approach, as well as equations pertaining to the (more involved) case of mixtures of multivariate Student-t distributions, which are used in the simulations presented in Sect. III.
As discussed in Appendix A, the main theoretical appeal of this particular update rule is that, as tends to infinity, the corresponding Kullback divergence is guaranteed to be less than .
The above update process can be repeated a number of times, and although there is no need for a formal stopping rule, some measures of performance against the target density can be used as a guide. As the objective of importance function adaptations is to minimise the Kullback divergence between the target density and the importance function, we can stop the process when further adaptations do not result in significant improvements in . To this end, it can be shown that may be estimated by the normalised perplexity
is the Shannon entropy of the normalised weights, a frequently used measure of the quality of an importance sample. Thus, minimization of the Kullback divergence can be approximately connected with the maximization of the perplexity (18). Values of this criterion close to 1 will therefore indicate good agreement between the importance function and the target density.
Another frequently used criterion for importance sampling is the so-called effective sample size (ESS),
which lies in the interval and can be interpreted as the number of sample points with non-zero weight (Liu and Chen, 1995b). Both measures (18, 20) are interconnected, as an importance function which is close to the target density will have both a high normalised perplexity and a relatively large number of points with non-zero weight, compared to an ill-fitting importance function. Given a real-valued function of interest one can also estimate the asymptotic variance of the self-normalised importance sampling estimator (cf. eq. 7) using the importance sample itself, as
Beware that this formula (which is derived from theorem 2 of Geweke (1989)) is only valid with normalised weights, and that it is a variance conditional on the current importance proposal , i.e. it does not take into account the adaptation. This measure can be related to the so-called integrated autocorrelation time (IAT) used for Markov chain Monte Carlo simulations, which, in this case, takes into consideration the level of autocorrelation present in the chain (Robert and Casella, 2004; Roberts and Rosenthal, 2001; Dunkley et al., 2005).
A first illustration
To illustrate the PMC approach, we explore a banana-like target density presented Fig. 1. The same target distribution will be studied in greater depth in the next section. The results of the first eleven (11) iterations of the PMC algorithm using a mixture of Student-t densities are shown Fig. 2 (see also Appendix A for details of the update procedure).
While this target density shows slightly more pronounced curvature for an example of a posterior density in practice, it serves to illustrate the process of adaptation of the importance function. The initial importance function is a mixture of multivariate Student-t’s, consisting of nine components placed around the centre of the range for the first two variables, each with a relatively large variance (for the first dimension , second ) and degrees of freedom . The different coloured circles in Fig. 1 indicate the location of the component means, and the circle size is proportional to the weight associated with the component. At the fourth iteration (), we see that the importance function starts to resemble the shape of the target density, with components becoming more separated and moving into the tails of the target. By the sixth iteration () the importance function has further adapted to the shape of the target banana density. For this target density and importance function, Fig. 3 presents estimates of the normalised perplexity and normalised effective sample size (ESS/) for the first 10 iterations over 500 simulation runs. As shown, the estimates of the normalised perplexity improve rapidly from approximately 0.14 for the second importance function (Iteration 2) to approximately 0.81 for the last importance function (Iteration 10), with a similar increase in estimates of the normalised effective sample size (ESS/, increasing from 0.10 to 0.60). For this importance function and target density, the normalised perplexity starts to level off after the 10th iteration (around ), indicating that there is no need for further adaptation of the proposal density. As mentioned previously however, in general, one does not need to observe the convergence of the proposal (as for MCMC) in order to stop the sampling process.
An important consideration, and a choice that needs to be made at the start of the algorithm are the parameter values , and for the initial importance function , including the degrees of freedom in the case of the Student-t mixture, and the sample size . A poor initial importance function, such as one that is tightly centred around only one mode in the case of a multimodal posterior or a narrow importance function with light tails, may take a long time to adapt or may miss important parts of the posterior. For importance sampling the choice of requires both fat tails and a reasonable match between and the target in regions of high density. Such an importance function can be more easily constructed in the presence of a well informed guess about the parameters and possibly the shape of the posterior density. Sample size considerations also play an important role – smaller samples can adapt quite quickly with less computational time but may provide less reliable information about the posterior density relative to larger samples. Such considerations are important as we look at posterior densities of increasingly high dimensions, and thus we can expect to take a larger sample size as the dimension of the problem increases. We will discuss these issues further in the context of simulated and actual data, and also in Sect. V.
In this section, we test the performance of PMC using simulated data, and compare the results to an adaptive MCMC procedure.
In order to provide a good test for both approaches we use the target density considered in (Haario et al., 2001), which is difficult to explore but which also provides a realistic scenario for many problems encountered in cosmology. The target density is based on a centred -multivariate normal, with covariance , which is slightly twisted by changing the second co-ordinate to . Other co-ordinates are unchanged. We obtain a twisted density which is centered with uncorrelated components. Since the Jacobian of twist is equal to 1, the target density is:
For the target density that we will consider, we set , , , which results in a banana-shaped density in the first two dimensions (see Fig. 1).
For the target density considered, interest is in how well PMC and MCMC are able to approximate the tails of the target. Whilst the curvature present in the first two dimensions of this target density is slightly more pronounced than what is typically seen in practice for cosmology it serves to highlight the difficulties faced by both PMC and MCMC in covering the parameter space. In particular, little accurate information is available in order to guide the choice of importance function (for PMC) and proposal distribution (for MCMC) and so both approaches are forced to learn about the parameter space.
Test run proposal for PMC
For PMC, and in the absence of any detailed a priori information about the target density, except the possible range for each variable, we have chosen the first importance function to be a mixture of multivariate Student-t distributions with components displaced randomly in different directions slightly away from the centre of the range for each variable: the mean of the components is drawn from a -multivariate Gaussian with mean and covariance equal to where is some positive-definite matrix; the variance for components was chosen to be . We choose a mixture of components of Student-t distribution with degrees of freedom; and is a diagonal matrix with diagonal entries . This choice of ensures adequate coverage, albeit somewhat overdispersed, of the feasible parameter region. In this simulated example, Student-t distributions are preferable to Gaussian distributions because the range of the variables is unbounded (in contrast to the cosmology examples to be discussed in Sect. IV).
A representation of the first importance function for the first two dimensions is shown in the top left-hand box in Fig. 2, with a typical evolution over the next few iterations in the other panels. In pilot runs of various importance functions against the target density, the best fitting importance function required at least seven components in order to adequately represent the coverage of the entire density.
For PMC, an important issue is the sample size for each iteration. A poor initial importance function with a relatively small sample size will take a long time to adapt or it may even be unable to recover sufficiently to provide reasonable parameter estimates. Such problems are more likely to occur as the dimension of the parameter space increases, the so-called curse of dimensionality. For the simulation exercise each iteration is based on a sample of k points. To prevent numerical instabilities in the updating of the parameters, components with a very small weight () or containing less than sample points are discarded.
Test run proposal for MCMC
As little information is available for the target density, an adaptive MCMC approach is used which can allow for faster learning of the target density than using either independent or non-adaptive random walk proposals Roberts et al. (1997). For MCMC, the proposal distribution is a centred Gaussian with covariance matrix which is updated along the iterations. An important choice for adaptive MCMC concerns the scaling of the proposal and the rate of adaptation. There has been much research on this Roberts et al. (1997); Atchadé and Fort (2008), and a common choice for the covariance of the Gaussian proposal is to consider where is an estimate of the covariance of the target density, at update . The choice is considered to be optimal when the chain is in its stationary regime, and for target densities that have a product form Roberts et al. (1997). However for the target density we consider this does not hold: the first two components are not independent despite being uncorrelated and dependence is not linear but quadratic. However, with no other theoretical results to follow we start with a scaling factor of that form and for the simulation results to follow assess the effect on convergence and results using alternative values. We update the covariance matrix by the recursive formula
where is the sample covariance of the previous update, and is the covariance of the sampled estimates from the previous update to the current iteration. The value of is with chosen suitably to allow for a cooling of the update, which is a necessary condition to ensure convergence of this adaptive MCMC to the target density as well as convergence of the empirical averages Roberts and Rosenthal (2007); Atchadé and Fort (2008).
In pilot runs, we explored the effect of this schedule for various values of in and we observe that the choice of plays a role on the time to convergence (for the estimation of the quantities of interest, see below) and on the acceptance rate of the chain.
To ensure a fair comparison with PMC, we start the chain at a random point drawn from the same Gaussian distribution as for PMC (i.e , using the same values for as used for PMC). We also explored in pilot runs the role of the initial value of the chain: despite it being known that MCMC is sensitive to the choice of the initial position of the chain - which has no real counterpart in PMC - this hasn’t been found to have a major impact on performance (for reasonable choice of the initial value at least) in this particular study. We also fixed the update schedule to be every k points and we assessed the effect on the results from using less or more points before updating. For the simulation results to follow, has been set to which ensures convergence after the burn-in period (see Sect. III.1), and a mean acceptance rate at convergence of about . The proposal distribution is updated for the entire length of the chain and is not stopped after the burn-in period.
iii.1 Test runs
For PMC and the proposal outlined, the perplexity appeared to level off at around the th iteration, so for the results to follow for PMC we ran the PMC algorithm for iterations (k points per iteration) and used a final draw of k points. To assess MCMC for the same number of points we used a chain length of k points with a burn-in of k points. Results for both approaches at successive intervals before k points are also provided. To assess the performance of the approaches, each simulation was replicated times.
iii.2 Results of the simulations
For the results of the simulations, we are interested in both the mean estimates of the parameters (in particular for and ) and also estimates of the confidence region coverage ( and ) which will provide an indication of how well both approaches are covering the tails of the target density. For each run , we provide the results for various functions of interest:
We note here as the indicator function for the region. and are indicators only for the first dimension, while and are dealing with the first 2. In all cases, the remaining dimensions are marginalised over.
Table 1 shows the results for estimation of for functions and ( and respectively). The results provided show the mean and standard deviation of estimates calculated over 500 runs. Although the performance is quite similar for both methods, PMC does display a two-fold reduction in standard deviation compared to MCMC for both functions. A closer look at the results reveals that for the empirical distributions of the estimates (see Figure 4) are quite similar for both methods, except for the variance which is much reduced for PMC. For , on the other hand, the empirical distribution of the estimates for PMC are quite skewed, resulting in a slight positive bias for the majority of the runs (second panel of Figure 4). The difference between and can be explained from Figure 1 which shows that failure to visit sufficiently the downward low-probability tails does indeed imply a positively biased estimates for the mean of the second component. PMC does appear to be more sensitive to this issue than MCMC, despite the fact that the estimates for MCMC display a larger overall variability.
Figure 5 provides the results for the confidence region coverage. To depict the variability of the data, the results are displayed by using whisker plots. The results from both PMC and MCMC against all of the performance measures are similar, with both showing good coverage of the target density. The distribution of this estimator is again more skewed for PMC than it is for MCMC, particularly for the 95% regions in the bottom panel of Figure 5. Nevertheless, the variability of the estimates obtained with PMC also is significantly reduced compared to MCMC.
Figure 4 shows the evolution of the results for and from k points to k points for both PMC (left panels) and MCMC (right panels). The results from Fig. 4, in general, highlight the reduction in variance of the Monte Carlo estimates for PMC in comparison to MCMC. In particular, it is interesting to note that the variance of the estimates, either for or for 100k posterior evaluations under MCMC are comparable to estimates obtained using PMC at the second iteration (20k points).
Simulating from this target distribution is a challenging problem for both methods. In particular, the use of a vague initial importance function in a multidimensional space represents a challenge to PMC and it has been observed that the importance function takes some time to properly adapt to the target density (about 10 iterations). The choice of the initial importance function in PMC is more crucial than is the choice of the initial proposal distribution in adaptive MCMC. Although different variations for updating the covariance matrix for the MCMC approach are possible we did not see a significant improvement in the results presented from using alternative covariance structures. For most of the simulation results, the proposal covariance matrix was observed to adapt relatively quickly to the true covariance matrix. Changes to the covariance structure considered included changes to the update frequency, the starting proposal , the scaling of the proposal (value of ) and adaptation of the covariance (value of ). Hence, the PMC approach may require more precise a priori knowledge of the target density than MCMC.
In the next section, we apply the PMC approach to typical cosmological examples, and provide results in comparison to MCMC.
Iv Application to cosmology
We apply our new adaptive importance sampling method to the posterior of cosmological parameters. Flat CDM models with either a cosmological constant (CDM) or a constant dark-energy equation-of-state parameter (CDM) are explored and tested with recent observational data of CMB anisotropies, supernovae type Ia and cosmic shear, as described in the next section.
The three data sets and likelihood functions used here are the same as in Kilbinger et al. (2008); the CMB measurements and likelihood are based on the five-year WMAP data release Hinshaw et al. (2009), the SNIa data set is the first-year SNLS survey Astier et al. (2006), while the cosmic shear is from the CFHTLS-Wide third release (T0003, Fu et al., 2008). The results presented in the following sections can be compared to the MCMC analysis in Kilbinger et al. (2008).
iv.1 Data sets
To obtain theoretical predictions of the CMB temperature and polarization power- and cross-spectra we use the publicly available package CAMB (Lewis et al., 2000). The likelihood is calculated using the public WMAP5 code (Dunkley et al., 2009).
The WMAP5 likelihood takes as input the TT,TE,EE and BB theoretical power spectra calculated by CAMB, and returns a likelihood computed from a sum of different parts. It computes a pixel-based Gaussian likelihood based on template-cleaned maps (Gold et al., 2008) and their associated inverse covariance matrices (see Page et al. (2007) for details) at large angular scales ( for TT, for TE, EE and BB). At small angular scales, it computes an approximate likelihood based on pseudo-spectra and their associated covariance for TT and TE (Hinshaw et al., 2007), based respectively on the (Q,V) and (V,W) channel pairs for TE and TT.
In addition, the likelihood computation takes into account analytic marginalisations on nuisance parameters such as the beam transfer function and point-sources uncertainties (Hinshaw et al., 2007; Nolta et al., 2008). We ignore corrections due to SZ and impose a larger (flat) prior on the Hubble constant. Indeed, CMB data alone exhibits a degeneracy between the Hubble constant and e.g. the cosmological constant (Efstathiou and Bond, 1999) which is removed by adding other cosmological probes.
The acoustic oscillation peaks in the CMB anisotropy spectrum are a standard ruler at a redshift of about . CMB therefore measures the angular diameter distance to that redshift which depends mainly on the total matter-energy density () and weakly on the Hubble constant . The overall anisotropy amplitude is determined by the large-scale normalization . The relative height of the peaks is sensitive to the baryonic and dark matter densities. On large scales, secondary anisotropies are generated at late times () due to reionisation, which is parametrised by the optical depth , and the integrated Sachs-Wolfe (ISW) effect, which is a probe of .
The SNLS data set is described in detail in Astier et al. (2006). We use their results from the SNIa light-curve fits which, for each supernova, provides the rest-frame -band magnitude , the shape or stretch parameter and the colour . We use the standard likelihood analysis described in Kilbinger et al. (2008), adopted from Astier et al. (2006).
Under the assumption that supernovae of type Ia are standard candles we can fit the luminosity distance to the SNIa data. The luminosity distance is a function of , and . Three additional parameters are the universal SNIa magnitude and the linear response parameters to stretch and colour, and , respectively. Those three parameters are specific to our choice of distance estimator, and can be regarded as nuisance parameters. The Hubble constant is integrated into the parameter , so there is no explicit dependence on in the SNIa posterior.
The CFHTLS-Wide year data release (T0003), the data and weak lensing analysis as well as cosmological results are described in Fu et al. (2008). As in Fu et al. (2008) we use the aperture-mass dispersion between 2 and 230 arc minutes as second-order lensing observable (Schneider et al., 1998). We assume a multivariate Gaussian likelihood function and take into account the correlation between angular scales. The theoretical aperture-mass dispersion is obtained by non-linear models of the large-scale structure Smith et al. (2003).
The galaxy redshift distribution is obtained by using the CFHTLS-Deep redshift distribution Ilbert et al. (2006) and rescaling it according to the magnitude distribution of CFHTLS-Wide galaxies. We fit the resulting histogram with eq. (14) from Fu et al. (2008), introducing the three fit parameters . The histogram data is modeled as multivariate, uncorrelated Gaussian, the corresponding likelihood is included, independent of the lensing likelihood, in the analysis.
Weak gravitational lensing by the large-scale structure is sensitive to the angular diameter distance and the amount of structure in the Universe. It is an important probe to measure the normalisation on small scales. With the current data, this parameter is however largely degenerate with . This degeneracy is likely to be lifted by future surveys which will include the measurement of higher-order statistics Takada and Jain (2004); Kilbinger and Schneider (2005) and shear tomography Hu (1999). In particular from the latter a great improvement on the determination of is to be expected, a parameter which is only weakly constrained by lensing up to now Hoekstra et al. (2006); Jarvis et al. (2006).
iv.2 Cosmological parameter and priors
We sample a hypercube in parameter space which corresponds to flat priors for all parameters, see Table 2 for more details. Additional priors exist, both in explicit and implicit form, which represent regions of parameter space which are unphysical or where numerical fitting formulae break down. For example, we exclude extremely high baryon fractions () because of numerical problems in the computation of the transfer function. Further, for very low values of both and the pivot scale for the non-linear power spectrum is outside the allowed range. Very rarely, the calculation of the likelihood for individual points in parameter space is unsuccessful because of numerical errors or limitations of the likelihood code. Since these points cannot be taken into account, a pragmatic solution is to formally modify the prior to exclude those points. Note that these rare cases occur mainly in regions of very low likelihood.
|Total matter density||0.01||1.2||C||S||L|
|Dark-energy eq. of state||-3.0||0.5||C||S||L|
|Primordial spectral index||0.7||1.4||C||L|
|Normalization (large scales)||C|
|Normalization (small scales)
|Absolute SNIa magnitude||S|
|galaxy -distribution fit||L|
iv.3 Initial choice of the importance function
As described earlier in Sect. II.3.3, it is important to have a good guess for the initial importance function. In all cases considered here, we rely on an estimate of the maximum likelihood point and the Hessian at that point (Fisher matrix) to build our initial proposals. We use the conjugate-gradient approach Press et al. (1992) to find the maximum likelihood point at which to calculate the Fisher matrix using the theoretical model. We construct a mixture model consisting of Gaussian components. Student-t mixtures with small degrees of freedom were tested and turned out to be a bad approximation to the posterior under study, resulting in low perplexities. Each mixture component is shifted randomly from the maximum by a small amount. A random scaling is applied to the covariance of each component; i.e. the eigenvectors and ratios between the eigenvalues of the covariance are the same as the ones of the Fisher matrix.
We obtain good results for shifts of about 0.5% to 2% of the box size. Here, a trade-off between too large shifts (resulting in low importance weights) and too small shifts (components stay near the maximum, the posterior tails do not get sampled) has to be found. The stretch factor is chosen randomly between typical values of 1 and 2. In some cases, in particular with high dimensionality, the derivation of the Fisher matrix is not stable and the matrix is numerically singular. In such cases we set the off-diagonal elements of to zero.
We found a sample size between 7 500 and 10 000 points to be adequate for most cases. The number of components of the initial importance function was chosen between 5 to 10. For the final iteration we used a sample size five times that of the initial sample size.
The PMC algorithm is reliable and very efficient in sampling and exploring the parameter space. Both the perplexity as well as the effective sample size increase quickly with each iteration (Fig. 6). The perplexity reaches values of 0.95 or more in many cases, although in particular for higher dimensional posteriors the final values are lower. Satisfactory results (i.e. yielding consistent mean and marginals compared to MCMC, see below) are obtained for perplexities larger than about 0.6.
The distribution of importance weights gets narrower from iteration to iteration (Fig. 7). Initially, many sampled points exhibit very low weights. After a few iterations, the importance function has moved towards the posterior increasing the efficiency of the sampling.
Our initial mixture model starts with all mixture components close to the maximum likelihood point. With consecutive iterations the components spread out to better cover the region where the posterior is significant. This can be seen in Fig. 8.
Compared to traditional MCMC, our new PMC method is faster by orders of magnitude. The time-consuming calculation of the posterior can be performed in parallel and therefore a speed-up by a factor of the number of CPUs is obtained. In times where clusters of multi-core processors are readily available, this speed-up is easily of the order of 100. In addition, MCMC has a low efficiency with typical acceptance rates of 0.25 – 0.3. The PMC normalised effective sample size in the WMAP5 case is 0.7 which results in a much larger final sample for the same number of posterior calculations of around .
We emphasise again that with MCMC one can make only limited use of parallel computing since one has to wait for each Markov chain to converge, and because it is not straightforward to combine chains, as mentioned earlier.
Comparison with MCMC
The MCMC results we present here are either obtained using the adaptive MCMC algorithm or a classical one. Indeed, as we show in the following, adaptive MCMC can have some issues that a less efficient classical MCMC algorithm can avoid. Apart from those special cases, the MCMC and adaptive MCMC gave very similar results, the latter usually reaching a better acceptance rate, and thus a better efficiency.
We find excellent agreement between using our respective implementations of MCMC (adaptive or not) and PMC. Mean, confidence intervals and 2d-marginals are very similar using both methods. The performance of PMC is superior to MCMC in some cases, which is illustrated by the following examples.
An inherent problem of MCMC is that even for a long run there can be regions in parameter space that are not sampled in an unbiased way. This is illustrated in Fig. 10. The feature at the 99.7%-level of MCMC (left panel, for large values of and ) is due to an “excursion” of the chain into a low-likelihood region at step 130 500, lasting for 300 steps. We ran the chain for 300 000 steps and the feature was still visible. A second run of the chain did not exhibit this anomaly. This kind of sample “noise” can be prevented by running a chain for a very long time or by combining several (converged) chains. Such features are much less likely to occur in an importance sample which consists of uncorrelated points.
A second issue are parameters which are nearly unconstrained by the data with the result that the marginalised posterior in that dimension is flat. To illustrate this we choose weak lensing alone which can not constrain (Fig. 11). Using the Fisher matrix as initial Gaussian proposal for adaptive MCMC, the chain stays in a small region in the -direction; the covariance being very flat, most jumps ends up out of the prior distribution. This results in an update variance for this parameter which is much too small, and in a bad exploration of the posterior in this flat direction as shown Fig. 11. The classical MCMC algorithm, with the same proposal yields better results, but with a very low acceptance rate and needing 500 000 steps to reach the result presented Fig. 11. Alternatively, modifying the initial proposal to be smaller and better adapted to the prior, or increasing the covariance stretch factor from the optimal value of (see Sect. III.0.3) to helps the chain to explore more of the parameter space in the latter steps of the adaptation. These modifications to the algorithm also result in very low acceptance rate, and somehow go against the very idea of an adaptive algorithm, since they require very fine tuning of the initial proposal!
With PMC we obtain a much better performance and recover very well the flat posterior.
In Tables 3 and 4 we show the mean and 68% confidence intervals for CMB alone for the CDM model and for lensing+SNIa+CMB using CDM, respectively. The differences in mean and 68%-confidence intervals is less than a few percent in most cases. Figure 10 shows that the lower-confidence regions and the correlation between parameters agrees very well between (non-adaptive) MCMC and PMC.
In this paper, we have introduced and assessed an adaptive importance sampling approach, called Population Monte Carlo (PMC), which aims to overcome the main difficulty in using importance sampling, namely the reliance on a single efficient importance function. PMC achieves this goal by iteratively adapting the importance function towards the target density of interest. A significant appeal of the approach, when compared to alternatives such as MCMC, lies in the possibility to use (massive) parallel sampling which considerably reduces the computational time involved in the estimation of parameters for many astrophysical and cosmological problems. Simulated and actual data have been used in this work to assess the performance of PMC for estimation of parameters in a Bayesian inference with features approaching classical cosmological parameter posteriors.
The PMC approach is, in essence, an iterated importance sampling scheme that simultaneously produces, at each iteration, a sample approximately simulated from a target distribution and an approximation of in the form of the current proposal distribution. As such, the samples produced by the PMC approach can be exploited as regular importance sampling outputs at any iteration . Samples from previous iterations can be combined (Cornuet et al., 2009), and approximations like can be updated dynamically, without necessarily requiring the storage of samples.
Although adaptation of the importance function has the explicit aim of improving the coverage of the posterior density there are instances where this objective may not be met. In some cases, successive updates of the importance function may result in: (a) an importance function which is too peaked and which has light tails (invalid importance function); (b) an importance function which fits only one mode (in the case of a multimodal posterior); (c) numerical problems due to the adaptation procedure (usually involving poor conditioning of some of the covariance matrices). Such cases are likely to produce a poor approximation to the integral of interest, or alternatively lead to highly variable parameter estimates over iterations. These problems can be quickly discovered or signalled by observing a poor ESS, and parameter estimates or normalised perplexity which do not stabilise after a few iterations.
Such cases of poor performance as outlined can be significantly reduced by choosing a reasonably well informed initial importance function with a large enough sample size at each iteration, especially on the initial iteration that requires many points to counter-weight a potentially poor importance function. In general, the initial importance function should be chosen to cover a region of the parameter space that has support larger than the posterior. In the absence of reliable prior information, finding such an importance function may be difficult to do. One approach may be to locate the components in the centre of the feasible range (if available) for each variable, with reasonably large variances to ensure some coverage of the parameter space. We found this approach to be reasonably successful for the simulated data case discussed in Sect. III. In the presence of some prior information, for example an estimate of the maximum-likelihood point and an approximation of the covariance matrix (via the Hessian), components can be placed around these points with variance comparable to the approximation. Another approach may be to perform a singular value decomposition of the covariance matrix, and make use of the eigenvectors and eigenvalues to place components along the most likely directions of interest. Alternatively and in the same spirit, components can be placed according to the principal points of the resulting sample, using a k-means clustering approach (Tarpey, 2007). Both approaches have been reasonably successful for a range of posterior densities examined, and by placing components in regions of high posterior support in addition to the mode have the potential to further significantly reduce the number of iterations for difficult posterior densities.
The main appeals or advantages of the PMC method are worth re-emphasising at this point:
Parallelisation of the posterior calculations
Low variance of Monte Carlo estimates
Simple diagnostics of ‘convergence’ (perplexity)
We address these three points in more detail now.
(1.) The first advantage, namely the ability to parallelise the computational
task, is becoming increasingly useful through the availability of
cheap multi-CPU computers and the standardization of clusters of
computers. Software to implement the parallelisation task, such as
Message Parsing Interface (MPI)
(2.) In general, for PMC and an importance function that is closely matched to the target density, significant reductions in the variance of the Monte Carlo estimates are possible in comparison to estimates obtained using MCMC (Robert and Casella, 2004). For example, for the posterior estimates for the WMAP5 data we observed a 10-fold reduction in variance for the same number of sample points as observed for MCMC. Such reductions suggest that the computational time savings extend not only to the number of CPU’s available but to smaller sample requirements for PMC in total compared to MCMC to achieve similar variability of estimates. For cosmological applications, this observation is valuable as we observed, e.g., in Fig. 6 for CMB data, that the fit between the adapted importance function and the target posterior is sometimes quite good. By combining samples across iterations further computational savings are also possible. The absence of construction of a Markov chain for PMC can also have the desirable attribute of reducing sample noise, as observed for the SNIa data in Sect. IV.4.2.
(3.) As shown in Sect. II.3.2, the perplexity (eq. 18) is a relatively simple measure of sampling adequacy to the target density of interest. For MCMC and other approaches which rely on formal measures of convergence, assessment of convergence can be very difficult with users facing a potential array of associated diagnostic tests.
In addition to the above points, a further appeal of PMC is the ability to provide a very good approximation to the marginal posterior or evidence, which naturally follows as a byproduct of the approach. To demonstrate this appeal, further research is underway to explore the use of PMC in the context of model selection problems in cosmology.
Acknowledgements.We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is provided by the NASA Office of Space Science. We thank the Planck group at IAP and the Terapix group for support and computational facilities. DW and MK are supported by the CNRS ANR “ECOSSTAT”, contract number ANR-05-BLAN-0283-04 ANR ECOSSTAT. The authors would like to thank F. Bouchet, S. Bridle, J.-M. Marin, Y. Mellier and I. Tereno for helpful discussions.
Appendix A Details of the importance function updates for PMC
The method proposed in Cappé et al. (2007) to adaptively update the parameters of the importance function is based on a variant of the EM (Expectation-Maximization) algorithm Dempster et al. (1977), which is the standard tool for the estimation of the parameters of mixture densities. We describe below the principle underlying the algorithm of Cappé et al. (2007), showing in particular that each iteration decreases, up to the importance sampling approximation errors, the Kullback divergence between the target and the importance function .
Remember that our goal is to minimise (10), as increases, by iteratively tuning the parameters of the mixture importance function defined in (11). Developing the logarithm in (10), this objective can be equivalently formulated in terms of the maximization of the following quantity
with respect to . Using Bayes’ rule, we denote by
the posterior probability that belongs to the th component of the mixture (for the mixture parameters ). The EM principle consists of evaluating, at iteration , the following intermediate quantity
Using the concavity of the log as well as the expression of in (25), it is easily checked that
and hence that . Thus, any value of which increases the intermediate quantity above the level also results in, at least, an equivalent increase of the actual objective function . In the EM algorithm, one sets and to the values where the intermediate quantity is maximal, thus satisfying the previous requirement. Furthermore, the maximization of leads to a closed form solution whenever belongs to a so-called exponential family of probability densities.
In the example of the multivariate Gaussian density recalled in (12), the parameter consists of the mean and the covariance matrix and the intermediate quantity may be written as
up to terms that do not depend on , or . Routine calculations show that the maximum of (28) is achieved for
In practice, both the numerator and denominator of each of the above expressions are integrals under which must be approximated. The solution proposed in Cappé et al. (2007) is based on self-normalised importance sampling using the weighted sample simulated at the previous iteration . The corresponding empirical update equations are given in eqs. (14)–(16) of Sect. II.3.1.
The Student-t distribution provides a family of multivariate densities with parameters and which have the same interpretation as in the Gaussian case (except for the fact that the covariance is equal to rater than ) but with an additional shape factor which allow for heavier tails: letting yields back the Gaussian but for , one obtains a density with polynomially decreasing tails whose only finite moments are the two first ones (note that it is also possible to extend the family to the case where ). Using mixtures of Student-t distributions will thus be mostly useful in cases where the target posterior distribution itself has heavy tails. The parameter update corresponding to mixtures of Student-t distributions is a bit more involved but follows the same general pattern. For the sake of completeness, we just recall below the formulas given in Cappé et al. (2007):
and denotes the -dimensional Student-t probability density function,
Sampling from a multivariate Student-t distribution is most easily undertaken by using its derivation in terms of a multivariate Gaussian () and chi-squared distribution (),
and taking advantage of the fact that sampling from Y and Z is straightforward.
- For WMAP5, is a deduced quantity that depends on the other parameters
- A. Lewis and S. Bridle, Physical Review D 66, 103511 (2002).
- Ø. Rudjord, N. E. Groeneboom, H. K. Eriksen, G. Huey, K. M. Górski, and J. B. Jewell, submitted to ApJ (2008), also arXiv:0809.4624.
- J. F. Taylor, M. A. J. Ashdown, and M. P. Hobson, MNRAS 389, 1284 (2008).
- J. Dunkley, E. Komatsu, M. R. Nolta, D. N. Spergel, D. Larson, G. Hinshaw, L. Page, C. L. Bennett, B. Gold, N. Jarosik, et al., ApJS 180, 306 (2009).
- C. Robert and G. Casella, Monte Carlo Statistical Methods (Springer-Verlag, New York, 2004), 2nd ed.
- A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000).
- W. A. Fendt and B. D. Wandelt, submitted to ApJ (2007), also arXiv:0712.0194.
- T. Auld, M. Bridges, M. P. Hobson, and S. F. Gull, MNRAS 376, L11 (2007).
- A. Hajian, Physical Review D 75, 83525 (2007).
- S. Geman and D. Geman, IEEE Transactions on Pattern Analysis and Machine Intelligence 6, 721 (1984).
- D. L. Larson, H. K. Eriksen, B. D. Wandelt, K. M. Gorski, G. Huey, J. B. Jewell, and I. J. O’Dwyer, The Astrophysical Journal 656, 653Ã±660 (2007).
- R. Shaw, M. Bridges, and M. P. Hobson, Mon. Not. R. Astron. Soc. 378, 1365 (2007).
- F. Feroz and M. P. Hobson, submitted to MNRAS (2008).
- N. Chopin and C. P. Robert, arxiv.org/abs/0801.3887 (2008).
- J. S. Rosenthal, Far East J. Theor. Stat. 4, 207 (2000).
- O. Cappé, R. Douc, A. Guillin, J.-M. Marin, and C. Robert, Statist. Comput. in press (2007), also arXiv:0710.4242.
- C. Robert, The Bayesian Choice (Springer-Verlag, New York, 2001), 2nd ed.
- J.-M. Marin, K. Mengersen, and C. Robert, in Handbook of Statistics, edited by C. Rao and D. Dey (Springer-Verlag, New York, 2005), vol. 25.
- H. Haario, E. Saksman, and J. Tamminen, Bernoulli 7, 223 (2001).
- G. O. Roberts and J. S. Rosenthal, Statistical Science 16, 351 (2001).
- O. Cappé, A. Guillin, J.-M. Marin, and C. Robert, J. Comput. Graph. Statist. 13, 907 (2004).
- J. Von Neumann, J. Resources of the National Bureau of Standards – Applied Mathematics Series 12, 36 (1951).
- R. Rubinstein, Simulation and the Monte Carlo Method (John Wiley, New York, 1981).
- E. Veach and L. J. Guibas, in SIGGRAPH 95 Proceedings (Addison-Wesley, 1995), pp. 419–428.
- J.-M. Cornuet, J.-M. Marin, A. Mira, and C. P. Robert, in prep (2009).
- J. Liu and R. Chen, Journal of the American Statistical Association 90, 567 (1995a).
- D. Whitley, Statistics and Computing 4, 65 (1994).
- J. Liu and R. Chen, J. American Statist. Assoc. 90, 567 (1995b).
- J. Geweke, Econometrica 57, 1317 (1989).
- J. Dunkley, M. Bucher, P. G. Ferreira, K. Moodley, and C. Skordis, MNRAS 356, 925 (2005).
- G. O. Roberts, A. Gelman, and W. R. Gilks, Annals of Applied Probability 7, 110 (1997).
- Y. Atchadé and G. Fort, Tech. Rep. (2008), arXiv math.PR/0807.2952.
- G. O. Roberts and J. S. Rosenthal, Journal of Applied Probability 44, 458 (2007).
- M. Kilbinger, K. Benabed, J. Guy, P. Astier, I. Tereno, L. Fu, D. Wraith, J. Coupon, Y. Mellier, C. Balland, et al., Submitted to A&A (2008), also arXiv:0810.5129.
- G. Hinshaw, J. L. Weiland, R. S. Hill, N. Odegard, D. Larson, C. L. Bennett, J. Dunkley, B. Gold, M. R. Greason, N. Jarosik, et al., ApJS 180, 225 (2009), eprint 0803.0732.
- P. Astier, J. Guy, N. Regnault, R. Pain, E. Aubourg, D. Balam, S. Basa, R. G. Carlberg, S. Fabbro, D. Fouchez, et al., A&A 447, 31 (2006).
- L. Fu, E. Semboloni, H. Hoekstra, M. Kilbinger, L. van Waerbeke, I. Tereno, Y. Mellier, C. Heymans, J. Coupon, K. Benabed, et al., A&A 479, 9 (2008).
- B. Gold, C. L. Bennett, R. S. Hill, G. Hinshaw, N. Odegard, L. Page, D. N. Spergel, J. L. Weiland, J. Dunkley, M. Halpern, et al., ApJS in press (2008), also arXiv:0803.0715.
- L. Page, G. Hinshaw, E. Komatsu, M. R. Nolta, D. N. Spergel, C. L. Bennett, C. Barnes, R. Bean, O. Doré, J. Dunkley, et al., ApJS 170, 335 (2007).
- G. Hinshaw, M. R. Nolta, C. L. Bennett, R. Bean, O. Doré, M. R. Greason, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, et al., ApJS 170, 288 (2007).
- M. R. Nolta, J. Dunkley, R. S. Hill, G. Hinshaw, E. Komatsu, D. Larson, L. Page, D. N. Spergel, C. L. Bennett, B. Gold, et al., ApJS in press (2008), also arXiv:0803.0593.
- G. Efstathiou and J. R. Bond, MNRAS 304, 75 (1999).
- P. Schneider, L. van Waerbeke, B. Jain, and G. Kruse, MNRAS 296, 873 (1998).
- R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, MNRAS 341, 1311 (2003).
- O. Ilbert, S. Arnouts, H. J. McCracken, M. Bolzonella, E. Bertin, O. Le Fèvre, Y. Mellier, G. Zamorani, R. Pellò, A. Iovino, et al., A&A 457, 841 (2006).
- M. Takada and B. Jain, MNRAS 348, 897 (2004).
- M. Kilbinger and P. Schneider, A&A 442, 69 (2005).
- W. Hu, ApJ 522, L21 (1999).
- H. Hoekstra, Y. Mellier, L. van Waerbeke, E. Semboloni, L. Fu, M. Hudson, L. Parker, I. Tereno, and K. Benabed, ApJ 457, 116 (2006).
- M. Jarvis, B. Jain, G. Bernstein, and D. Dolney, ApJ 644, 71 (2006).
- W. H. Press, S. A. Teukolsky, B. P. Flannery, and W. T. Vetterling, Numerical Recipes in C (Cambridge University Press, 1992).
- T. Tarpey, Computational Statistics 22, 71 (2007).
- A. Dempster, N. Laird, and D. Rubin, J. Royal Statist. Society Series B 39, 1 (1977).