Adaptive Gaussian Copula ABC
Yanzhi Chen Michael U. Gutmann























apalike



























       






bayesnet
math
arrows.meta











[2][]
  #1  #2



‖‖

 


 





































 



 



theorem[2][Theorem]

  #1#2.
definition[2][Definition]

  #1#2.
proposition[2][Proposition]

  #1#2.
lemma[2][Lemma]

  #1#2.
exercise[2][Exercise]

  #1#2.
solution[2][Solution]

  #1
problem[2][Problem]

  #1#2.
question[2][Question]

  #1#2.
corollary[2][Corollary]

  #1#2.



thmTheorem
TheoremTheorem
LemmaLemma
lemLemma
CorollaryCorollary
PropositionProposition
RemarkRemark
ProofProof




[

Adaptive Gaussian Copula ABC

 Yanzhi Chen Michael U. Gutmann
The University of Edinburgh The University of Edinburgh  ]


   is a set of techniques for Bayesian inference when the likelihood is intractable but sampling from the model is possible. This work presents a simple yet effective ABC algorithm based on the combination of two classical ABC approaches — regression ABC and sequential ABC. The key idea is that rather than learning the posterior directly, we first target another auxiliary distribution that can be learned accurately by existing methods, through which we then subsequently learn the desired posterior with the help of a Gaussian copula. During this process, the complexity of the model changes adaptively according to the data at hand. Experiments on a synthetic dataset as well as three real-world inference tasks demonstrates that the proposed method is fast, accurate, and easy to use. 





§ INTRODUCTION

Many parametric statistical models are specified in terms of a parametrised data generating process. We can sample or simulate data from this kind of models but their likelihood function is typically too costly to evaluate. The models are called implicit <cit.> or simulator-based models <cit.> and are widely used in scientific domains including ecology <cit.>,  epidemiology <cit.>, psychology <cit.> and cosmology <cit.>. For example, the demographic evolution of two species can be simulated by a set of stochastic differential equations controlled by birth/predation rates but computation of the likelihood is intractable. Often, the true interest is rarely in simulating data, but in the inverse problem of identifying model parameters that could have plausibly produced the observed data. This usually not only includes point estimates of the model parameters but also a quantification of their uncertainty. 

The above task can, in principle, be solved by the well-known framework of Bayesian inference. However, due to the absence of likelihood functions for implicit models, exact inference becomes difficult or even impossible. The technique of approximate Bayesian computation (ABC) enables inference in such circumstances <cit.>. The most basic ABC algorithm works by repeatedly simulating data with different parameters and only accepting those parameters whose simulated data 'resemble' the observed ones. This generally requires hitting a small ϵ-ball for good accuracy, which is unfortunately impractical in high-dimensional spaces due to the well-known curse of dimensionality <cit.>. 

Several efforts have been made to alleviate the above accuracy-efficiency dilemma in ABC <cit.>. The first line of work, regression ABC <cit.>, accelerates inference by first allowing a much larger ϵ-ball, and then capturing the relationship between the simulated data and the parameters inside this ϵ-ball through a regression function, with which accepted samples are post-adjusted so as to compensate for the accuracy loss due to the larger ϵ used. Another line of work improves the efficiency of ABC by learning the posterior iteratively <cit.>. In these methods, preliminary models of the posterior are used to guide further simulations, leading to a significant drop of overall simulation costs. Both lines of methods are shown to be effective in a wide range of inference tasks. 

In this paper, we show how these two lines of works can refine each other, yielding a powerful yet simple algorithm. We bridge the two lines of work using Gaussian copula modelling <cit.>. The resulting proposed approach, which we call adaptive Gaussian copula ABC (AGC-ABC), has similar power as recent flexible machine-learning methods <cit.> but requires less simulated data. This naturally leads to performance gains for smaller computational budgets, as we confirm empirically.

The rest of the paper is organized as follows. Section 2 gives an overview of regression ABC and sequential ABC. Section 3 details the proposed AGC-ABC algorithm. Section 4 compares the proposed method to existing works and Section 5 contains the experimental results.  Section 6 concludes the paper. 


§ BACKGROUND

Given observed data , the task of approximate Bayesian inference for implicit models is to estimate the posterior density of the  parameters  without evaluating the likelihood directly. To achieve this goal, methods known as approximate Bayesian computation (ABC) typically reduce the observed data to summary statistics ^o = () and approximate the likelihood by ( - ^o < ϵ |  ) where ϵ is chosen sufficiently small. The approximate posterior is then

    ( |^o)    ∝π() ·( - ^o < ϵ |).

ABC algorithms typically only yield samples from the above approximate posterior. The most basic algorithm, ABC rejection sampling, first proposes  from the prior π() and then only accepts those  whose simulated summary statistics fall within an ϵ-ball (^o) = {:  - ^o ≤ϵ} around ^o. While simple and robust, the algorithm suffers from a poor trade-off between accuracy and computational efficiency <cit.>. Two lines of work aim at improving the trade-off from different perspectives: regression ABC and sequential ABC. 



 §.§ Regression ABC

This line of work first allows a much larger ϵ-ball in rejection ABC, and then uses the accepted samples to estimate the relationship between  and  through regression. The learned regression model is then used to adjust the samples so that, if the regression model is correct, the adjusted samples correspond to samples from ( | ^o) in (<ref>) with ϵ = 0.

Conditional on the summary statistics , an additive noise model is typically used to describe the relationship between  and  as

    | =   𝔼[|] _mean    +    |  _residual.

Here, we used the notation | and | to highlight that the above relationship is conditional on . In particular, the distribution of the residual may depend on the value of  so that p( | ) ≠ p(). Under the above model, to sample from the posterior ( | ^o), we need to evaluate 𝔼[|^o] and sample from the conditional distribution of  at ^o.

To learn 𝔼[|], regression ABC fits a regression function g(·) on the samples (^(i), ^(i)), i=1, …, I, collected in the preliminary run of rejection ABC. The fitting is typically done by the method of least squares, i.e. by minimising

    J(g) = 1/I∑_i=1^I  ^(i) - g(^(i)) ^2_2

since it is simple and, for sufficiently many samples (training data), the optimal g() is guaranteed to approach 𝔼[|]. A critical aspect of this approach is the choice of the function family among which one searches for the optimal g. For small amounts of training data, linear regression <cit.> is a safe choice while more flexible neural networks <cit.> may be used when more training data are available.

On obtaining g(·) we get an approximation for the posterior mean at ^o via g(^o). What remains unknown is the distribution p(|^o). To estimate this, regression ABC makes the homogeneity assumption that the distribution of the residual remains roughly unchanged for all  in the ϵ-ball (^o) around ^o:

    p(|) ≈ p()     ∀∈(^o)

With this assumption, the empirical residuals ^(i) - g(^(i)) are used as samples from p( | ^o), and the adjustment

    '^(i) = g(^o)_est.mean +   ^(i) - g(^(i))_est.residual

yields samples from ( | ^o) in (<ref>) with ϵ = 0. However, this result hinges on (a) the mean function being learned well and (b) the homogeneity assumption holding so that the estimated residuals do correspond to samples from p(|^o).
For a theoretical analysis of regression ABC, see e.g <cit.>.



 §.§ Sequential ABC

This second line of work aims at accelerating the inference by learning the posterior in an iterative way <cit.>. Using the prior () as proposal distribution in the generation of the parameter-data tuples (^(i), ^(i)), the methods first learn a rough but computationally cheaper approximation ( | ^o) to the posterior, e.g. by using a larger ϵ value <cit.>, or a small fraction of the overall simulation budget <cit.>. The learned posterior ( | ^o) is then used as a new proposal distribution (new 'prior') () in the next round. From newly generated tuples (^(i), ^(i)), with ^(i)∼(), the methods then construct an auxiliary posterior 
(|^o). In order to correct for using another proposal distribution than the prior, the auxiliary posterior (|^o) and samples from it need to be reweighted by

    w() ∝π()/()

see e.g. <cit.>. A reason why this method accelerates inference is that the new, more informed, proposal () generates more often summary statistics  that fall inside the ϵ-ball (^o). For a given ϵ, this then results in a higher acceptance rate and a reduction of the simulation cost. 

The concrete algorithms of sequential ABC may differ across methods. Earlier sequential Monte Carlo ABC methods <cit.> adopt a non-parametric posterior modeling, whereas the recent sequential neural posterior ABC ones <cit.> take a parametric approach using mixture density networks. However, the underlying basic principles remain roughly the same. 





§ METHODOLOGY

We here explain the proposed approach that makes use of the basic principles of both regression ABC and sequential ABC.

The performance of regression ABC is likely to suffer if the homogeneity assumption in (<ref>) does not hold. If ABC is run with a broad proposal distribution (e.g. the prior), only very few of the simulated  are close to ^o, so that deviation from the homogeneity is generally likely (Figure <ref>.a). However, if we use a more informed proposal distribution, more samples of  around ^o can be collected (Figure <ref>.b). In other words, the shape of the proposal distribution strongly affects whether the homogeneity assumption in regression ABC reasonably holds. This motivates us to take a two-stage approach where, in a first coarse-grained phase, we learn a rough approximation to the posterior and then, in a subsequent fine-grained-phase, use it as proposal distribution in order to generate a sufficient amount of (, ) pairs for which the homogeneity assumption in regression ABC is better satisfied. 

Regression function. We will be using regression-ABC in both the coarse- and fine-grained phase, for which we need to choose a family of regression functions g. Simple regression models with linear functions g are easy to learn but are not able to capture possible nonlinear relationships between  and . On the other hand, neural networks are more flexible; however they also require more training data. Here, we thus perform model selection by comparing the validation error between linear regression and a neural network trained with early stopping. The split into training/validation set is done on a 80/20 basis. 
Before regression, each  is preprocessed such that they have the same scale along each dimension.

Coarse-grained phase. In this stage we learn a rough but robust approximation ( | ^o) to the posterior using only a small fraction λ of the overall simulation budget N. We do this by performing ordinary regression ABC based on parameters generated from the prior (). We then pick a subset of the adjusted samples {'^(1),...,  '^(m)} whose corresponding summary statistics  are among the top m closest to ^o and model the adjusted samples with a Gaussian distribution so that

    ( | ^o) = 𝒩(; g(^o), ).

Here, g(·) is the learned regression function and  is the inflated sample covariance matrix computed as = α·1/m∑_i^m  ['^(i) - g(^o) ]['^(i) - g(^o) ]^⊤, α=1.5. By using a Gaussian, we work with the maximum entropy distribution given mean and covariance of the adjusted samples. The motivation for inflating the covariance matrix is that a slightly broader proposal distribution is more robust to e.g. estimation errors in the learned regression function (Figure <ref>.c).

For the value of m, we choose m so that the top 20% summary statistics closest to ^o are retained, which is a conservative but robust strategy given the small simulation budget in this phase. The whole coarse-grained learning procedure is summarized in Algorithm 1, where, for simplicity, we incorporated the reduction to summary statistics directly into the data generating process (simulator). 

Fine-grained phase. After learning ( | ^o) we use it as proposal distribution () in another round of regression ABC with the remaining  1-λ simulation budget. We adjust the top n samples whose summary statistics  are closest to ^o and then model their distribution as a Gaussian copula. 

Gaussian copula are a generalisation of the Gaussian distribution. They can be seen to model (the adjusted) ' = (θ'_1, …, θ'_K)^⊤ in terms of latent variables  that are element-wise non-linearly transformed to give θ'_k = F_k^-1(Φ(z_k)) where F_k is the marginal cumulative distribution function (CDF) of θ'_k and ∼𝒩(; 0, ). Under this construction, the auxiliary posterior ( | ^o) of the adjusted samples equals

    ( | ^o) = 𝒩(; 0, )/∏𝒩(z_k; 0, 1)_c_GC(u_1,...,u_k; Λ)·[ ∏^K_k f_k(θ_k) ],

where  is the correlation matrix whose diagonal elements are all ones and c_GC(·) is called the Gaussian copula density. 

The learning of Gaussian copula can be done in a semi-parametric way. Denote the adjusted samples by {'^(1),...,  '^(n)}:

    
  * Marginal distribution. Each f_k(θ_k) is learned by kernel density estimation (KDE), which is well-known to be typically accurate in 1D <cit.>.
    
  * Correlation matrix. The matrix  is learned by converting each sample ^'(j) to its corresponding latent representation ^(j) as z^(j)_k = Φ^-1(r_k(θ^'(j)_k)) where  r_k(x) = 1/n∑^n_j=11[θ^'(j)_k < x], followed by estimating  as = 1/n∑_j^n ^(j)^(j)⊤.


As for sequential ABC in (<ref>), we need to convert the auxiliary posterior (|^o) to an estimate of the actual posterior ( | ^o) by reweighing it by the ratio of the prior to proposal distribution. This gives

    ( | ^o) ∝()/( | ^o)( | ^o),

where the proposal distribution ( | ^o) was defined in (<ref>) and ( | ^o) was defined in (<ref>). Note that while the auxiliary posterior (|^o) is a Gaussian copula, the re-weighted ( | ^o) is generally not (because of e.g. the presence of the prior). 

Unlike in the coarse-grained phase, in this phase, we select n such that the top 2000 summary statistics closest to ^o are retained. We adopt this more aggressive strategy because here we have a larger simulation budget and better satisfy the homogeneity condition. We are therefore safe to keep a larger number of samples so as to reduce the Monte Carlo error. In addition, by retaining a fixed number of samples we effectively shrink the value of ϵ when the overall simulation budget N is increased. The whole fine-grained phase is summarized in Algorithm 2. 







    

    




























    

    



    






































 


§ COMPARISONS TO OTHER METHODS

Here we compare AGC-ABC with three methods in the field: neural network regression ABC <cit.>, Gaussian copula ABC <cit.>, and sequential neural posterior ABC <cit.>. 

AGC-ABC vs. NN-ABC. Both AGC-ABC and NN-ABC address the homogeneity issue (<ref>) in regression ABC. NN-ABC aims to model possible heterogeneity of the residuals by taking into account the possibly changing scale of the residuals. In contrast to this, AGC-ABC uses a two-stage approach to stay in the homogeneous regime, where not only the scale but also the dependency structure of the residuals is modelled.

AGC-ABC vs. GC-ABC. Both AGC-ABC and GC-ABC make use of Gaussian copula in posterior modeling. GC-ABC can be viewed as a special case of AGC-ABC without the adapting coarse-grained stage. The posterior in GC-ABC is a Gaussian copula whereas ours is generally not.





AGC-ABC vs. SNP-ABC. Both AGC-ABC and SNP-ABC learn the posterior in a coarse-to-fine fashion. The difference lies in the model complexity. SNP-ABC models the relationship between  and  with a mixture of Gaussians whose parameters are computed by a neural network. While such constructions can approximate the posterior more and more accurately as the number of mixture components increases, the number of neural network parameters is generally much larger than in AGC-ABC where we only need to estimate the posterior mean and a correlation matrix. Additionally, the marginal distributions in our copula-based method can generally be modelled more accurately than in SNP-ABC, since we learn them nonparametrically for each dimension. 























§ EXPERIMENTS



 §.§ Setup

We compare the proposed adaptive Gaussian Copula ABC (AGC-ABC) with five ABC methods in the field: rejection ABC <cit.>, regression ABC <cit.>, GC-ABC <cit.>, NN-ABC <cit.> and SNP-ABC <cit.>. Comparisons are done by computing the JSD between the true and the approximate posterior, and when the true posterior is not available, we approximate it by running a ABC rejection sampling algorithm with (a) well chosen summary statistics  and (b) an extremely small ϵ value e.g the 0.00001 quantile of the population of - ^o, and then estimate it by KDE. Similarly, when the posterior approximated in each method is not available in analytical form (e.g for REJ-ABC, REG-ABC and NN-ABC), we estimate it by KDE. For two distributions P and Q, the JSD is given by

    JSD(P, Q)    = 1/2KL[P P+Q/2] + 1/2KL[Q P+Q/2].





The calculation of JSD is done numerically by a Riemann sum over 30^K equally spaced grid points with K being the dimensionality of . The region of the grids is defined by the minimal and maximal values of the samples from P and Q jointly. Before the JSD calculation, all distributions are further re-normalized so that they sum to unity over the grid.

To make the comparison sensible, each method has the same total number of simulations (budget) N available. Note, however, that the different methods may use the simulation budget in different ways:

    
  * For rejection ABC we pick the top 2000 samples closest to ^o for posterior estimation;
    
  * For the two linear regression-based methods, regression ABC and GC-ABC, we pick the top 2000 samples for posterior learning; 
    
  * For NN-ABC and SNP-ABC we use all samples for learning, in line with their goal of modelling the relationship between  and  without restriction to an ϵ-ball. 
    
  * For AGC-ABC, as mentioned before, we use the top 20% of the samples in the coarse-grained phase while the top 2000 samples in the fine-grained phase. 

Both AGC-ABC and the SNP-ABC learn the posteriors in a coarse-to-fine fashion. For these two methods, we assign 20%/80% of the overall simulation budget N to the coarse/fine phase respectively. Following the original literature <cit.>, the neural networks in SNP-ABC have one/two hidden layers to model one/eight mixture of Gaussian in the two phases respectively, with 50 tanh units in each layer. For AGC-ABC, the neural network, which is considered in model selection, has two hidden layers with 128 and 16 sigmoid units each. This is also the same network used in NN-ABC. All neural networks in the experiments are trained by Adam <cit.> with its default setting.








 §.§ Results




  
A toy problem The first problem we study is a 2D Gaussian copula model:

    (|)    = c_CG(
             u_1, u_2;
             [   1,  θ_3; θ_3,    1;      ])
             ∏_k=1^2 f_k(x_k;θ_k) 
     
        f_1(x_1;θ_1)    = Beta(x_1; θ_1, 2) 
    
        f_2(x_2;θ_2)    = θ_2 𝒩(x_2; 1, 1) + (1-θ_2) 𝒩(x_2; 4, 0.25)

where u_k=F_k(x_k) is the value of the marginal CDF at x_k and c_CG(·) is the aforementioned Gaussian copula density. The parameters of interest are = (θ_1, θ_2, θ_3 ) and the true parameters are ^* = (6.0, 0.5, 0.6). We place a flat, independent prior on : θ_1 ∼𝒰(0.5, 12.5), θ_2 ∼𝒰(0, 1), θ_3 ∼𝒰(0.4, 0.8). The summary statistics  are taken to be the union of (a) the 20 equally-spaced marginal quantiles and (b) the correlation between z_1 = Φ^-1(u_1), z_2 = Φ^-1(u_2) of the generated data = {^(1),…, ^(200)}. In this problem, the true posterior is available analytically.

Figure <ref> shows the JSD between the true and estimation posterior for the different methods (on log-scale, vertical lines indicate standard errors, each JSD is obtained by calculating the average of 15 runs for different observed data, the results shown in the figures below have the same setup). Expectedly given the model specification, AGC-ABC yields a lower JSD than the other methods for all simulation budgets. We further see that there is only a small difference in performance between AGC-ABC and NN-ABC, especially in large budget settings. A possible explanation is that the distribution of the residuals p(|) does only weakly depend on . To further investigate this, we compared the distribution of the residuals at ^o with the distribution p_ϵ() of the residuals within ϵ-balls of increasing radius around ^o. 
Table 1 shows the JSD values between the two distributions for different values of ϵ.[For ϵ = 50%, the JSD is 0.061.] It can be seen that the JSD between p(|^o) and p_ϵ() increases as ϵ increases but the changes are not large. This is in line with the small gain of AGC-ABC over NN-ABC. The supplementary material further compares the contour plots of p(|^o) and p_ϵ(). 





  
The M/G/1 queueing problem. 
The M/G/1 queueing model is a real-world problem that describes the processing of incoming jobs in a single server system. Jobs in this model are served on a first-in-first-out (FIFO) basis. The whole data generating process in M/G/1 can be described as:

    s_i     ∼𝒰(θ_1, θ_1 + θ_2) 
    
        v_i - v_i-1   ∼Exp(θ_3) 
    
        d_i - d_i-1   = s_i + max(0, v_i - d_i-1)

where s_i, v_i, d_i are the service time, visiting time and finish time of job i respectively. The parameters of interest in this problem are = ( θ_1, θ_2, θ_3 ) and the true parameter values are ^* = ( 1.0, 4.0, 0.2 ). Again, we place a flat, independent prior on : θ_1 ∼𝒰(0, 10), θ_2 ∼𝒰(2, 6), θ_3 ∼𝒰(0, 0.33). The observed data  in this problem are given by the time intervals x_i = d_i - d_i-1 between the finish time of two consecutive jobs. 200 such intervals are observed so that = {x_1, …, x_200}. The summary statistics ^o are taken as the 16-equally spaced quantiles of . These statistics are further preprocessed to have zero mean and unit variance on each dimension. 

Figure <ref> shows the JSD of the different methods under various simulation budgets. We see that there is a larger gap between AGC-ABC and traditional regression ABC methods for small simulation budgets. This is because the underlying residual strongly depends on  in this problem; it is highly heterogeneous (see Table <ref> and the contour plots in the supplementary material). This also points to why NN-ABC does here not perform as good as AGC-ABC — only modelling the dependency between the summary statistics and the scale (variance) of the residuals does not seem enough. SNP-ABC achieves very good performance compared to traditional ABC methods due to its heterogeneous modelling but it is still less accurate than AGC-ABC, which is natural since it requires more data to fit.









  
MA(2) time series problem The third problem we consider is the second order moving average model in time series analysis, known as MA(2). In this model, data are generated as follows

    x^(t)   = w^(t) + θ_1 w^(t-1) + θ_2 w^(t-2),

where x^(t) is the observation at time t and w^(t) is some unobservable standard normal noise. The parameters of interest are = ( θ_1, θ_2 ) and the true parameters are  ^* = (0.6, 0.2). As in the previous problems, we place a flat and uniform prior on : θ_1 ∼𝒰(0, 1), θ_2 ∼𝒰(0, 1). Following the ABC literature <cit.>, we adopt the autocovariance with lag 1 and lag 2 as the summary statistics. These statistics are computed from a time series = (x^(1), …, x^(200)) of length 200. The statistics are further preprocessed to have zero mean and unit variance for each dimension. 



Figure <ref> reports the JSD of the different methods under various simulation budgets. While AGC-ABC performs best for small budgets, the simpler REG-ABC method performs better in the larger budget case. This might be due to Gaussian copula being inadequate for the modelling of the residuals. However, the difference between AGC-ABC and REG-ABC is not large as evidenced by the contour plots in the supplementary materials. We also see that the 8-mixture of Gaussian construction in SNP-ABC might be not powerful enough for this problem, so that it incurs a bias which is not
eliminated as the simulation budget increases.





  
Lotka-Volterra problem The last problem we consider is a stochastic dynamical system in biology that describes predator-prey dynamics. There are four possible events in the system: (a) a predator being born, (b) a predator dying, (c) a prey being born, (d) a prey being eaten by the predator. The events can be described probabilistically as

    (X → X+1 | )    ∝exp(θ_1) · XY,   
    (X → X-1 | )    ∝exp(θ_2) · X, 
    (Y → Y+1 | )    ∝exp(θ_3) · Y,  
    (Y → Y-1 | )    ∝exp(θ_1) · XY,

where X, Y are the numbers of predators and preys respectively. For a more rigorous formulation see e.g. <cit.>. The parameters of interest are = ( θ_1, θ_2, θ_3 ) and the true parameters are ^* = ( log(0.01), log(0.5), log(0.1) ). The initial sizes of the populations are ( X^0, Y^0 ) =  (50, 100). While the likelihood is intractable, sampling from the model is possible <cit.>. We place a flat and independent prior on :
θ_1 ∼𝒰(-5, -1), θ_2 ∼𝒰(-1, 1), θ_3 ∼𝒰(-1, 1). We simulate the system for a total of 8 time units and record the values of X and Y after every 0.2 time units. This yields two time series of length 41 being our observed data . The summary statistics are  taken as the records between every two consecutive time points.

Figure <ref> shows the JSD of the different methods. The results show that AGC-ABC offers a better efficiency-accuracy trade-off than the other methods, with a clear performance gain for small simulation budgets. This is due to the fact that the underlying distribution p(|) is highly heterogeneous (see Table <ref> and the the contour plots of p(|) in the supplementary material).












§ CONCLUSION

We considered the problem of posterior density estimation when the likelihood is intractable but sampling from the model is possible. We proposed a new method for approximate Bayesian computation (ABC) that combines the basic ideas from two different lines of ABC research, namely regression ABC and sequential ABC. We found that the resulting algorithm strikes a good trade-off between accuracy and computational cost, being particularly effective in the regime of smaller simulation budgets. 

The motivation behind the proposed algorithm was the homogeneity assumption on the residuals that is required for regression ABC to work well. The proposed method takes a sequential approach by first generating training data so that the homogeneity assumption is better satisfied, and then models the data with the aid of a Gaussian copula and existing techniques from regression ABC. The method — adaptive Gaussian copula ABC —  can thus either be viewed as an adaptive version of classical regression ABC methods <cit.>, or a computationally cheaper version of recent sequential neural posterior approaches <cit.>.

While Gaussian copula are powerful they are not silver bullets. Extending the Gaussian copula model to other statistical models such as mixture of Gaussian copula <cit.> or vine copulas <cit.> may be a future research direction worth exploring. 

 



§ ACKNOWLEDGEMENT

We would like to thank the anonymous reviewers for their insightful and fruitful suggestions. 











§ CONTOUR PLOTS OF THE RESIDUALS IN EACH PROBLEM






 §.§ Gaussian copula toy problem












 §.§ M/G/1 problem









 §.§ MA(2) problem







 §.§ Lotka-Volterra problem









§ CONTOUR PLOTS OF APPROXIMATE POSTERIOR LEARNED BY EACH METHOD IN MA(2)









