Transformations in Semi-Parametric Bayesian Synthetic Likelihood

by   Jacob W. Priddle, et al.

Bayesian synthetic likelihood (BSL) is a popular method for performing approximate Bayesian inference when the likelihood function is intractable. In synthetic likelihood methods, the likelihood function is approximated parametrically via model simulations, and then standard likelihood-based techniques are used to perform inference. The Gaussian synthetic likelihood estimator has become ubiquitous in BSL literature, primarily for its simplicity and ease of implementation. However, it is often too restrictive and may lead to poor posterior approximations. Recently, a more flexible semi-parametric Bayesian synthetic likelihood (semiBSL) estimator has been introduced, which is significantly more robust to irregularly distributed summary statistics. In this work, we propose a number of extensions to semiBSL. First, we consider even more flexible estimators of the marginal distributions using transformation kernel density estimation. Second, we propose whitening semiBSL (wsemiBSL) – a method to significantly improve the computational efficiency of semiBSL. wsemiBSL uses an approximate whitening transformation to decorrelate summary statistics at each algorithm iteration. The methods developed herein significantly improve the versatility and efficiency of BSL algorithms.



There are no comments yet.


page 15

page 16

page 18


Efficient Bayesian synthetic likelihood with whitening transformations

Likelihood-free methods are an established approach for performing appro...

Adaptive MCMC for synthetic likelihoods and correlated synthetic likelihoods

Approximate Bayesian computation (ABC) and synthetic likelihood (SL) are...

BSL: An R Package for Efficient Parameter Estimation for Simulation-Based Models via Bayesian Synthetic Likelihood

Bayesian synthetic likelihood (BSL) is a popular method for estimating t...

Robust Bayesian Synthetic Likelihood via a Semi-Parametric Approach

Bayesian synthetic likelihood (BSL) is now a well established method for...

Finding our Way in the Dark: Approximate MCMC for Approximate Bayesian Methods

With larger data at their disposal, scientists are emboldened to tackle ...

Towards constraining warm dark matter with stellar streams through neural simulation-based inference

A statistical analysis of the observed perturbations in the density of s...

Generalized Bayesian Likelihood-Free Inference Using Scoring Rules Estimators

We propose a framework for Bayesian Likelihood-Free Inference (LFI) base...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Simulator models are a type of stochastic model that is often used to approximate a real-life process. Unfortunately, the likelihood function for simulator models is generally computationally intractable, and so obtaining Bayesian inferences is challenging. Approximate Bayesian computation (ABC) (Sisson et al., 2018a) and Bayesian synthetic likelihood (BSL) (Price et al., 2018; Wood, 2010) are two methods for approximate Bayesian inference in this setting. Both methods eschew evaluation of the likelihood by repeatedly generating pseudo-observations from the simulator, given an input parameter value. ABC and BSL methods have been applied in many different fields; recently, in biology, to model the spread of the Banana Bunchy Top Virus (Varghese et al., 2020); in epidemiology, to model the transmission of HIV (McKinley et al., 2018) and tuberculosis (Lintusaari et al., 2019), and, in ecology, to model the dispersal of little owls (Hauenstein et al., 2019)

. ABC is a more mature and established technique than BSL, and so it is more prevalent in applied fields. However, ABC can suffer from the curse of dimensionality with respect to the dimension of the summary statistic, requires a large number of model simulations, and the results can be highly dependent on a set of tuning parameters. BSL methods can be used to overcome many of these limitations.

Synthetic likelihood methods approximate the likelihood function with a tractable distribution; in contrast, ABC methods are effectively non-parametric (Blum and François, 2010). The original synthetic likelihood method of Wood (2010)

approximates the summary statistic likelihood with a Gaussian distribution and then uses a Markov chain Monte Carlo (MCMC) sampler for maximum likelihood estimation. Later,

Price et al. (2018)

consider the Gaussian synthetic likelihood in the Bayesian setting, and refer to their method as Bayesian synthetic likelihood. In practice, the Gaussian assumption of the summary statistic vector may be too restrictive, leading to a poor estimate of the likelihood, and then a poor estimate of the posterior. Herein, we refer to the Gaussian BSL method as standard BSL, denoted sBSL.

A few authors have considered more flexible density estimators to improve the robustness of sBSL to irregular summary statistic distributions (e.g. Papamakarios et al., 2018; An et al., 2020; Fasiolo et al., 2018). In particular, the semi-parametric Bayesian synthetic likelihood (semiBSL) method of An et al. (2020), estimates the intractable summary statistic likelihood semi-parametrically – non-parametrically estimating the marginal distributions using kernel density estimation (KDE), and parametrically estimating the dependence structure using a Gaussian copula. An et al. (2020) show empirically that semiBSL performs favourably to sBSL when summary statistics are distributed irregularly. semiBSL maintains many of the attractive properties of sBSL, including its scalability to a high dimensional summary statistic and ease of tuning.

Despite the appeal of semiBSL, the number of model simulations required to accurately estimate the correlation matrix scales poorly with the dimension of the summary statistic. The equivalent problem for sBSL, the scaling of the estimation of covariance matrix with the number of model simulations, has been explored by An et al. (2019), Ong et al. (2018a), Ong et al. (2018b), Everitt (2017), Frazier et al. (2019) and Priddle et al. (2020). However, there are currently no methods designed specifically for the semi-parametric estimator, which, in practice, may preclude its application to problems where model simulation is computationally expensive. The first contribution of this article adapts and extends the methodology presented in Priddle et al. (2020), which combines a whitening transformation with shrinkage covariance matrix estimation, to the semiBSL context.

SemiBSL provides additional robustness over sBSL when the summary statistic marginals deviate from normality. However, as we demonstrate in subsequent sections, for some distributions the KDE will fail. For instance, when a marginal summary statistic distribution has extremely heavy tails, the KDE will allocate essentially no density to the center of the distribution, and all weight to the tails (see Figure 2). In addition, it is well-known that the global bandwidth KDE rarely provides adequate smoothing over all features of the underlying distribution (Wand et al., 1991; Yang et al., 2003). Our second contribution addresses this problem with a procedure that draws upon and extends the vast body of literature on density estimation. Specifically, we consider transformation kernel density estimation (TKDE, Wand et al., 1991) to estimate the marginal distributions of the summary statistic. The idea is to transform the distribution so that the standard global bandwidth KDE is accurate, and then transform back to the original domain to estimate the density. We adapt the hyperbolic power transformation of Tsai et al. (2017), and propose a procedure to effectively apply TKDEs in a semiBSL algorithm.

The remainder of this article is structured as follows. In sections 2 and 3, we provide an overview of sBSL and semiBSL, respectively. In section 4, we present our method to significantly improve the computational efficiency of semiBSL. In section 5, we propose a new estimator of the marginal summary statistic distributions for semiBSL using TKDE. We assess the accuracy of the TKDEs on a number of test densities with known distribution. In section 6, we apply our new methods to four different examples. Last, we conclude in section 7.

2 Bsl

Synthetic likelihood algorithms are applicable in settings where the likelihood function is intractable but simulation from the model is straightforward, where (with ) is the set of observed data and is an unknown parameter. Here, our target is the posterior distribution , where is the prior distribution on the parameter. In synthetic likelihood, among other likelihood-free algorithms, such as approximate Bayesian computation (ABC) (see Sisson et al., 2018b), it is standard practice to degrade the data to a vector of informative summary statistics to help mitigate problems associated with dimensionality. Specifically, let be the summary statistic function that maps an -dimensional dataset to a -dimensional summary statistic. For , the implied target conditional on the summary statistic, often referred to as the partial posterior, is then ; depending (to a large extent) on the informativeness of the summary statistic, . However, since is intractable, it is generally the case that is also intractable, which leads us to consider sampling based methods that do not require evaluation of to obtain approximate inferences from the partial posterior.

In essence, synthetic likelihood methods assume a parametric form of the likelihood, which acts as a surrogate for the true likelihood and may be used directly in an MCMC (Markov chain Monte Carlo) sampler. In sBSL (see Price et al., 2018), the summary statistic likelihood is approximated with a Gaussian distribution, . The synthetic likelihood parameters and are typically unknown, but a series of independent and identically distributed simulations from the model with corresponding summary statistics can be used to construct the Monte Carlo estimates:


These may be used to yield the Gaussian synthetic likelihood estimator, , and the corresponding sBSL posterior approximation:

There are two main appeals of BSL: (1) that it can handle a relatively high dimensional summary statistic, and (2) that it can be more computationally efficient than competing likelihood-free Bayesian methods (Price et al., 2018; Frazier et al., 2019). These are both direct benefits of specifying a parametric form of the summary statistic likelihood. However, as demonstrated by An et al. (2020)

, in cases where the marginal summary statistic distributions deviate greatly from Gaussian, with, for example, heavy skewness, heavy tails or multiple modes, sBSL methods begin to break down. Often the posterior distribution will fail to adequately approximate the true partial posterior. In particularly challenging cases, the variance of the log synthetic likelihood estimator may be so large that the MCMC chain will become stuck within only a few iterations, and no discernible posterior distribution may be recovered (see Figure


3 semiBSL

In this section, we provide an overview of the semiBSL method of An et al. (2020). semiBSL provides additional robustness for a non-Gaussian distributed summary statistic. In semiBSL, the semi-parametric likelihood estimator is constructed as follows. Denote

the random variable corresponding to the

summary statistic. Given the set of

model simulations, the true PDF (probability density function)

is approximated using the kernel density estimate:


where and is the bandwidth. The kernel function may be any symmetric PDF; in semiBSL, the Gaussian kernel is used due to its simplicity and unbounded support. The above kernel density estimator uses a global (constant) bandwidth, selected according to the rule of Silverman (1986). It is straightforward to obtain the corresponding estimate of the CDF (cumulative density function) from the above equation.

Following estimation of the marginal summary statistic distributions, the dependence between the summaries is modelled via the Gaussian copula. Essentially, copula modelling allows the dependence structure and the marginal distributions to be estimated independently, allowing the user to consider alternative and more flexible marginal density estimators than the Gaussian distribution, as the case is in sBSL. For an introduction to copula models, we refer the reader to Trivedi et al. (2007). The Gaussian copula density,

is parameterised by the correlation matrix

and the vector of standard Gaussian quantiles

, where

is the inverse CDF of the standard normal distribution and

for . Replacing with its kernel density estimate evaluated at the observed summary , and with the estimated correlation matrix , we obtain the semiBSL posterior:

In the above equation, where for and is estimated using a collection of simulated summary statistics . In practice, An et al. (2020) advocate to estimate with the Gaussian rank correlation (GRC) (see Boudt et al., 2012), which provides additional robustness to the potential lack of fit of the KDEs.

We highlight two main limitations of semiBSL. First, the number of model simulations required to accurately estimate scales poorly with . This may be problematic for applications where model simulation is computationally expensive, especially if a relatively low dimensional and informative summary statistic is unavailable. Furthermore, the KDE is unreliable for distributions with extremely heavy tails, which may induce unduly high variance in the estimator and cause semiBSL to fail. In subsequent sections, we propose methods to overcome each of these limitations.

4 Whitening semiBSL

We now propose a method to improve the computational efficiency of semiBSL. Namely, we extend the whitening BSL (wBSL) methodology proposed by Priddle et al. (2020) to the semiBSL context. The motivation behind wBSL is articulated in Theorem 1 of Priddle et al. (2020). The main consequence of the theorem is that for a Gaussian log synthetic likelihood estimator with diagonal covariance structure, must scale linearly with to control the variance of the estimator. On the other hand, to control the variance of the traditional Gaussian log synthetic likelihood estimator (that estimates the full covariance structure), must scale quadratically with . This result suggests that there are significant computational benefits possible in BSL algorithms if the summary statistics are uncorrelated.

Despite such a compelling result, it is a challenging problem to find a summary statistic vector that is both independent across its dimensions and retains a large proportion of the information content intrinsic to the observed data. The main idea of wBSL is that an approximate whitening or decorrelation transformation may be applied to the summary statistic at each algorithm iteration. In doing so, the covariance shrinkage estimator of Warton (2008) may be applied with a high penalty, producing an accurate, low variance estimate of the likelihood function for a relatively small number of model simulations. If the full penalty is applied, this coincides with the Gaussian synthetic likelihood estimate with a diagonal covariance structure, and thus the desired computational gains may be achieved. In several empirical examples, Priddle et al. (2020) demonstrate that wBSL is able to produce an accurate partial posterior approximation, with an order of magnitude less model simulations than sBSL. Given the semi-parametric synthetic likelihood estimator uses the Gaussian copula, it is likely that it will inherit similar computational gains to the classical Gaussian estimator, particularly in cases where the marginal distributions are close to Gaussian. However, the extension of these concepts to semiBSL is not yet clear; here we provide an outline of our methodology, which we refer to as wsemiBSL.

Consider the Gaussian approximation of the summary statistic likelihood:

where the dependence of and on has been suppressed for notational convenience. It is straightforward to show that:

where and .

The main disparity between wBSL and wsemiBSL, is that in wsemiBSL the whitening transformation is applied to the standard Gaussian quantiles, and not directly to the summary statistics. We find that in the context of semiBSL, the latter approach does not produce as accurate posterior approximations (results not shown). Specifically, we propose to apply the whitening transformation to convert the random vector of arbitrary distribution with covaraince matrix into the transformed vector

for some whitening matrix , such that the covariance

is the identity matrix. Like in wBSL, we estimate the whitening matrix off-line using

independent model simulations such that given some carefully chosen parameter value with reasonable posterior support. Picchini et al. (2020) detail how Bayesian optimization may be used to rapidly generate a that has reasonable support under the posterior. This method, or the techniques described in Priddle et al. (2020), may be employed to find a suitable for our methods. is set high (much higher than ) to ensure an accurate estimate of is obtained. Of course, for the transformation to be exact, must evolve as a function of ; however, like in wBSL, we hold constant to preserve the target partial posterior obtained using semiBSL (when no penalty is applied), and so generally . Given the inverse transformation and Jacobian term , the summary statistic likelihood under the transformed variable is

where is the covariance matrix of the transformed quantiles . Of course, in semiBSL, we replace each marginal with the kernel density estimate and with a sample estimate. That is,

where and for . is estimated using simulated quantiles which constitute the rows of the matrix such that and for and . Given the whitening transformation approximately decorrelates the summary statistic quantiles, the Warton (2008) covariance shrinkage estimator

may be applied accurately with a high degree of shrinkage, where , is an estimate of the correlation matrix and is the shrinkage parameter. Effectively, is a constant that is multiplied by the off-diagonal elements of the sample covariance. Thus, shrinks the pairwise covariance elements to , assuming independent summary statistic quantiles. The heavier the shrinkage, the lower the value of required to precisely estimate the likelihood.

The choice of whitening matrix was considered carefully in Priddle et al. (2020). Any that satisfies will whiten the data at ; however, as the current parameter value deviates further from , the transformation will become less accurate. The most suitable for BSL is the one that most effectively decorrelates summary statistics generated by parameter values that reside in regions of the parameter space with non-negligible posterior density. Priddle et al. (2020) consider the five optimal whitening matrices of Kessy et al. (2018). Priddle et al. (2020)

find that in the context of BSL, principal components analysis (PCA) whitening produces the most accurate partial posterior approximations upon the application of heavy shrinkage. Thus, in wsemiBSL we also use the PCA whitening matrix,

where and

are the eigenvalue and eigenvector matrices of the covariance matrix

such that .

5 Transformation KDE in semiBSL

Our second contribution significantly improves the robustness of the semi-parametric estimator proposed in An et al. (2020) in the context of BSL. As demonstrated in Figure 2, if a given marginal summary statistic distribution has extremely heavy tails, as is common in financial applications for example (see Section 6.4), the standard KDE does not accurately approximate the true marginal distribution for reasonable sample sizes (number of model simulations in our context). We propose a new semi-parametric estimator that uses transformation kernel density estimation (see Wand et al., 1991) to model each marginal summary statistic. Like in the classic semiBSL estimator, we model the dependence between the summary statistic dimensions using the Gaussian copula. By doing so, the whitening method proposed in the previous section may be applied in conjunction with the new estimator, to achieve computational gains on top of the improved robustness. In this section, we provide details of our TKDE method for semiBSL.

Transformation kernel density estimation was introduced by Wand et al. (1991); although, the general ideas have been applied in many different contexts (see, for example, Kingma et al., 2016; Parno and Marzouk, 2018). In brief, the idea is to transform a sample of data so that the standard global bandwidth kernel density estimator (as in (3)) is more accurate, and then transform back to the original domain to obtain the estimate of the desired density.

Recall we are interested in estimating the marginal distributions of the summary statistic vector. That is, for the marginal , we wish to provide an estimate of the true density with support given access to our sample . Hereafter we suppress the notation for simplicity, and emphasise that we are considering a univariate distribution. Denote a family of bijective and differentiable transformations indexed by the parameter that map to the real line. The PDF of the transformed random variable is given by:

The value of is chosen so that is approximately Gaussian. Given this, KDE should provide an accurate approximation of the PDF on the transformed domain according to

An estimate of the density on the original domain is then obtained via the inverse transformation:

The above estimator can be thought of as using a location adaptive bandwidth on the original domain. This allows more appropriate smoothing over all features of the density, and often leads to a more accurate density approximation. Variable bandwidth methods, such as those proposed in Loftsgaarden et al. (1965) and Breiman et al. (1977) explicitly model the bandwidth as a function of the data. We find (results not shown) for the test distributions considered in this paper, that TKDE performs better.

A non-trivial aspect of applying TKDE to semiBSL is choosing an appropriate family of transformations, and then finding a method of efficiently estimating . The most suitable family of transformations is highly dependent on the shape of the data. Wand et al. (1991) focus on right-skewed data and use the shifted power transformation; Yang et al. (2003) use sequential transformations from the Johnson family to estimate the density of a wide range of distributions and Buch-Larsen et al. (2005) use the Champernowne transformation for heavy-tailed data. Our method can be extended to use any of these transformations (among others), but, due to its flexibility, we focus on the hyperbolic power transformation (HPT) introduced by Tsai et al. (2017), which has not previously been used in the TKDE context. The HPT is given by:

where is median centered, , and . are the power parameters; are the scale parameters, and is the normalising constant. By splitting the data either side of the median, the transformation is able to handle bimodal distributions, provided the modes are not well separated. As demonstrated by Tsai et al. (2017), the HPT outperforms other relevant normality transformations for a wide range of distributions.

There are many different optimality criteria possible to determine . Wand et al. (1991) and Yang et al. (2003) use asymptotic results based on minimising the mean integrated square error. Here we follow the approach used in Tsai et al. (2017) and use maximum likelihood estimation. That is, given we wish to transform the summary statistics such that the global bandwidth KDE will perform well, we target the standard normal distribution in our transformation. It can be shown that the objective function is given by:

where is the PDF of the standard normal distribution, and the Jacobian term is:

In practice, there are often several solutions to the score equation, however, only one is the global maximum. Tsai et al. (2017) employ the simplex method (see Nelder and Mead, 1965) to approximate the MLEs by iteratively optimising each split of the data separately and then perturbing the estimate of the slope parameter according to its MLE:

In our implementation, we take a similar approach to Tsai et al. (2017) in splitting the data and estimating each of the pairs and separately. However, we update the value of using the MLE (using the relevant split of the data) for each evaluation of the likelihood, due to its dependence on the other parameters. We find that this approach works well without having to iteratively maximise the parameters and perturb . This is crucial in the context of semiBSL as each iteration of MCMC will involve an estimate of the synthetic likelihood at the proposed parameter value. Of course, our method only serves as an approximation of the true maximum, but as we shall demonstrate, this is sufficient to significantly improve the accuracy of the density estimate over standard KDE. Each marginal summary statistic distribution may be estimated in parallel, meaning the overall additional computational time is small. We use the quantile approach outlined in Tsai et al. (2017) to initialise for each optimisation problem. Alternatively, optimal parameters found at previous iterations may be used to inform initial parameter values at subsequent iterations.

Despite the appeal of the HPT, we find that for very heavy tailed data, the transformation is not numerically stable. It is also a non-trivial task to reparameterise the transformation such that it is numerically stable. Therefore, we propose an extension of the HPT that uses a series of transformations to first reduce the heaviness of tails, allowing the HPT to subsequently be applied more effectively. The

transformations do not require any estimation of parameters, and so they add negligible computation time. For positively skewed data with heavy kurtosis, we use the transformation:

where if , otherwise . Analogously, we use the following transformation for negatively skewed data with heavy kurtosis:

where if , otherwise . Lastly, for symmetric data with heavy kurtosis we use

When one of the above three transformations is applied concurrently with the HPT, we refer to each method as semiBSL TKDE1, semiBSL TKDE2, or semiBSL TKDE3, respectively. SemiBSL TKDE0 refers to semiBSL TKDE without an initial transformation, and semiBSL TKDE is the general method of using transformation kernel density estimation for semiBSL. We emphasise that the main purpose of the transformations is to transform the data such that the HPT can be accurately computed, not to transform the data to Gaussian. It is the HPTs job to Gaussianise the transformed data. Figure 1 shows the estimated density after each step of the TKDE procedure for several test densities with known PDF. Each test density is close to standard Gaussian after applying the HPT (row 3, Figure 1), and the final density estimate is close to the true density (row 4, Figure 1). For our semiBSL TKDE method, we recommend the user perform a number of model simulations at , visualise the marginal summary statistic distributions and then decide whether or not (and which) transformation is necessary.

To illustrate the efficacy of the proposed density estimation procedure, we perform a simulation study using a wide a range of distributions with known density. This follows directly from the work in An et al. (2020). Specifically, we assume the observed data is drawn from the standard Gaussian distribution , and the summary statistic is given by (this is the sinh-archsinh transformation of Jones and Pewsey, 2009). and control the skewness and kurtosis respectively. Here we choose the values of and to reflect the shapes of densities that arise in practice, for example, in the models of Section 6. We also consider an observed dataset drawn directly from a bimodal Gaussian distribution, such that and take . For each test density, we estimate the PDF using KDE and TKDE for , and . For TKDE, we show the results using the most appropriate

transformation (or lack thereof), see Figure 2. Furthermore, we estimate the total variation distance between the true and estimated PDFs using numerical integration over a grid of parameter values based on 1000 independent replicates of the above procedure. We report the sample mean and standard deviation of the 1000 total variation distances (see Table 1). The total variation between two PDFs

and is given by .

Figure 1: Intermediate densities of TKDE procedure for various test densities. Each row corresponds to a step in the density estimation: row 1 is a histogram of the original data; row 2 is a KDE after the log transformation; row 3 is a KDE after the HPT (with the standard normal distribution in black) and row 4 is the final density estimate on the original domain (with the true PDF shown in black). Columns correspond to each test density: skewness and kurtosis (; left), kurtosis only (), skewness only (), bimodal (), heavy skewness () and skewness with heavy kurtosis (; right).
Figure 2: Comparison of density estimators (KDE and TKDE) with the true density. Rows correspond to (top row), and (bottom row) model simulations. Columns correspond to the same test densities listed in Figure 1.
Skewness and Kurtosis 0.201 (0.027) 0.101 (0.033) 0.138 (0.014) 0.053 (0.012) 0.116 (0.010) 0.041 (0.008)
Kurtosis                   0.162 (0.022) 0.095 (0.030) 0.099 (0.011) 0.050 (0.011) 0.079 (0.008) 0.039 (0.008)
Skewness                   0.136 (0.023) 0.072 (0.025) 0.094 (0.011) 0.038 (0.011) 0.080 (0.008) 0.030 (0.007)
Bimodal                   0.253 (0.028) 0.175 (0.032) 0.189 (0.010) 0.121 (0.015) 0.159 (0.007) 0.100 (0.011)
Heavy Kurtosis                   0.166 (0.013) 0.058 (0.033) 0.163 (0.008) 0.026 (0.011) 0.165 (0.006) 0.019 (0.007)
Skewness and Heavy Kurtosis 0.044 (0.011) 0.014 (0.009) 0.023 (0.005) 0.007 (0.004) 0.018 (0.004) 0.006 (0.003)
Table 1: Total variation distances between density estimators (KDE and TKDE) and the true density. The mean is reported for each of the four test densities using , and model simulations. The corresponding standard deviations are given in parentheses.

Figure 2 demonstrates that the proposed TKDE scheme is able get much closer to the true PDF than standard KDE, even with a small number of model simulations. The TKDE nicely captures the peaks of each distribution, and provides adequate smoothing over the tails. For , the KDE appears completely flat due to the extremely heavy tails, whereas the TKDE is very accurate. TKDE also outperforms KDE for the bimodal test density, with a noticeably better performance for . Interestingly, in some cases, we find that the transformation is detrimental to the TKDE, and so the user must carefully decide whether or not the transformation is needed. The simulation results in Table 1 support the above findings, with all TKDEs having a lower total variation distance than the corresponding KDE. The benefits of TKDE are most apparent for heavy tailed distributions.

6 Results

In this section, we apply our methods to four examples. The examples, and what they are designed to demonstrate are listed below:

  1. MA example: demonstrates the potential efficiency gains of wsemiBSL; the robustness of semiBSL TKDE, and the simultaneous use of whitening and TKDE in semiBSL (wsemiBSL TKDE hereafter) for improved effiency and robustness.

  2. Fowler’s Toads example: demonstrates the potential efficiency gains of wsemiBSL.

  3. M/G/1 example: demonstrates the improved robustness of semiBSL TKDE.

  4. -stable stochastic volatility model: demonstrates the improved robustness of semiBSL TKDE.

The likelihood for the MA example is known, allowing us to compare the result of our methods to the output of a Metropolis-Hastings sampler that uses the true likelihood. Each of the remaining three models have an intractable likelihood and are representative of a real-life modelling scenario.

In all cases, we use the Metropolis-Hastings algorithm with a Gaussian random walk. The random walk covariance matrix is set to be roughly the (approximate) posterior covariance obtained from pilot runs. Unless stated otherwise, the value of is tuned such that the standard deviation of the log synthetic likelihood evaluated at is in the range , as Price et al. (2018) find that this maximises the computational efficiency of sBSL. We compare posterior approximations using the total variation distance, as described in Section 5. For wsemiBSL, we use to accurately estimate . Each sampler is run for MCMC iterations.

6.1 Ma(2)

The observation in a moving average process of order 2, denoted MA, evolves according to:

subject to the constraints , and . It is straightforward to show that the likelihood is Gaussian with zero mean vector and pentadiagonal covariance matrix, with entries given by: , and , where . The MA model is commonly used as a toy example to demonstrate likelihood-free methods (see, Chiachio et al., 2014; Marin et al., 2012; Frazier et al., 2019).

We simulate 50 observations from the MA process at and set this to be our observed data, such that . We assume that is known, and equal to 1. For semiBSL, we are interested in cases where the marginal summary statistic distributions deviate from Gaussian. As in Section 5, we use the sinh-archsinh transformation of Jones and Pewsey (2009) to transform and generate a summary statistic with non-Gaussian marginals; thus, , where is the sinh-archsinh transformation applied elementwise. We use a uniform prior over the parameter space.

We first test our methods with a summary statistic generated with 4 different and combinations. We consider , which corresponds to no transformation; , which creates negative skewness; , which creates positive kurtosis and , which creates negative kurtosis and positive skewness. For each of these summary statistics, we consider the following methods: semiBSL (equivalent to wsemiBSL with ); wsemiBSL with ; semiBSL TKDE () and wsemiBSL TKDE with . We compare the results to the ‘true’ posterior, which is obtained using an MCMC sampler with the true likelihood.

Posterior approximations are shown in Figure 3. Comparing columns 1 with 3 (no shrinkage, ), and columns 2 with 4 (complete shrinkage, ), it can be seen that the posterior approximations obtained with TKDE are generally more accurate in terms of the total variation distance to the ‘true’ posterior. The only case where the posterior approximation obtained using TKDE is less accurate than the corresponding estimate that uses KDE (albeit slightly, with tv distances of 0.2 and 0.16, respectively), is when and the marginal summary statistics have negative kurtosis (; row 3 of Figure 3).

The bivariate posterior approximations obtained using wsemiBSL with complete shrinkage are good approximations of the ‘true’ posterior in all cases (Figure 3). Interestingly, we find that there is a quite a strong dependence between the regularity of the marginal summary statistic distributions and the capacity of wsemiBSL to significantly reduce the number of model simulations. For no summary statistic transformation, the skewness transformation and the skewness and kurtosis transformation, wsemiBSL is extremely effective – allowing us to reduce by about an order of magnitude. However, for the kurtosis transformation, we are only able to reduce by a factor of three, while accurately estimating the log synthetic likelihood. In addition, we find that wsemiBSL is generally not as effective in reducing the number model simulations when TKDE is used compared to when KDE is used, to estimate the marginal summary statistic distributions. This is the case for all summary statistics except for the kurtosis transformation, where could be reduced further (from to ) when TKDE is used compared to standard KDE.

We consider two additional summary statistics with extremely heavy kurtosis. Specifically, we set and , which creates negative kurtosis, and also and , which creates heavy negative kurtosis and positive skewness. This presents an extremely challenging example for standard semiBSL. We find that is required for semiBSL TKDE for both datasets and is required for semiBSL when and . However, we are unable to find an that can accurately estimate the synthetic likelihood when and , since the KDE completely fails even for a huge number of model simulations (due to the heaviness of the tails). We also consider for semiBSL, representing the same number of model simulations used for semiBSL TKDE. For these examples, wsemiBSL is ineffective at reducing the required number of model simulations as the marginal summary statistics deviate too far from Gaussian and the pairwise correlation is low.

Bivariate posterior approximations are shown in Figure 4. For , it can be seen that standard semiBSL completely fails, while semiBSL TKDE produces an accurate posterior approximation, for both summary statistics. From Figure 5, it can be seen that the acceptance rates are much higher for semiBSL TKDE than standard semiBSL. For and for standard semiBSL, the variance of the log synthetic likelihood is so high that no samples are accepted. When model simulations are used for semiBSL, the parameter space appears to be explored well (Figure 5), but the posterior approximation is far less accurate than the semiBSL TKDE method ( compared to ), which only used model simulations.

Figure 3: Bivariate posterior approximations and true marginal summary statistic distributions for the MA(2) example – plot 1. Columns denote (left to right) the true marginal summary statistic distribution, semiBSL KDE, wsemiBSL with , semiBSL TKDE and wsemiBSL TKDE with . Each row uses the same marginal summary statistics. Black contours correspond to be the output of an M-H sampler using the known likelihood, green contours are for and blue contours are for . The parameter used to generate the observed data is shown in red. tv denotes the total variation distance between the approximate and true bivariate posterior distributions.
Figure 4: Bivariate posterior approximations and true marginal summary statistic distributions for the MA(2) example – plot 2. Similar to Figure 3, the columns (left to right) denote the true marginal summary statistic distributions, semiBSL with , semiBSL TKDE with and semiBSL with .
Figure 5: Trace plots corresponding to the results in Figure 4 (in the respective order). is green and is red.

6.2 Fowler’s Toads

The next example we consider is the individual-based movement model of Fowler’s Toads (Anaxyrus fowleri) developed by Marchand et al. (2017). The model has since been considered as a test example in likelihood-free literature by several authors (see An et al., 2020; Frazier and Drovandi, 2020; Priddle et al., 2020). Marchand et al. (2017) consider three models, each assuming that toads take refuge during the day and forage throughout the night. The models differ in their returning behaviour; here we expressly consider the random return model. We provide only a brief overview of the model herein, and refer the reader to Marchand et al. (2017) for more details.

To simulate from the model, we draw an overnight displacement from the Levy alpha-stable distribution , where and . At the end of the night, toads return to their previous refuge site with probability , or take refuge at their current overnight displacement. In the event of a return on day , the refuge site is chosen randomly from the set of previous refuge sites, thereby giving higher weighting to sites that have been visited multiple times. Here is the refuge locations of toads over days, generated at .

The summary statistic is 48-dimensional, and is constructed as follows. For each toad, we split the observed data in two, corresponding to displacements less than or greater than 10m. We count the number of absolute displacements less than 10m. For the latter dataset, we find the distance moved distribution at time lags 1, 2, 4 and 8 days, and compute both the log of the differences in the quantiles and the median. For this example, the marginal summary statistic distributions are roughly Gaussian (see Appendix A, Figure 10), meaning sBSL or wBSL would likely perform well. However, semiBSL (and wsemiBSL) will provide additional robustness over their Gaussian counterparts with little additional computation and so we would generally advocate to use these methods even for such models. Of course, TKDE is not necessary for this example.

We find that model simulations is adequate for standard semiBSL. We compare the output of standard semiBSL to wsemiBSL with (), () and () – results are shown in Figure 6. For all cases, the wsemiBSL posterior approximation is close to the standard semiBSL posterior approximation. With complete shrinkage (), we are able to reduce the number of model simulations by an order of magnitude.

Figure 6: Contour plots of the bivariate posterior approximations for the toad model. Rows (top to bottom) correspond to the and combination, and columns correspond to each pair of parameters. is shown as a red dot. The total variation distance between each bivariate semiBSL posterior approximation and each bivariate wsemiBSL posterior approximation is shown in the bottom right of each panel.

6.3 M/g/1

The M/G/1 queueing model is a stochastic single-server queue model whereby ‘customers’ arrive according to a Poisson process and service times have a general distribution. Here we expressly consider the case where service times are , as this has been a popular choice in other likelihood-free literature (see e.g. An et al., 2020; Blum and François, 2010). The time between arrivals is distributed. We assume that only the inter-departure times are known, and take this to be the observed data . We observe 50 inter-departure times (corresponding to 51 customers) and set , generated at . The prior is on .

The marginal summary statistic distributions are right skewed with moderate kurtosis (see Appendix A, Figure 11). Thus, for our TKDE method, it would be reasonable to use semiBSL TKDE0 or semiBSL TKDE1. wsemiBSL does not provide additional benefit for this example since the summary statistics have very low correlation. We run semiBSL TKDE0, semiBSL TKDE1 and standard semiBSL. We compare the results of each sampler to the ‘true’ posterior, obtained using the MCMC scheme for the M/G/1 queue model of Shestopaloff and Neal (2014). We use to estimate the summary statistic likelihood for semiBSL.

Bivariate posterior approximations are shown in Figure 7. Both semiBSL TKDE methods produce more accurate posterior approximations than standard semiBSL. semiBSL TKDE estimates the marginal distribution more accurately than semiBSL; the and marginals are similar. The additional transformation (in semiBSL TKDE1 compared to semiBSL TKDE0) slightly improves the accuracy of the posterior approximation for this example, as evidenced by the total variation distance.

Figure 7: Contour plots of the bivariate posterior approximations for the M/G/1 example. Columns (left to right) correspond to each bivariate marginal , and , respectively. Rows correspond to the method used. ‘Exact’ denotes the posterior approximation obtained using the method of Shestopaloff and Neal (2014). is shown as a red dot.

6.4 -Stable Stochastic Volatility Model

Stochastic volatility models (SVMs) are commonly used in econometric applications, such as the modelling of financial returns (see Vankov et al., 2019). In SVMs, the observed data are assumed to follow a latent stochastic process in evenly spaced discrete time. The return process is given by:

where is the observed data at time , which directly depends on the log volatility and the shock ; is the modal instantaneous volatility; is the persistence parameter, and is the variance of . The typical model formulation uses a Gaussian shock parameter, (Kim et al., 1998); however, due to the heavy tailedness of asset returns, more recent studies have found the stable distribution to be more appropriate (Casarin, 2004). That is, we assume , where , , and control the tail heaviness (with a lower having heavier tails), the skewness, the scale and the location, respectively. Despite the additional flexibility inherited by this family of SVMs, the PDF of the stable distribution is unavailable in closed form for most parameter values. This motivates the development of likelihood-free algorithms such as ABC and BSL for heavy tailed distributions (see, for example, Ebert et al., 2019; Vankov et al., 2019). The extremely heavy tails of may cause sBSL and standard semiBSL to fail.

We test our methods on two datasets. We infer and assume the remaining parameters are known and fixed, such that: , , , and for each dataset. We set and for datasets 1 and 2, respectively and generate 50 observations from the stable SVM and set this to be in each case. We take . Given the marginal summary statistic distributions are symmetric and heavily skewed (see Figure 8), we use semiBSL TKDE3. We do not consider wsemiBSL for this model, since there is only a low degree of correlation between the pairwise statistics. The results are compared directly to standard semiBSL. We find is sufficient to control the variance for semiBSL TKDE for each dataset, and is required for semiBSL for the dataset. We are unable to find a large enough to accurately estimate the standard semi-parametric synthetic likelihood for the dataset. Similar to the MA example, we also consider for semiBSL for each dataset – the same value of we use for semiBSL TKDE.

Marginal posterior approximations are shown in Figure 8. The corresponding trace plots are shown in Figure 9. We observe similar results to the MA example. For , for standard semiBSL, the acceptance rate is low (extremely low for dataset 2), while we observe high acceptance rates and good mixing for semiBSL TKDE for both datasets. The posterior approximations for dataset 1 obtained using semiBSL are reasonable, but are poor for dataset 2. On the other hand, the posterior approximations for semiBSL TKDE for each dataset are smooth and have reasonable support for the true parameter value. When is increased to 20000 for standard semiBSL, the posterior approximation is smoother; however there is less support for the true parameter value than the posterior approximation generated using semiBSL TKDE.

Figure 8: Posterior approximations and marginal summary statistic distributions for the stable SVM. Top row corresponds to dataset 1, and bottom row corresponds to dataset 2. The columns (left to right) correspond to the marginal summary statistic distributions, and the parameters and , respectively.
Figure 9: Trace plots corresponding to Figure 8 results. Rows correspond to each dataset – top row when and bottom row when . Columns correspond to the method and number of model simulations combination.

7 Discussion

In this article, we proposed two extensions to semiBSL. First, we extended the wBSL method of Priddle et al. (2020) to the semiBSL context. We demonstrated in a number of empirical examples that our new method, wsemiBSL, is able to produce accurate posterior approximations with up to an order of magnitude less model simulations than standard semiBSL, even when the summary statistic deviates from normality. We also proposed a new method to estimate the marginal summary statistic distributions in semiBSL using TKDE. We show several examples where standard semiBSL will fail due to heavy kurtosis in the marginal summary statistic distributions, whereas our semiBSL TKDE method produces accurate posterior approximations in each case. Furthermore, we showed that wsemiBSL can be used in conjunction with TKDE for both improved computational efficiency and robustness to irregular summary statistic distributions.

There are a few limitations to the proposed methods. For wsemiBSL, we find that there is a rather strong dependence between the regularity of the marginal summary statistic distributions and the potential for large reductions in . That is, the efficiency gain appears to diminish as the marginal summary statistic distributions become increasingly non-Gaussian.

In future work, it may be of interest to consider ways to further increase the robustness of semiBSL to non-linear dependence structures. One way of overcoming such a problem may be via more advanced multivariate transformations such as normalising flows (see Rezende and Mohamed, 2015; Papamakarios et al., 2017). In addition, sBSL is known to be adversely affected in the setting of model misspecification, or more specifically, summary statistic incompatibility (see Frazier and Drovandi, 2020). Future work may investigate the equivalent problem in the context of semiBSL.


  • An et al. (2020) An, Z., Nott, D. J., and Drovandi, C. (2020). Robust Bayesian synthetic likelihood via a semi-parametric approach. Statistics and Computing, 30(3):543–557.
  • An et al. (2019) An, Z., South, L. F., Nott, D. J., and Drovandi, C. C. (2019). Accelerating Bayesian synthetic likelihood with the graphical lasso. Journal of Computational and Graphical Statistics, 28(2):471–475.
  • Blum and François (2010) Blum, M. G. and François, O. (2010).

    Non-linear regression models for approximate Bayesian computation.

    Statistics and Computing, 20(1):63–73.
  • Boudt et al. (2012) Boudt, K., Cornelissen, J., and Croux, C. (2012). The Gaussian rank correlation estimator: robustness properties. Statistics and Computing, 22(2):471–483.
  • Breiman et al. (1977) Breiman, L., Meisel, W., and Purcell, E. (1977). Variable kernel estimates of multivariate densities. Technometrics, 19(2):135–144.
  • Buch-Larsen et al. (2005) Buch-Larsen, T., Nielsen, J. P., Guillén, M., and Bolancé, C. (2005). Kernel density estimation for heavy-tailed distributions using the Champernowne transformation. Statistics, 39(6):503–516.
  • Casarin (2004) Casarin, R. (2004). Bayesian inference for generalised Markov switching stochastic volatility models. CEREMADE Journal Working Paper, (0414).
  • Chiachio et al. (2014) Chiachio, M., Beck, J. L., Chiachio, J., and Rus, G. (2014). Approximate Bayesian computation by subset simulation. SIAM Journal on Scientific Computing, 36(3):A1339–A1358.
  • Ebert et al. (2019) Ebert, A., Pudlo, P., Mengersen, K., and Wu, P. (2019). Combined parameter and state inference with automatically calibrated ABC. arXiv preprint arXiv:1910.14227.
  • Everitt (2017) Everitt, R. G. (2017). Bootstrapped synthetic likelihood. arXiv preprint arXiv:1711.05825.
  • Fasiolo et al. (2018) Fasiolo, M., Wood, S. N., Hartig, F., Bravington, M. V., et al. (2018). An extended empirical saddlepoint approximation for intractable likelihoods. Electronic Journal of Statistics, 12(1):1544–1578.
  • Frazier and Drovandi (2020) Frazier, D. T. and Drovandi, C. (2020). Robust approximate Bayesian inference with synthetic likelihood. arXiv preprint arXiv:1904.04551.
  • Frazier et al. (2019) Frazier, D. T., Nott, D. J., Drovandi, C., and Kohn, R. (2019). Bayesian inference using synthetic likelihood: asymptotics and adjustments. arXiv preprint arXiv:1902.04827.
  • Hauenstein et al. (2019) Hauenstein, S., Fattebert, J., Grüebler, M. U., Naef-Daenzer, B., Pe’er, G., and Hartig, F. (2019). Calibrating an individual-based movement model to predict functional connectivity for little owls. Ecological Applications, 29(4):e01873.
  • Jones and Pewsey (2009) Jones, M. C. and Pewsey, A. (2009). Sinh-arcsinh distributions. Biometrika, 96(4):761–780.
  • Kessy et al. (2018) Kessy, A., Lewin, A., and Strimmer, K. (2018). Optimal whitening and decorrelation. The American Statistician, 72(4):309–314.
  • Kim et al. (1998) Kim, S., Shephard, N., and Chib, S. (1998). Stochastic volatility: likelihood inference and comparison with ARCH models. The Review of Economic Studies, 65(3):361–393.
  • Kingma et al. (2016) Kingma, D. P., Salimans, T., Jozefowicz, R., Chen, X., Sutskever, I., and Welling, M. (2016). Improved variational inference with inverse autoregressive flow. In Advances in Neural Information Processing Systems, pages 4743–4751.
  • Lintusaari et al. (2019) Lintusaari, J., Blomstedt, P., Rose, B., Sivula, T., Gutmann, M. U., Kaski, S., and Corander, J. (2019). Resolving outbreak dynamics using approximate Bayesian computation for stochastic birth–death models. Wellcome Open Research, 4(14):14.
  • Loftsgaarden et al. (1965) Loftsgaarden, D. O., Quesenberry, C. P., et al. (1965). A nonparametric estimate of a multivariate density function. The Annals of Mathematical Statistics, 36(3):1049–1051.
  • Marchand et al. (2017) Marchand, P., Boenke, M., and Green, D. M. (2017). A stochastic movement model reproduces patterns of site fidelity and long-distance dispersal in a population of Fowler’s toads (Anaxyrus fowleri). Ecological Modelling, 360:63–69.
  • Marin et al. (2012) Marin, J.-M., Pudlo, P., Robert, C. P., and Ryder, R. J. (2012). Approximate Bayesian computational methods. Statistics and Computing, 22(6):1167–1180.
  • McKinley et al. (2018) McKinley, T. J., Vernon, I., Andrianakis, I., McCreesh, N., Oakley, J. E., Nsubuga, R. N., Goldstein, M., White, R. G., et al. (2018). Approximate Bayesian computation and simulation-based inference for complex stochastic epidemic models. Statistical Science, 33(1):4–18.
  • Nelder and Mead (1965) Nelder, J. A. and Mead, R. (1965). A simplex method for function minimization. The Computer Journal, 7(4):308–313.
  • Ong et al. (2018a) Ong, V. M., Nott, D. J., Tran, M.-N., Sisson, S. A., and Drovandi, C. C. (2018a). Likelihood-free inference in high dimensions with synthetic likelihood. Computational Statistics & Data Analysis, 128:271–291.
  • Ong et al. (2018b) Ong, V. M., Nott, D. J., Tran, M.-N., Sisson, S. A., and Drovandi, C. C. (2018b). Variational Bayes with synthetic likelihood. Statistics and Computing, 28(4):971–988.
  • Papamakarios et al. (2017) Papamakarios, G., Pavlakou, T., and Murray, I. (2017). Masked autoregressive flow for density estimation. In Advances in Neural Information Processing Systems, pages 2338–2347.
  • Papamakarios et al. (2018) Papamakarios, G., Sterratt, D. C., and Murray, I. (2018). Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. arXiv preprint arXiv:1805.07226.
  • Parno and Marzouk (2018) Parno, M. D. and Marzouk, Y. M. (2018). Transport map accelerated Markov chain Monte Carlo. SIAM/ASA Journal on Uncertainty Quantification, 6(2):645–682.
  • Picchini et al. (2020) Picchini, U., Simola, U., and Corander, J. (2020). Adaptive MCMC for synthetic likelihoods and correlated synthetic likelihoods. arXiv preprint arXiv:2004.04558.
  • Price et al. (2018) Price, L. F., Drovandi, C. C., Lee, A., and Nott, D. J. (2018). Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics, 27(1):1–11.
  • Priddle et al. (2020) Priddle, J. W., Sisson, S. A., Frazier, D. T., and Drovandi, C. (2020). Efficient Bayesian synthetic likelihood with whitening transformations. arXiv preprint arXiv:1909.04857.
  • Rezende and Mohamed (2015) Rezende, D. J. and Mohamed, S. (2015). Variational inference with normalizing flows. arXiv preprint arXiv:1505.05770.
  • Shestopaloff and Neal (2014) Shestopaloff, A. Y. and Neal, R. M. (2014). On Bayesian inference for the M/G/1 queue with efficient MCMC sampling. arXiv preprint arXiv:1401.5548.
  • Silverman (1986) Silverman, B. W. (1986). Density estimation for statistics and data analysis, volume 26. CRC press.
  • Sisson et al. (2018a) Sisson, S. A., Fan, Y., and Beaumont, M. A. (2018a). Handbook of Approximate Bayesian Computation. Chapman and Hall.
  • Sisson et al. (2018b) Sisson, S. A., Fan, Y., and Beaumont, M. A. (2018b). Overview of approximate Bayesian computation. In Sisson, S. A., Fan, Y., and Beaumont, M. A., editors, Handbook of Approximate Bayesian Computation, pages 3–54. Chapman and Hall/CRC Press.
  • Trivedi et al. (2007) Trivedi, P. K., Zimmer, D. M., et al. (2007). Copula modeling: an introduction for practitioners. Foundations and Trends® in Econometrics, 1(1):1–111.
  • Tsai et al. (2017) Tsai, A. C., Liou, M., Simak, M., and Cheng, P. E. (2017). On hyperbolic transformations to normality. Computational Statistics & Data Analysis, 115:250–266.
  • Vankov et al. (2019) Vankov, E. R., Guindani, M., Ensor, K. B., et al. (2019). Filtering and estimation for a class of stochastic volatility models with intractable likelihoods. Bayesian Analysis, 14(1):29–52.
  • Varghese et al. (2020) Varghese, A., Drovandi, C., Mira, A., and Mengersen, K. (2020). Estimating a novel stochastic model for within-field disease dynamics of banana bunchy top virus via approximate Bayesian computation. PLOS Computational Biology, 16(5):e1007878.
  • Wand et al. (1991) Wand, M. P., Marron, J. S., and Ruppert, D. (1991). Transformations in density estimation. Journal of the American Statistical Association, 86(414):343–353.
  • Warton (2008) Warton, D. I. (2008). Penalized normal likelihood and ridge regularization of correlation and covariance matrices. Journal of the American Statistical Association, 103(481):340–349.
  • Wood (2010) Wood, S. N. (2010). Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310):1102.
  • Yang et al. (2003) Yang, C., Duraiswami, R., Gumerov, N. A., and Davis, L. (2003). Improved fast Gauss transform and efficient kernel density estimation. IEEE.

Appendix A Summary Statistic Distributions

Figure 10: Marginal summary statistic distributions for the Fowler’s toads example.
Figure 11: Marginal summary statistic distributions for the M/G/1 example.