Reconstructing Probability Distributions with Gaussian Processes

TL;DR

This paper introduces a Gaussian process-based method to reconstruct and resample probability distributions from existing MCMC samples, significantly reducing computational costs for cosmological data analysis.

Contribution

The authors present a novel Gaussian process regression approach to reconstruct log-probability distributions from MCMC samples, enabling efficient resampling and analysis.

Findings

Reconstructed the Planck 2018 posterior distribution.

Generated a new MCMC chain with forty times more points in thirty minutes.

Demonstrated the method's utility for importance sampling and convergence testing.

Abstract

Modern cosmological analyses constrain physical parameters using Markov Chain Monte Carlo (MCMC) or similar sampling techniques. Oftentimes, these techniques are computationally expensive to run and require up to thousands of CPU hours to complete. Here we present a method for reconstructing the log-probability distributions of completed experiments from an existing MCMC chain (or any set of posterior samples). The reconstruction is performed using Gaussian process regression for interpolating the log-probability. This allows for easy resampling, importance sampling, marginalization, testing different samplers, investigating chain convergence, and other operations. As an example use-case, we reconstruct the posterior distribution of the most recent Planck 2018 analysis. We then resample the posterior, and generate a new MCMC chain with forty times as many points in only thirty minutes. Our likelihood reconstruction tool can be found online at https://github.com/tmcclintock/AReconstructionTool.

Equations

$$⟨y^{\ast }|{\theta }^{\ast },\{{\mathbf{\theta }}_i,y_i{\}}_{i=1}^N⟩=\mu ({\mathbf{\theta }}^{\ast })+{\Sigma }_{\ast }^T{\Sigma }^{-1}\mathbf{y}.$$

$$k({\mathbf{\theta }}_i,{\mathbf{\theta }}_j)=\exp \left[-\frac{({\mathbf{\theta }}_i-{\mathbf{\theta }}_j{)}^2}{2l^2}\right].$$

$$\ln \mathcal{L}(\mathbf{y}|l,\mathbf{\theta })=-\frac{1}{2}\left[(\mathbf{y}-\mu {)}^T{\Sigma }^{-1}(l,\mathbf{\theta })(\mathbf{y}-\mu )+\log \det (\Sigma (l,\mathbf{\theta }))+d\log 2\pi \right].$$

$${\hat{\bar{p}}}_i$$

$$\hat{\mathbf{C}}(p_i,p_j)$$

$$\hat{\mathbf{C}}=\mathbf{Q\Lambda }{\mathbf{Q}}^T,$$

$$\mathbf{x}=s{\mathbf{\Lambda }}^{\frac{1}{2}}{\mathbf{Q}}^T{\mathbf{x}}^{\prime }+\widehat{\bar{\mathbf{p}}},$$

$$d^2=(\mathbf{x}-\mathbf{p}{)}^T{\hat{\mathbf{C}}}^{-1}(\mathbf{x}-\mathbf{p}).$$

$${\mathbf{p}}_{\mathbf{x}}=\underset{\mathbf{p}}{\mathrm{argmin}}{d}^2(\mathbf{x},\mathbf{p})|\forall \mathbf{p}\in \{{\mathbf{p}}_i{\}}_{i=1}^{N_c}.$$

$${\mathbf{q}}^{\prime }=\mathbf{Qp}.$$

$${\mathbf{q}}_i=\frac{{\mathbf{q}}_i^{\prime }-{\hat{\mathbf{q}}}^{\prime }}{{\hat{\mathbf{\sigma }}}^{\prime }}\forall {\mathbf{q}}_i^{\prime }\in \{{\mathbf{q}}_i^{\prime }{\}}_{i=1}^{N_c},$$

Full text

††pubyear: 2019††pagerange: Reconstructing Probability Distributions with Gaussian Processes–Reconstructing Probability Distributions with Gaussian Processes

1 Introduction

Cosmology has entered the “precision era” where experiments aim to put ever-tighter constraints on a small number of physical parameters that describe our Universe. Since the discovery of the accelerated expansion of the universe using Type Ia supernovae (Riess et al., 1998; Perlmutter et al., 1999) up through recent analyses of the large scale structure of the universe (e.g. Beutler et al., 2011; Alam et al., 2017; Abbott et al., 2018a; Planck Collaboration et al., 2018; Hikage et al., 2019; Hildebrandt et al., 2018), the primary means of making such constraints has stayed the same. Using Markov Chain Monte Carlo (MCMC) or similar techniques, most analyses end up producing samples of the posterior probability distribution of the parameters of interest given the data under some model. These samples are used to compute statistical moments of the marginalized distributions of the parameters, which are reported as final results.

The plethora of cosmic probes means that one can combine experiments in order to enhance overall constraining power (Abbott et al., 2018b; Baxter et al., 2019, for recent examples). Unfortunately, the technical and detailed nature of individual analyses can make it difficult, if not impossible, to make certain combinations. Indeed, even reproducing past results can be laborious. So much so, that a cottage industry has sprouted just to create software for performing joint analyses and enabling reproducibility (e.g. Cosmosis111https://bitbucket.org/joezuntz/cosmosis/wiki/Home Zuntz et al. 2015, Cobaya222https://github.com/CobayaSampler/cobaya, CosmoMC333https://cosmologist.info/cosmomc/ Lewis & Bridle 2002, and MontePython444https://github.com/brinckmann/montepython_public Audren et al. 2013). As indicated by the breadth of detailed instructions for installing and configuring these packages, resampling a posterior from a completed analysis is not as simple as pushing a button. Compounding this problem is the shear computational power required. For example, an MCMC chain from the DES Year 1 key project combining cosmic shear, galaxy-galaxy lensing, and galaxy clustering took approximately three thousand total CPU hours. Combining DES with an analysis that has a different computational bottleneck would exacerbate the problem, potentially making a full joint analysis prohibitively expensive. A final hurdle to resampling is the simple fact that analysis code is often not released until long after chains have been made public. For instance, the posterior chains from the Planck Collaboration et al. (2018) analysis have been released, but software for evaluating their likelihood has not. Nevertheless, performing joint analyses is essential for discovering new physics. Thus, we seek a solution to the related problems of reproducibility, expense, and inconvenience.

In an effort to tackle these problems simultaneously, we propose using Gaussian process regression to reconstruct the posterior probability distribution of an experiment. Gaussian processes are probabilistic interpolators that can be used to model complicated functions in many dimensions. They have been successfully used in the past in cosmology for emulation (Lawrence et al., 2017; Nishimichi et al., 2018; Zhai et al., 2018; McClintock et al., 2019), and have many applications in machine learning (Rasmussen & Williams, 2005). Our method allows for reconstructing and resampling probability distributions given a small set of samples taken from an MCMC chain. Sampling from the reconstructed distribution can be orders of magnitude faster than evaluating the likelihood of a given experiment. For instance, large scale structure likelihoods are bottlenecked by Boltzmann codes, which do not need to be evaluated using our technique.

As a demonstration, we take the MCMC chains released in Planck Collaboration et al. (2018) and reconstruct the (unreleased) posterior probability of that analysis. By sampling from our reconstructed distribution, we recover constraints on all free parameters. Our MCMC chain has forty times as many samples even after accounting for burn-in, and can be produced in minutes on a laptop computer.

The paper is structured as follows: section 2 covers the basics of Gaussian process regression. In section 3 we discuss how to select training points from an MCMC chain and how we construct the Gaussian process for interpolation. section 4 shows some examples of our method, including reconstructing the Planck Collaboration et al. (2018) posterior. Finally, section 5 concludes with a discussion of some applications of our technique as well as potential extensions.

2 Regression with Gaussian Processes

Here we discuss the basics of regression with Gaussian processes. Gaussian processes have many uses and a thorough discussion of the topic is provided in Rasmussen & Williams (2005). In this work, we use the Gaussian process implementation found in the Python package george555http://george.readthedocs.io/en/latest (Ambikasaran et al., 2015).

In our case, we are have $N$ samples of the log-probability of some distribution $\ln\mathcal{P}$ at locations in parameter space given by $\theta$. In this section, we will write $\ln\mathcal{P}=y$, so that our training data takes the form of a series of $N$ distinct points in parameter space along with the corresponding log-probability, i.e. $\{\bm{\theta}_{i},y_{i}\}_{i=1}^{N}$. We know the values of $y$ exactly, since they are taken directly from an MCMC chain or an equivalent set of samples as the result of an experiment.

We want to build an interpolator in order to be able to predict the log-probability of a new location in parameter space $y^{*}(\theta^{*})$ with no assumptions about the specific function form of the log-probability. Gaussian processes are a tool that is ideally suited for this task. Gaussian processes assume that the covariance between two samples $y_{i}(\bm{\theta}_{i})$ and $y_{j}(\bm{\theta}_{j})$ depends only on the values of $\bm{\theta}_{i}$ and $\bm{\theta}_{j}$, and that together $\{\bm{\theta},y\}$ form a multivariate normal distribution. Our samples $\{\bm{\theta}_{i},y_{i}\}_{i=1}^{N}$ represent $N$ draws from this distribution. Therefore, a prediction for $y^{*}(\bm{\theta}^{*})$ is the mean of the conditional multivariate normal given our samples according to

[TABLE]

In Equation 1, $\mu(\bm{\theta}^{*})$ is the so-called mean function, and $\mathbf{y}$ is a column vector containing $\{y_{i}\}_{i=1}^{N}$. $\Sigma_{*}^{T}$ denotes the $1\times N$ row vector of the covariance evaluated at all pairs of $\bm{\theta}^{*}$ with $\{\bm{\theta}_{i}\}_{i=1}^{N}$, while $\Sigma$ is the $N\times N$ matrix containing the covariance between all pairs of $\{\bm{\theta}_{i}\}_{i=1}^{N}$. Conventionally, the mean function is set to 0 assuming that the training data has been standardized such that the sample mean is $\frac{1}{N}\sum y_{i}=0$, although it need not be. This mean function can also be understood as the estimate for $y^{*}$ when $\bm{\theta}^{*}$ is far away from the known samples, or the value that the Gaussian process extrapolates to. Since we are building an interpolator for log-probabilities, we require this mean function to be small. We set it to a value less than the minimum log-probability in our samples, $\mu<\min(\{y_{i}\})$, and assume that our samples represent the “region of interest” that the original MCMC chain derives from. That is to say, if the MCMC chain does not well sample the parameter space, one cannot expect the interpolator to do any better. We found that the exact value of $\mu$ does not affect our results, as long as its absolute value is not too large compared to the training points. In other words, if $|\mu|\gg\max(\{|y_{i}|\})$ the interpolator will not perform well.

The elements of the covariance matrices are given by the kernel function $[\Sigma]_{ij}=k(\bm{\theta}_{i},\bm{\theta}_{j})$. Specifying this kernel function completes the Gaussian process. We found that the squared exponential kernel function is sufficient. This kernel is given by

[TABLE]

Equation 2 states that the covariance between two points in parameter space follows a Gaussian, so that as points become farther separated they influence each other less. This range of influence is controlled by the hyperparameter $l$, which is the characteristic length scale of correlations within the data. The value of $l$ is found by maximizing the likelihood

[TABLE]

The left term is the $\chi^{2}$ of the data fit by the model while the middle term penalizes excessively large covariance. The right most term is a normalization constant that depends on the dimensionality $d$ of the data (i.e. $\bm{\theta}$ is a vector in a $d$-dimensional parameter space). Upon maximizing Equation 3, $l$ is fixed and we can predict new values of the log-probability according to Equation 1.

Notably, we use only a single hyperparameter in Equation 2. Usually, when using a Gaussian process to interpolate in many dimensions, the characteristic length scale in each of those dimensions can be different, and one requires a different $l$ for each dimension. For instance, the characteristic length scale for interpolating over $\Omega_{m}$ should be very different than the length scale required to interpolate over $A_{s}$, since these parameters differ by many orders of magnitude. In the following section we detail how we transform the training samples in our model so that we require only a single hyperparameter.

3 Modeling Probability Distributions with Gaussian Process Regression

MCMC chains represent probability distributions via the number density of points from the chain at any given location in parameter space. However, samplers also record the log-probability of each entry in the chain, although this quantity is oftentimes not used. In contrast, we reconstruct the log-probability as a function of the parameters in the chain while only using a special selection of samples. In section 3.1, we detail our method of selecting training points from an MCMC chain. These training points are then used to optimize a Gaussian process. Next, in section 3.2 we detail how we construct the interpolator and perform regression.

3.1 Selecting Training Samples

Each sample in an MCMC chain is a vector in a $d$-dimensional parameter space $\mathbf{p}=\{p_{i}\}_{i=1}^{d}$. The sample mean and sample covariance matrix of the parameters can be computed from the chain according to

[TABLE]

where there are $N_{c}$ samples in the chain indexed by $k$, $\hat{\bar{p}}_{i}$ is the sample mean of the $i$-th parameter, and $\hat{\mathbf{C}}$ is a $d\times d$ dimensional matrix. We eigendecompose $\hat{\mathbf{C}}$ into components

[TABLE]

where $\mathbf{\Lambda}$ is a diagonal matrix containing the eigenvalues, and $\mathbf{Q}$ is an orthogonal matrix with columns containing the eigenvectors. $\mathbf{Q}$ can be understood as a rotation matrix that rotates any sample in parameter space into a space in which the parameters are all orthogonal.

Next, we select $N_{s}$ samples in a $d$-dimensional unit hypercube, $\mathbf{x}^{\prime}=\{x_{i}^{\prime}\}_{i=1}^{d}$ centered on the origin. Our results are not strongly sensitive to how these samples are selected. By default, we use a simple Latin-Hypercube sampling as implemented in the package pyDOE666https://github.com/clicumu/pyDOE2(Surowiec et al., 2017). It is unlikely that this is an optimal choice, and selecting training samples for Gaussian processes remains an open topic of research (Heitmann et al., 2009; Nishimichi et al., 2018). Since our results are not sensitive to this choice of sampling, we defer further study of this topic to a future work.

The $N_{s}$ samples are then rotated into parameter space, scaled by a constant factor $s$, and mean centered according to

[TABLE]

where $\hat{\bar{\mathbf{p}}}$ is a vector containing the sample mean of the parameters in the chain, and $\mathbf{\Lambda}^{\frac{1}{2}}$ is a diagonal matrix containing the square roots of the eigenvalues. While $\mathbf{x}^{\prime}$ was a vector in a unit hypercube, $\mathbf{x}$ is a vector in parameter space. The constant factor $s$ “spreads out” the random samples. Without it, the samples are contained approximately within the $1\sigma$ level of the chain. Setting $s\sim 8$ results in points reaching the $\sim 4\sigma$ level. Figure 1 show an example of a selection of training points from a two parameter MCMC chain that sampled a multivariate Gaussian likelihood. The black points are spread evenly across the sample space, while the blue points of the chain are clustered as one might expect.

For a given $\mathbf{x}\in\{\mathbf{x}_{i}\}_{i=1}^{N_{s}}$ the distance between that point and a sample in the chain is given by

[TABLE]

Thus, to select our final set of training points $\mathbf{p}_{\mathbf{x}}$ we find the nearest point in the chain for each $\mathbf{x}$ according to

[TABLE]

This set of points $\{\mathbf{p}_{\mathbf{x},i}\}_{i=1}^{N_{s}}$ is a subset of the original MCMC chain that follows the Latin-Hypercube design we chose. There is a chance that multiple $\mathbf{x}$ select the same sample $\mathbf{p}$, however this can be checked and corrected for by sampling additional training points.

This method for selecting training points is well suited for subsampling points from a smooth, nearly Gaussian distribution. If the points in an MCMC chain have a complicated distribution, one may need to tailor the selection of training points to achieve a desired result. Additionally, one may choose specific points from the chain to include in the subsample, such as the maximum a-posteriori probability point (i.e. the most likely point in the chain), or points in the tails.

3.2 Performing Regression

Constructing the Gaussian process and performing regression requires additional work. Given a new location in parameter space $\mathbf{p}^{*}$, simply applying the inverse transformation of Equation 7 is insufficient. This is because in this basis, the characteristic length in each dimension is different, meaning we require one hyperparameter for each free parameter in the MCMC chain. This becomes infeasible for large parameter spaces such as the 27 dimensional Planck analysis.

Instead, we construct the Gaussian process in a special basis. First, we transform each point in the original MCMC chain according to

[TABLE]

We then compute the sample mean $\hat{\bar{\mathbf{q}}}^{\prime}$ and sample covariance $\hat{\mathbf{C}}_{\mathbf{q}^{\prime}}$ for these points. In this basis, the new samples $\{\mathbf{q}_{i}^{\prime}\}_{i=1}^{N_{c}}$ have a diagonal covariance matrix, however each dimension of $\mathbf{q}^{\prime}$ has a different variance. Thus, our final transformation is to standardize these points by subtracting the sample mean and dividing each dimension by the sample standard deviation

[TABLE]

where $\hat{\bm{\sigma}}^{\prime}$ is a vector containing the sample standard deviations for each dimension. In this parameter space, the set of points $\{\mathbf{q}_{i}\}_{i=1}^{N_{c}}$ has zero mean, unit variance, and no correlation between the parameters in each dimension.

Each of the training points $\{\mathbf{p}_{\mathbf{x},i}\}_{i=1}^{N_{s}}$ selected in section 3.1 correspond to individual $\mathbf{q}$ points in this space. We select this set $\{\mathbf{q}_{i}\}_{i=1}^{N_{s}}$ in order to train the Gaussian process. Since each dimension of parameter space in this basis has unit variance, the Gaussian process requires only a single hyperparameter to successfully model the distribution.

When predicting the log-probability at a new point in parameter space $\mathbf{p}^{*}$, the process is straightforward. First, rotate the parameter vector by multiplying by $\mathbf{Q}$. Second, subtract the sample mean $\hat{\bar{\mathbf{q}}}^{\prime}$ and divide by the sample standard deviation for each dimension $\hat{\bm{\sigma}}^{\prime}$. Third, estimate the log-probability according to Equation 1. This completes the reconstruction of the probability distribution sampled by the input MCMC chain. We can now use the Gaussian process to resample the distribution in order to perform calculations including marginalization, importance sampling, joint sampling, and computing Bayesian evidence. Select examples are presented in the following section.

4 Example Use Cases

The ability to quickly resample a probability distribution has many potential applications. In this section we highlight two examples that have immediate utility.

4.1 Sampling Disjoint Parameter Spaces

Oftentimes two experiments share parameters of interest but contain their own set of nuisance parameters. Using our method, we can reconstruct the individual probability distributions of two or more experiments in order to jointly sample all parameters at once. Figure 2 demonstrates this using a simple example involving three parameters. Suppose we have one experiment that samples $\{x,y\}$ and another that samples $\{y,z\}$, and MCMC samples for both (blue and green contours in the figure, respectively). We reconstruct the probability distributions of both experiments. Assuming that the two experiments are independent and have flat priors, we can simply multiply the two reconstructed probabilities to determine the joint likelihood of the two experiments, leading to the red contours in Figure 2. As an aside, we note that if either experiment had informative priors, we would need to remove one set of priors in order to construct a valid posterior. This demonstrates the utility of our method for quickly combining any set of experiments just using their MCMC chains.

4.2 Resampling the Planck 2018 Posterior

As a cosmologically relevant example, we illustrate how our code can be utilized to reconstruct the Planck Collaboration et al. (2018) TT, TE, EE+lowE posterior.

Figure 3 shows constraints on the first seven parameters of the Planck chain in blue. The entire distribution has 27 free parameters, however all parameters past the first six are nuisance parameters. We checked that our reconstruction successfully recovered those 20 nuisance parameters, but present only the cosmologically interesting parameters here.

Note that the blue contours are less smooth than those found in Planck Collaboration et al. (2018). This is because we have not applied a kernel density estimate for smoothing to highlight the ability of our reconstruction to achieve good convergence. The green contours come from sampling the reconstructed probability distribution using our method. This new MCMC chain is forty times longer than the original Planck chain after removing burn-in.

Resampling our reconstructed posteior yields posteriors that are in excellent agreement with the original chains. Visually inspecting the two distributions appears to show differences in the peaks of the distribution, but this is merely a result of shot noise in the original chains. We verify that our reconstruction accurately recovers the log-posterior probability of the experiment $\ln\mathcal{P}_{\rm Planck}$ as shown in Figure 4. This figure shows a histogram of the fractional error in $\ln\mathcal{P}_{\rm Planck}$ predicted by our method for all points in the original MCMC chain that were not used in the reconstruction. The vast majority of points are predicted with better than $0.2\%$ accuracy. Assuming that $\ln\mathcal{P}_{\rm Planck}\propto\chi^{2}$, this demonstrates that our reconstruction method accurately predicts the $\chi^{2}$ figure of merit. For reference, we find that the largest shift in the recovered mean and variance among the 6 cosmologically interesting parameters in Planck was 0.2% and 6% respectively.

To put the power of this new tool into perspective, we note that our reconstructed posterior was created using only 1200 training samples from the original MCMC chain. Sampling from the reconstructed distribution took approximately thirty minutes on a Macbook Pro laptop with a 6-core i7 processor. The speed offered by our method can be used for computationally demanding tasks such as marginalization and computing Bayesian evidences.

5 Conclusions

We present a method for reconstructing probability distributions given a small, carefully selected sample of points in parameter space at which the posterior is evaluated. This method can be utilized to reconstruct the posterior (or likelihood) function of an experiment given publicly available MCMC data. We demonstrate its utility by showcasing how this technique can be used to jointly sample distributions with partially overlapping parameter spaces. We also showed that we can reconstruct the Planck posterior distribution of the TT, TE, EE+lowE analysis presented in Planck Collaboration et al. (2018) using only 1200 samples from the original chain. We sample our reconstructed distribution and obtain an MCMC chain with forty times as many samples as the original chain in only thirty minutes.

The utility of being able to reconstruct posteriors from publicly available chains to enable joint analyses was the primary motivator for this work. However, perhaps the most interesting aspect of the tool we have developed is its potential for significantly reducing the computational effort necessary for generating MCMC samples in the first place. For instance, while our reconstruction of the Planck posterior was based on an existing Planck chain, one could have easily adopted a fiducial model in conjunction with a Fisher matrix analysis to determine a few thousand points at which to evaluate the likelihood. The reconstructed posterior can be used to run an MCMC, which can be tested for accuracy through sampling of the Latin Hypercube, and comparing the predicted and actual posterior at the sampled points. In this way, one can determine the regions of parameter space that require further evaluations in order to achieve accurate interpolation. Through iteration, we envision using this method for constructing the full posterior distribution with many fewer likelihood evaluations than brute-force MCMC sampling. This type of approach could be used in conjunction with other proposals such as likelihood-free inference and expectation optimization (Alsing et al., 2019; Seljak & Yu, 2019). We postpone an investigation of the feasibility of such an approach to future work.

Acknowledgements

TM thanks V. Miranda, Niall MacCrann, and Anže Slosar for helpful discussion. ER was supported by the DOE grant DE-SC0015975, grant FG-2016-6443, and the Cottrell Scholar program of the Research Corporation for Science Advancement. ER would like to thank Ann Lee for useful conversations that helped spark interest in this work.

Software: Python, Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), SciPy (Jones et al., 2001), emcee (Foreman-Mackey et al., 2013), george (Ambikasaran et al., 2015), ChainConsumer (Hinton, 2016).