Local Kernel Dimension Reduction in Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) is a popular sampling method in applications involving intractable likelihood functions. Without evaluating the likelihood function, ABC approximates the posterior distribution by the set of accepted samples which are simulated with parameters drawn from the prior distribution, where acceptance is determined by the distance between the summary statistics of the sample and the observation. The sufficiency and dimensionality of the summary statistics play a central role in the application of ABC. This paper proposes Local Gradient Kernel Dimension Reduction (LGKDR) to construct low dimensional summary statistics for ABC. The proposed method identifies a sufficient subspace of the original summary statistics by implicitly considers all nonlinear transforms therein, and a weighting kernel is used for the concentration of the projections. No strong assumptions are made on the marginal distributions nor the regression model, permitting usage in a wide range of applications. Experiments are done with both simple rejection ABC and sequential Monte Carlo ABC methods. Results are reported as competitive in the former and substantially better in the latter cases in which Monte Carlo errors are compressed as much as possible.


Introduction
Monte Carlo methods are popular tools in sampling and inference problems.
While the Markov Chain Monte Carlo methods find successes in applications where likelihood functions are known up to a unknown constant, MCMC can not be used in scenarios where likelihoods are intractable.For these cases, if the problem can be characterized by a generative model, Approximate Bayesian Computation (ABC) is often a candidate approach.ABC is a Monte Carlo method that approximates the posterior distribution by jointly generating simulated data and parameters and does the sampling based on the distance between the simulated data and the observation, without evaluating the likelihoods.ABC was first introduced in population genetics [1] [2] and then been introduced to a range of complex applications including dynamical systems [3], ecology [4], Gibbs random fields [5] and demography [6].
The accuracy of ABC posterior depends on the Monte Carlo errors in-duced in the samplings.Given the generative model p(y|θ) of observation y obs with parameter θ, consider summary statistics s obs = G s (y obs ) and s = G s (y), where G s : Y → S is the mapping from the original sample space Y to low dimensional summary statistics S. The posterior distribution, p(θ|y obs ), is approximated by p(θ|s obs ), which is constructed as p(θ|y obs ) = p ABC (θ, s|s obs )ds, with p ABC (θ, s|s obs ) ∝ p(θ)p(s|θ)K( s − s obs /ǫ), where K is a smoothing kernel with bandwidth ǫ.In the case of simple rejection ABC, K is often chosen as an indicator function I( s − s obs < ǫ).
If the summary statistics s are sufficient, it can be shown that (1) would converges to the posterior p(θ|s obs ) as ǫ goes to zero [7].
As can be seen above, sampling is based on the distance between the summary statistics of the simulated sample s and the observation s obs .Approximation errors are induced by the distance measure and proportional to the distance threshold ǫ.It is desirable to set ǫ as small as possible, but a small threshold will increase the simulation time.This is a trade-off between the accuracy and the efficiency (simulation time) determined by the choices of thresholds and summary statistics.According to recent results on asymptotic properties of ABC [?] [?], assuming that the summary statistics follow the central limit theorem, the convergence rate of ABC when accepted sample size N → ∞ is depended on the behavior of µ = ǫd N , where ǫ is the threshold above and the d N is defined as of same magnitude of eigen(Σ N ), the eigenvalue of the covariance matrix of the summary statistics as the function of N. In practice, if a specific sampling method is chosen, the threshold ǫ is constrained by the computing resources and time, thus can be accordingly determined.The design of summary statistics then remains the most versatile and difficult part in developing an efficient ABC algorithm.To avoid the "curse of dimensionality", summary statistics should be also chosen to be low dimensional in addition of sufficiency.
A vast body of literature of ABC have been published.Many are devoted to reduce the sampling error by using more advanced sampling method, from simple Rejection method [8], Markov Chain Monte Carlo(MCMC) [9] to more sophisticated methods like sequential Monte Carlo [10][3] and adaptive sequential Monte Carlo methods [?].
In this paper, we focus on the problem of summary statistics.In early works of ABC, summary statistics are chosen by domain experts in an adhoc manner.It is manageable if the dimensionality is small; which usually means that the model is well understood by the experts.But choosing a set of appropriate summary statistics is much more difficult in complex models.To address this problem, a set of redundant and hopefully sufficient summary statistics are often constructed in the first place, as initial summary statistics; dimension reduction methods are then applied yielding a set of low dimensional summary statistics while persevering the sufficiency or information.
Many dimension reduction methods have been proposed for ABC.Entropy based subset selection [11], partial least square [12], neural network [13] and expected posterior mean [14] are a few of them.The entropy based subset selection method works well in instances where the set of low dimensional summary statistics is a subset of the initial summary statistics, but the computational complexity increases exponentially with the size of the initial summary statistics.The partial least square and neural network methods aim to capture the nonlinear relationships of the original summary statistics.
In both cases, a specific form of the regression function is assumed.A comprehensive review [15] discusses the methods mentioned above and compares the performances.While the results are a mixed bag, it is reported that the expected posterior mean method (Semi-automatic ABC) [14] produces relatively better results compared to the methods mentioned above in various experiments.It is a popular choice also due to its simplicity.
Semi-automatic ABC [14] uses the estimated posterior mean as summary statistics.A pilot run of ABC is first conducted to identify the regions of parameter space with non-negligible probability mass.The posterior mean is then estimated using the simulated data from that region and is used as the summary statistics in a following formal run of ABC.A linear model of the form: θ i = β (i) f (y) + ǫ i is used in the estimation, where f (y) are the possibly non-linear transforms of the data.For each application, the features f (y) are carefully designed to achieve a good estimation.In practice, it may be difficult to determine a good set of features given a particular application.
In these cases, a vector of powers of the data (y, y 2 , y 3 , y 4 , ...) is often used as noted in [14].
To provide a principled way to design the regression function, capture the higher order non-linearity and realize an automatic construction of summary statistics, we introduce the kernel based sufficient dimension reduction method as an extension of the linear projection based Semi-automatic ABC.
The dimension reduction used here is a localized version of Gradient based kernel dimension reduction (GKDR) [18].GKDR, as a general dimension reduction method, estimates the projection matrix onto the sufficient low dimensional subspace by extracting the eigenvectors of the kernel derivatives matrices in the reproducing kernel Hilbert spaces (RKHS).We give a brief review of the method in Section 2. In addition to the GKDR, in which the estimation averages over all data points to reduce variance, a localized GKDR is proposed by averaging over a small neighborhood around the observation in ABC.Each point is weighted using a distance metric measuring the difference between the simulated data and the observation, similar to the distance kernel function in (1), to concentrate on the observation point.
Another proposal is to use different summary statistics for different parameters.Note that sufficient subspace for different parameters can be different, depending on the particular problem.In these cases, as will be shown later by experiments, applying separated dimension reduction procedure yields a better estimation.
The proposed method gives competitive results in comparison with Semi-automatic ABC [14] when using simple rejection sampling.Substantial improvements are reported in the sequential Monte Carlo cases, where threshold ǫ is pushed to as small as possible to isolate the performance of summary statistics from the Monte Carlo error.
The paper is organized as follows.In Section 2, we review GKDR and introduce its localized modification followed by discussions of computation considerations.In Section 3, we show simulation results for various commonly conducted ABC experiments, and compare the proposed method with the Semi-automatic ABC.

Local Kernel Dimension Reduction
In this section, we review the Gradient based Kernel Dimension Reduction (GKDR) and propose the modified Local GKDR (LGKDR).Discussions are given at the end of this section.

Gradient based kernel Dimension Reduction
Given observation (s, θ), where s ∈ R m are initial summary statistics and θ ∈ R is the parameter to be estimated in a specific ABC application.Assuming that there is a where B = (β 1 , ..., β d ) is the orthogonal projection matrix from R m to R d . The columns of B spans U and B T B = I d .Condition (2) shows that given B T s, θ is independent of the initial summary statistics s.It is then sufficient to use d dimensional constructed vector z = B T s as the summary statistics.This subspace U is called effective dimension reduction(EDR) space [17] in classical dimension reduction literatures.While there are a tremendous amount of published works about estimating the EDR space, in this paper, we propose to use GKDR in which no strong assumption of marginal distribution or variable type is made.The following is a brief review of GKDR, and for further details, we refer to [18] [19] [16].
Let B = (β 1 , ..., β d ) ∈ R m×d be the projection matrix to be estimated, and z = B T s.We assume ( 2) is true and p(θ|s) = p(θ|z).The gradient of the regression function is denoted by ∇ s as which shows that the gradients are contained in the EDR space.Given the The projection directions β lie in the subspace spanned by the eigenvectors of M. It is then possible to estimate the projection directions using eigenvalue decomposition.In GKDR, the matrix M is estimated by the kernel method described below.
Let Ω be an non-empty set, a real valued kernel k : x j ) ≥ 0 for any x i ∈ Ω and c i ∈ R. Given a positive definite kernel k, there exists a unique reproducing kernel Hilbert space (RKHS) H associated with it such that: (1)k(•, x) spans H;(2)H has the reproducing property [20]: ) be Gaussian kernels defined on R m and R, associated with RKHS H S and H Θ , respectively.With assumptions of boundedness of the conditional expectation E(θ|S = s) and the average gradient functional with respect to z, the functional can be estimated using cross-covariance operators defined in RKHS and the consistency of their empirical estimators are guaranteed [16].Using these estimators, we construct a covariance matrix of average gradients as where G S and G Θ are Gram matrices k S (s i , s j ) and k Θ (θ i , θ j ), respectively.
∇k S ∈ R n×m is the derivative of the kernel k S (•, s i ) with respect to s i , and ǫ n is a regularization coefficient.This matrix can be viewed as the straight forward extension of covariance matrix in principle component analysis (PCA); the data here are the features in RKHS representing the gradients instead of the gradients in their original real space.
The averaged estimator M = 1/n n i=1 M n (s i ) is calculated over the training sample (s 1 , θ 1 ), ..., (s n , θ n ).Finally, the projection matrix B is estimated by taking d eigenvectors corresponding to the d largest eigenvalues of M just like in PCA, where d is the dimension of the estimated subspace.

Local Modifications
As discussed above, the estimator M is obtained by averaging over the training sample s i .When applied to ABC, since only one observation sample is available, we propose to generate a set of training data using the generating model and introduce a weighting mechanism to concentrate on the local region around the observation and avoid regions with low probability density.
Given simulated data X 1 , ..., X N and a weight kernel K w : R m → R, we propose the local GKDR estimator where M is m × m matrix and K w (X i ) is the corresponding weight.K w (x) can be any weighting kernel.In the numerical experiments, a triweight kernel is used, which is written as where 1 u<1 is the indicator function, and X th is the threshold value which determines the bandwidth.The normalization term of the triweight kernel is omitted since it does not change the eigenvectors we are estimating.The bandwidth determined by X th is chosen by empirical experiments and will be described in 2.4.The Triweight kernel is chosen for its concentration in the central area than other "bell shaped" kernels and works well in our experiments.Other distance metrics could be used instead of squared distance.
Description of LGKDR algorithm are given in Algorithms 1. Procedure GenerateSample is the algorithm to generate sample with parameter as input.Procedure LGKDR is the algorithm to calculate matrix M(X i ) as given in ( 4) and ( 5).
Since the dimension reduction procedure is done before the sampling, it works as a pre-processing unit to the main ABC sampling procedure.It can be embodied in any ABC algorithm using different sampling algorithms.In this paper, the rejection sampling method is firstly employed for its simplicity and low computation complexity as a baseline.Further results on Sequential Monte Carlo ABC are also reported to illustrate the advantage of the purposed method.In these experiments, the distance thresholds are pushed to as small as possible to suppress the Monte Carlo errors and isolate the effects of summary statistics alone.

Separated Dimension Reduction
It is expected that separated construction of summary statistics for each parameter is beneficial for applications.Different information may be crucial for different parameters.If we estimate the projection directions separately for each parameter, the accuracy of the projection can be improved, and a input : weighting kernel K w , procedure GenerateSample, prior distribution D prior , number of accepted sample N, process LGKDR output: projection matrix B training sample generation; LGKDR input : projection matrix B, distance kernel K d , bandwidth ǫ, number of sample N ABC and observation X ob output: set of parameters {θ(j)} Algorithm 2: Rejection-ABC input : projection matrix B, distance kernel K d , target threshold ǫ t , number of particle N ABC , effective sample size threshold ess t output: set of parameters {θ(j)} Algorithm 3: Sequential-ABC lower dimensionality may be achievable.
The LGKDR incorporates information of θ in the calculation of gradient matrix M .If θ is a vector, the relation of different elements of θ are contained in the gram matrix G θ as in (4).Separate estimations concentrate on the information of the specific parameter rather than the whole vector.
As shown in the experiments in Section 3.2, it can construct significantly more informative summary statistics in some problems by means of reducing estimation error.
For Semi-automatic ABC [14], the summary statistics for each parameter is the estimated posterior mean, thus naturally separated.Using separated summary statistics means using a 1 dimensional summary statistics for each parameter.our experiments do not show good results using this setting.
For best subset selection methods [12][21], summary statistics are chosen as the best subset of the original summary statistics using mutual information or sufficiency criterion.It can also be extended to a separated selection procedure.In LGKDR, we simply construct different summary statistic by using only the particular parameter as the response variable.

Discussion on hyper parameters
Several parameters need to be set for achieving good results in LGKDR ABC.
Parameters for the sampling procedures will be discussed in the experiments section.In this section, the parameters of the LGKDR part is explained.
First, the bandwidth of the weighting kernel, which measures the degree of concentration of weightings of projection directions, affects the accuracy of LGKDR.By selecting a large bandwidth, the weights of directions spread out a larger region around the observation points.A small bandwidth concentrates the weights on the directions estimated close to the observation sample.In our experiments, a bandwidth corresponding to an acceptance rate of approximately 10% gives a good result and is used throughout the experiments.The same parameter is set for the Semi-automatic ABC as well for the similar purpose.A more principled method for choosing bandwidth, like cross validation, could be applied to select the acceptance rate if the corresponding computation complexity is affordable.
The bandwidth parameters of the Gaussian kernels σ S , σ Θ and the regularization parameter ǫ n are the configuring parameters that are crucial to all kernel based methods.The first two determine the function spaces associated with the positive definite kernels and the regularization parameter affects the convergence rate (see [?]).In our method, cross validation is adopted to select the proper parameters for different experiments.In the cross validation, for each set of candidate parameters, the summary statistics are constructed using a simulated observation θ obs , s obs , a training set (θ training , S training ) and a test set (θ test , S test ).A small pilot run of rejection ABC is performed and the estimation of parameters are calculated by kNN regression of θ test with the S test .K is set to 5 in all cases.The set of parameters that yield the smallest least error between the θ test and θ obs are chosen.The final summary statistics are constructed using these chosen parameters and are passed to the formal run of ABC subsequent to the summary construction procedure.

Computational Complexity
Computational complexity is one central concern of ABC methods.
LGKDR is more computationally demanding than the linear regression-type methods.
It requires matrix inversion and solving eigenvalue problems and above all, the cross validation procedure.The actual complexity depends on the training sample size used.For the experiments shown in this paper, the training sample size are fixed to 2×10 3 and 10 4 for LGKDR and Semi-automatic ABC, respectively.Under this setting, the total computational time of LGKDR are about 10 times over the linear regression.While the computational complexity is higher, it is a necessary price to pay if the non-linearity between the initial summary statistics is strong.For these cases, being unable to capture the non-linear information in summary statistics would induce a poor sampling performance, which leads to a biased estimation.Meanwhile, if the generating model itself is complex, the computational time used on the LGKDR will become less significant.Finally, although the cross validation procedure takes the majority of computation time in LGKDR, it needs to be performed only once for each problem to fix the parameters.Once the parameters are chosen, the computation of LGKDR is comparable to the linear-type algorithms.Overall the computational complexity of an ABC method depends on both the summary statistics calculation step and the following sampling step.For complex models like ones in population genetics, sampling is significantly more time consuming than the dimension reduction procedure whichever we use.

Experiments
In this section, we investigate three problems to demonstrate the performance of LGKDR.Our method is compared to the classical ABC using initial summary statistics and the Semi-automatic ABC [14] using estimated posterior means.In the first problem, we discuss a population genetics model, which was investigated in many ABC literatures.We adopt the initial summary statistics used in [22], and rejection ABC is used as the sampling algorithm.
In the second problem, a M/G/1 stochastic queue model which was used in [13] and [14] are discussed.While the model is very simple, the likelihood function could not be trivially computed.In the last experiment we explore the Ricker model as discussed in [23] and [14].The latter two problems are investigated by both Rejection ABC and sequential Monte Carlo ABC method (SMC ABC) [?], the first problem is omitted from SMC ABC because it involves repeated calling an outside program for simulation and is too time consuming for SMC ABC.

Implementation Details
The Rejection ABC is described in Algorithm 2 and the SMC ABC is shown in Algorithm 3. The hyper-parameters used in LGKDR is set as discussed in section 2.4.We use a modified code from [?] and R package "Easyabc" [?] in our SMC implementation and would like to thank the corresponding authors.Gaussian kernels are used in all the LGKDR algorithms.The detailed specifications of Semi-automatic ABC will be described in each experiment.
For evaluation of the experiments conducted using rejection ABC, a set of parameters θ j where j ∈ 1, ..., N obs and corresponding observation sample Y j obs are simulated from the prior and the conditional probability p(Y |θ), respectively, and used as the observations.For each experiment, we fix the total number of simulations N, and the number of accepted sample N acc .The sample used for rejection are then generated and fixed for all three methods.
In this case we can ignore the randomness in the simulation program and accurately determine the acceptance rate, which is the single most influential parameter for estimation accuracy.The Mean squared error (MSE) over the accepted parameters θj i and observation θ j are defined as The Averaged Mean Square Error (AMSE) is then computed as the average over MSE j of each observation pair (θ j , Y j obs ) as It is used as the benchmark for Rejection ABC.Because of the difference of computation complexity, for fairness of comparison, the acceptance rates are set differently.For LGKDR, the acceptance rate is set to 1%; while for Semiautomatic ABC and original ABC, the acceptance rates are set to 0.1%.The training sample and simulated sample are generated from the same prior and remain fixed.
For SMC ABC, to get to as small tolerance as possible, the simulation time is different for different method.AMSE is used as benchmark for accuracy and computational time are reported for each experiment.In the case of Ricker model, due to the extremely long simulation time, only one observation is used and AMSE is just MSE in this case.

Parameter Settings
Several parameters are necessary in running the simulations in ABC.For Rejection ABC, the total number of samples N and the accepted number of samples N acc are set before the simulation as mentioned above.For Semiautomatic ABC and LGKDR, a training set needs to be simulated to calculate the projection matrix.For LGKDR, a further testing set is also generated for cross validation purposes.The value of these parameters are reported in the corresponding experiments.The simulation time for generating these sample set are negligible compared to the main ABC, especially in SMC ABC.For LGKDR, another important parameter is the target dimensionality D. There is no theoretically sound methods available to determine the intrinsic dimensionality of the initial summary statistics; yet in practice, since the projection matrix consists of extracted eigenvectors of the matrix M as in (5), it is straightforward to test the performance of LGKDR using different projection matrices B with different setting dimensionality and determine a proper dimensionality directly.In our experiments, we run several rejection ABC using different B, and fix the dimensionality.A starting point can be set by preserving 70% of the largest eigenvalues in magnitude and it usually works well.There are a large collection of literatures on how to choose the number of principle components in PCA, which is similar to our problem, for example, see [?] and reference therein.

Population Genetics
Analysis of population genetics is often based on the coalescent model [24].A constant population model is used in simple situations, where the population is assumed unchanged across generations.The parameter of interests in this case is the scaled mutation rate θ, which controls the probability of mutation between each generation.The detailed introduction of coalescent models can be found in [25].Various studies [26] [9] [10] have been conducted in population genetics following different sampling algorithms.In this study, we adopt the setting of kernel ABC [22] and compare the performance with ABC and Semi-automatic ABC.
100 chromosomes are sampled from a constant population (N = 10000).
The summary statistics are defined using the spectrum of the numbers of segregating sites, s sf s , which is a coarse-grained spectrum consisting of 7 bins based on the Sturges formula (1 + log 2 S seg ).The frequencies were binned as follows: 0 − 8%, 8 − 16%, 16 − 24%, 24 − 32%, 32 − 40%, 40 − 48% and 48 − 100%, we use the uniform distribution θ ∼ [0, 30] in this study rather than the log-normal distribution in [22].As ABC is often used for exploratory researches, we believe that the performance based on an uninformative prior is important for evaluating summary statistics.The program package ms is used to generate the sample, which is of common choice in literature of coalescent model [27].
We test 3 typical scaled mutation rates 5, 8 and 10 rather than random draws from the prior.The results are averaged over 3 tests.A total number of 10 6 sample is generated; 10 5 sample is generated as the training sample for LGKDR and Semi-automatic ABC.Different acceptance rates are set for different methods as discussed above.We use s sf s as the summary statistics for both Semi-automatic ABC and LGKDR.Local linear regression is used as the regression function for the former.In LGKDR, the dimension is set to 2.
As shown in Table-1, the performance of both LGKDR and Semi-automatic ABC improve over original ABC method.
LGKDR and Semi-automatic ABC achieve very similar results suggesting that the linear construction of summary statistics are sufficient for this particular experiment.

M/G/1 Queue Model
The M/G/1 model is a stochastic queuing model that follows the first-comefirst-serve principle.The arrival of customers follows a Poisson process with However a simulation model with parameter (λ, µ) can be easily implemented to simulate the model.It has been analyzed by ABC using various different dimension reduction methods as in [14] and [13], with comparison to the indirect inference method.We only compare our method with Semi-automatic ABC, since it produces substantially better results then the other methods mentioned above.
The generative model of the M/G/1 model is specified by where Y n is the inter-departure time, U n is the service time for the nth customer, and W i is the inter-arrival time.The service time is uniformly distributed in interval [θ 1 , θ 2 ].The inter-arrival time follows an exponential distribution with rate θ 3 .These configurations stay the same as [13] and [14].
provement compared to Semi-automatic ABC.Separated estimations for θ 2 and θ 3 give no improvements.It suggests that the sufficient dimension reduction subspace for θ 1 is different from the others and a separated estimation of θ 1 is necessary.
For SMC ABC, a set of 10 pairs of parameters are generated, and results on SMC and LGKDR are reported.all other setting are same as rejection ABC.We omit the results of using Semi-automatic ABC since the sequential chain did not converge properly using these summary statistics and the in- LGKDR summary statistics is included in the total simulation time and listed in the bracket.The results show that LGKDR gives better results of parameter θ 1 and θ 2 , using less time compared to SMC ABC with set E2.
The estimation of θ 3 is worse but the difference is small (0.005).Focusing on θ 3 produces an estimation as good as in SMC ABC.

Ricker Model
Chaotic ecological dynamical systems are difficult for inference due to its dynamic nature and the noises presented in both the observations and the process.Wood [23] addresses this problem using a synthetic likelihood inference method.Fearnhead [14] tackles the same problem with a similar setting using the Semi-automatic ABC and reports a substantial improvement over other methods.In this experiment, we adopt the same setting and apply LGKDR with various configurations.
A prototypic ecological model with Richer map is used as the generating model in this experiment.A time course of a population N t is described by where e t is the independent noise term with variance σ 2 e , and r is the growth rate parameter controlling the model dynamics.A Poisson observation y is made with mean φN t .The parameters to infer are θ = (log(r), σ 2 e , φ).The initial state is N 0 = 1 and observations are y 51 , y 52 , • • • , y 100 .
The original summary statistics used by Wood [23] are the observation mean ȳ, auto-covariances up to lag 5, coefficients of a cubic regression of the ordered difference y t − y t−1 on the observation sample, estimated coefficients for the model y 0.3 t+1 = β 1 y 0.3 t + β 2 y 0,6 t + ǫ t and the number of zero observations 100 t=51 1(y t = 0).This set is denoted as E0 as in [14].Additional two sets of summary statistics are defined for Semi-automatic ABC.The smaller E1 contains E0 and 100 t=51 1(y t = j) for 1 ≤ j ≤ 4, logarithm of sample variance, log( 100 t=51 y j t ) for 2 ≤ j ≤ 6 and auto-correlation to lag 5. Set E2 further includes time-ordered observation y t , magnitude-ordered observation y (t) , y 2 t , y 2 (t) , {log(1+y t )}, {log(1+y (t) )}, time difference ∆y t and magnitude difference ∆y (t) .Additional statistics are added to explicitly explore the non-linear relationships of the original summary statistics and are carefully designed.
In Rejection ABC, we use set E0 for ABC without dimension reduction since the dimension of the larger sets induces severely decreased performance.
Sets E1 and E2 are used for Semi-automatic ABC as in [14].In LGKDR, we tested sets E0 and E1 in different experiments.The result on E2 is omitted as the result is similar with using the smaller set of statistics, indicating that manually designed non-linear features are unnecessary for LGKDR.The The results are shown in Table 4.The performance of Semi-automatic ABC using the bigger set E2 is similar to ABC but is substantially worsen with set E1, suggesting that the non-linear information are essential for an accurate estimation in this model.These features are needed to be explicitly designed and incorporated into the regression function for Semi-automatic

Conclusions
We proposed the LGKDR algorithm for automatically constructing summary statistics in ABC.The proposed method assumes no explicit functional forms of the regression functions nor the marginal distributions, and implicitly incorporates higher order moments up to infinity.As long as the initial summary statistics are sufficient, our method can guarantee to find a sufficient subspace with low dimensionality.While the involved computation is more expensive than the simple linear regression used in Semi-automatic ABC, the dimension reduction is conducted as the pre-processing step and the cost may not be dominant in comparison with a computationally demanding sampling procedure during ABC.Another advantage of LGKDR is the avoidance of manually designed features; only initial summary statistics are required.With the parameter selected by the cross validation, construction of low dimensional summary statistics can be performed as in a black box.For complex models in which the initial summary statistics are hard to identify, LGKDR can be applied directly to the raw data and identify the sufficient subspace.We also confirm that construction of different summary statistics for different parameter improve the accuracy significantly.
duced error was too large to be meaningful.In SMC ABC, two experiments are reported: SMC ABC1 and SMC ABC2.The number of particles are set to 2 × 10 4 and 10 5 , respectively.In LGKDR, the number of particles are set to 2×104 and the training sample size for the calculation of projection matrix is 2×10 3 , accepted from a training set of size 4×10 4 simulated samples using rejection ABC.The dimensionality is set to 5. Cross validation is conducted using a test set of size 2 × 10 4 .Results of SMC ABC are shown inTable-3.AMSEs are reported.The simulation time is shown as well.The computational time of constructing

sufficient dimension is set to 5 ;
a smaller value induces substantial worse results.We simulated a set of 30 parameters, a fixed simulated sample of size 10 7 for all the methods and a training sample of size 10 6 , a test sample of size 10 5 for LGKDR and Semi-automatic ABC.The values of log(r) and φ are fixed as in[14], and log(σ e ) are drawn from an uninformative uniform distribution on [log(0.1),0].