{'fraction_non_alphanumeric': 0.03463373194798405, 'fraction_numerical': 0.014310919427563223, 'mean_word_length': 4.6911491260691704, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Simple Approach to Maximum Intractable Likelihood Estimation

Approximate Bayesian Computation (ABC) can be viewed as an analytic approximation of an intractable likelihood coupled with an elementary simulation step. Such a view, combined with a suitable instrumental prior distribution permits maximum-likelihood (or maximum-a-posteriori) inference to be conducted, approximately, using essentially the same techniques. An elementary approach to this problem which simply obtains a nonparametric approximation of the likelihood surface which is then used as a smooth proxy for the likelihood in a subsequent maximisation step is developed here and the convergence of this class of algorithms is characterised theoretically. The use of non-sufficient summary statistics in this context is considered. Applying the proposed method to four problems demonstrates good performance. The proposed approach provides an alternative for approximating the maximum likelihood estimator (MLE) in complex scenarios.

It is well-known that ABC produces a sample approximation of the posterior distribution (Beaumont et al., 2002) in which there exists a deterministic approximation error in addition to Monte Carlo variability. The quality of the approximation to the posterior and theoretical properties of the estimators obtained with ABC have been studied in Wilkinson (2008); Marin et al. (2011); Dean et al. (2011) and Fearnhead and Prangle (2012). The use of ABC posterior samples for conducting model comparison was studied in Didelot et al. (2011) and Robert et al. (2011). Using this sample approximation to characterise the mode of the posterior would in principle allow (approximate) maximum a posteriori (MAP) estimation. Furthermore, using a uniform prior distribution, under the parameterisation of interest, over any set which contains the MLE will lead to a MAP estimate which coincides with the MLE. In lowdimensional problems if we have a sample from the posterior distribution of the parameters, we can estimate its mode by using either nonparametric estimators of the density or another mode-seeking technique such as the mean-shift algorithm (Fukunaga and Hostetler, 1975). Therefore, in contexts where the likelihood function is intractable we can use these results to obtain an approximation of the MLE.
We will denote the estimator obtained with this approximation AMLE.
Although Marjoram et al. (2003) noted that "It [ABC] can also be used in frequentist applications, in particular for maximum-likelihood estimation" this idea does not seem to have been developed. A method based around maximisation of a non-parametric estimate of the log likelihood function was proposed by Diggle and Gratton (1984) in the particular case of simple random samples; their approach involved sampling numerous replicates of the data for each parameter value and estimating the density in the data space. de Valpine (2004) proposes an importance sampling technique, rather closer in spirit to the approach developed here, by which a smoothed kernel estimation of the likelihood function up to a proportionality constant can be obtained in the particular case of state space models provided that techniques for sampling from the joint distribution of unknown parameters and latent states are available -not a requirement of the more general ABC technique developed below. The same idea was applied and analysed in the context of the estimation of location parameters, with particular emphasis on symmetric distributions, by Jaki and West (2008). The particular case of parameter estimation in hidden Markov models was also investigated by Dean et al. (2011), who relied upon the specific structure of that problem. To the best of our knowledge neither MAP estimation nor maximum likelihood estimation in general, implemented directly via the "ABC approximation" combined with maximisation of an estimated density, have been studied in the literature. However, there has been a lot of interest in this type of problem using different approaches (Cox and Kartsonaki, 2012;Fan et al., 2012;Mengersen et al., 2012) since we completed the first version of this work (Rubio and Johansen, 2012).
The estimation of the mode of nonparametric kernel density estimators which may seem, at first, to be a hopeless task has also received a lot of attention (see e.g. Parzen, 1962;Konakov, 1973;Romano, 1988;Abraham et al., 2003;Bickel and Früwirth, 2006). Alternative nonparametric density estimators which could also be considered within the AMLE context have been proposed recently in Cule et al. (2010); Jing et al. (2012).
The remainder of this paper is organised as follows. In Section 2, we present a brief description of ABC methods. In Section 3 we describe how to use these methods to approximate the MLE and present theoretical results to justify such use of ABC methods. In Section 4, we present simulated and real examples to illustrate the use of the proposed MLE approximation. Section 5 concludes with a discussion of both the developed techniques and the likelihood approximation obtained via ABC in general.

Approximate Bayesian Computation
We assume throughout this and the following section that all distributions of interest admit densities with respect to an appropriate version of Lebesgue measure, wherever this is possible, although this assumption can easily be relaxed. Let x = (x 1 , ..., x n ) ∈ R q×n be a simple random sample from a distribution f (·|θ) with support contained in R q , θ ∈ Θ ⊂ R d ; L(θ; x) be the corresponding likelihood function, π(θ) be a prior distribution over the parameter θ and π(θ|x) the corresponding posterior distribution. Consider the following approximation to the posterior where is an approximation of the likelihood function and K ε (x|y) is a normalised Markov kernel. K ε (·|y) is typically concentrated around y with ε acting as a scale parameter. It's clear that (2) is a smoothed version of the true likelihood and it has been argued that the maximisation of such an approximation can in some circumstances lead to better performance than the maximisation of the likelihood itself (Ionides, 2005), providing an additional motivation for the investigation of MLE via this approximation. The approximation can be further motivated by noting that under weak regularity conditions, the distribution π ε (θ|x) is close (in some sense) to the true posterior π(θ|x) when ε is sufficiently small. The simplest approach to ABC samples directly from (1) by the rejection sampling approach presented in Algorithm 1.

Algorithm 1
The basic ABC algorithm.
1: Simulate θ ′ from the prior distribution π(·). 2: Generate y from the model f (·|θ ′ ). 3: Accept θ ′ with probability ∝ K ε (x|y) otherwise return to step 1. Now, let η : R n·q → R m be a summary statistic, ρ : R m × R m → R + be a metric and ε > 0. The simplest ABC algorithm can be formulated in this way using the kernel The ABC rejection algorithm of Pritchard et al. (1999) can be obtained simply by setting η(x) = x.
Several improvements to the ABC method have been proposed in order to increase the acceptance rate, see Beaumont et al. (2002), Marjoram et al. (2003) and Sisson (2007) for good surveys of these. An exhaustive summary of these developments falls outside the scope of the present paper.

Algorithm
Point estimation of θ, by MLE and MAP estimation in particular, has been extensively studied (Lehmann and Casella, 1998). Recall that the MLE, θ, and the MAP estimatorθ are the values of θ which maximise the likelihood or posterior density for the realised data.
These two quantities coincide when the prior distribution is constant (e.g. a uniform prior π(θ) on a suitable (necessarily compact) set D which contains θ). Therefore, if we use a suitable uniform prior, it is possible to approximate the MLE by using ABC methods to generate an approximate sample from the posterior and then approximating the MAP using this sample. In a different context in which the likelihood can be evaluated pointwise, simulation-based MLEs which use a similar construction have been shown to perform well (see, e.g., Gaetan and Yao, 2003, Lele et al., 2007and Johansen et al., 2008. In the present setting the optimisation step can be implemented by estimating the posterior density of θ using a nonparametric estimator (e.g. a kernel density estimator) and then maximising this function: Algorithm 2.
We note that we have not here considered similar simulation-based approaches to the direct optimisation of the likelihood function for a number of reasons. One is computational cost (not having access to the likelihood even pointwise means that distributions concentrated around the mode could be constructed only by introducing several replicates of the data and the rejection or other mechanism used to produce samples from this distribution will become increasingly inefficient as the number of replicates increases); another is that the proposed method has the additional advantages that it fully characterises the likelihood surface and can be conducted concurrently with Bayesian analysis with no additional simulation effort.
Note that the first step of this algorithm can be implemented rather generally by using essentially any algorithm which can be used in the standard ABC context. It is not necessary to obtain an iid sample from the distribution π ε : provided the sample is appropriate for approximating that distribution it can in principle be employed in the AMLE context (although correlation between samples obtained using MCMC techniques and importance weights and dependence arising from the use of SMC can complicate density estimation, it is not as problematic as might be expected (Sköld et al., 2003)).
A still more general algorithm could be implemented: using any prior which has mass in some neighbourhood of the MLE and maximising the product of the estimated likelihood and the reciprocal of this prior (assuming that the likelihood estimate has lighter tails than the prior, not an onerous condition when density estimation is used to obtain that estimate) will also provide an estimate of the likelihood maximiser, an approach which was exploited by de Valpine (2004) (who provided also an analysis of the smoothing bias produced by this technique in their context). In the interests of parsimony we do not pursue this approach here, and throughout the remainder of this document we assume that a uniform prior over some set D which includes the MLE is used, although we note that such an extension eliminates the requirement that a compact set containing a maximiser of the likelihood be identified in advance.
One obvious concern is that the approach could not be expected to work well when the parame-ter space is of high dimension: it is well known that density estimators in high-dimensional settings converge very slowly. Three things mitigate this problem in the present context: • Many of the applications of ABC have been to problems with extremely complex likelihoods which have only a small number of parameters (such as the examples considered below).
• When the parameter space is of high dimension one could employ composite likelihood techniques with low-dimensional components estimated via AMLE. Provided appropriate parameter subsets are selected, the loss of efficiency will not be too severe in many cases. Alternatively, a different mode-seeking algorithm could be employed (Fukunaga and Hostetler, 1975).
• In certain contexts, as discussed below Proposition 2, it may not be necessary to employ the density estimation step at all.
Finally, we note that direct maximisation of the smoothed likelihood approximation (2) can be interpreted as a pseudo-likelihood technique (Besag, 1975), with the Monte Carlo component of the AMLE algorithm providing an approximation to this pseudo-likelihood.

Asymptotic Behaviour
In this section we provide some theoretical results which justify the approach presented in Section 3.1 under similar conditions to those used to motivate the standard ABC approach. We assume throughout that the MLE exists in the model under consideration but that the likelihood is intractable; in the case of non-compact parameter spaces, for example, this may require verification on a case-by-case basis.
We begin by showing pointwise convergence of the posterior (and hence likelihood) approximation under reasonable regularity conditions. It is convenient first to introduce the following concentration condition on the class of ABC kernels which are employed: Condition K A family of symmetric Markov kernels with densities K ε indexed by ε > 0 is said to satisfy the concentration condition provided that its members become increasingly concentrated as ε decreases such that As the user can freely specify K this is not a problematic condition. It serves only to control the degree of smoothing which the ABC approximation of precision ε can effect.
be a compact set that containsθ, the MLE of θ; and let K ε be the densities of a family of symmetric Markov kernels, which satisfies the concentration condition Proof. It follows from the concentration condition that: which converges to 0 as ε → 0 by continuity. Therefore f ε (x|θ) Now, by bounded convergence (noting that boundedness of f ε (x|θ), for ε < ǫ, follows from that of f itself), we have that: The result follows by combining (4) and (5), whenever π(θ|x) is itself well defined.
This result can be strengthened by noting that it is straightforward to obtain bounds on the error introduced at finite ε if we assume Lipschitz continuity of the true likelihood. Unfortunately, such conditions are not typically verifiable in problems of interest. The following result, in which we show that whenever a sufficient statistic is employed the ABC approximation converges pointwise to the posterior distribution, follows as a simple corollary to the previous proposition. However, we provide an explicit proof based on a slightly different argument in order to emphasize the role of sufficiency.
Proof. Using the integral Mean Value Theorem (as used in a similar context by Dean et al. (2011, Equa-tion 6)) we find that for θ ∈ D and any ε ∈ (0, ǫ): for some ξ(θ, x, ε) ∈ B ε , where λ is the Lebesgue measure. Then As this holds for any sufficiently small ε > 0, we have by ρ−continuity of f η (·|θ): Using the Dominated Convergence Theorem we have By the Fisher-Neyman factorization Theorem we have that there exists a function h : The result follows by combining (6), (7) and (8).
With only a slight strengthening of the conditions, Proposition 1 allows us to show convergence of the mode as ε → 0 to that of the true likelihood. It is known that pointwise convergence together with equicontinuity on a compact set implies uniform convergence (Rudin, 1976;Whitney, 1991). Therefore, if in addition to the conditions in Proposition 1 we assume equicontinuity of π ε (·|x) on D, a rather weak additional condition, then the convergence to π(·|x) is uniform and we have the following direct corollary to Proposition 1: Corollary 2. Let π ε (·|x) achieve its global maximum atθ ε for each ε > 0 and suppose that π(·|x) has unique maximiserθ. Under the conditions in Proposition 1; if π ε (·|x) is equicontinuous, then lim ε→0 π ε (θ ε |x) = π(θ|x).
That is, in the case of a sufficiently well-behaved density estimation procedure, using the simple form of the ABC estimator (Algorithm 1) we have that for any level of precision γ, the maximum of the AMLE approximation will, for large enough ABC samples, almost surely be γ−close to the maximum of the posterior distribution of interest, which coincides with the MLE under the given conditions. A simple continuity argument suffices to justify the use ofθ m,ε to approximateθ for large m and small ε.
The convergence shown in the above results depends on the use of a sufficient statistic. In contexts where the likelihood is intractable, this may not be available. In the ABC literature, it has become common to employ summary statistics which are not sufficient in this setting. Although it is possible to characterise the likelihood approximation in this setting, it is difficult to draw useful conclusions from such a characterisation. The construction of appropriate summary statistics remains an active research area (see e.g. Peters et al., 2010 andFearnhead andPrangle, 2012).
We finally provide one result which provides some support for the use of certain non-sufficient statistics when there is a sufficient quantity of data available. In particular we appeal to the large-sample limit in which it can be seen that for a class of summary statistics the AMLE can almost surely be made arbitrarily close to the true parameter value if a sufficiently small value of ε can be used. This is, of course, an idealisation, but provides some guidance on the properties required for summary statistics to be suitable for this purpose and it provides some reassurance that the use of such statistics can in principle lead to good estimation performance. In this result we assume that the AMLE algorithm is applied with the summary statistics filling the role of the data and hence the ABC kernel is defined directly on the space of the summary statistics.
In order to establish this result, we require that, allowing η n (x) = η n (x 1 , . . . , x n ) to denote a sequence of d η -dimensional summary statistics, the following four conditions hold: for π − a.e. θ S.ii g : Θ → R dη is an injective mapping. Letting H = g(D) ⊂ R dη denote the image of the feasible parameter space under g, g −1 : H → Θ is an α-Lipschitz continuous function for some α ∈ R + .

S.iii
The ABC kernels, defined in the space of the summary statistics, satisfy condition K, i.e. K η ε (·|η ′ ) it is concentrated within a ball of radius ε for all ε: S.iv The nonparametric estimator used always provides an estimate of the mode which lies within the convex hull of the sample.
Some interpretation of these conditions seems appropriate. The first tells us simply that the summary statistics converge to some function of the parameters in the large sample limit, a mild requirement which is clearly necessary to allow recovery of the parameters from the statistics. The second condition strengthens this slightly, requiring that the limiting values of the statistics and parameters exist in oneto-one correspondence and that this correspondence is regular in a Lipschitz-sense. The remaining conditions simply characterise the behaviour of the ABC approximation and the AMLE algorithm.

(b) The AMLE approximation of the MLE lies within B αε (θ ⋆ ) almost surely.
Proof. Allowing f ηn ε (η|θ) to denote the ABC approximation of the density of η n given θ, we have: where with the final equality following from S.i (noting that almost sure convergence of η n to g(θ) implies convergence in distribution of η n to a degenerate random variable taking the value g(θ)).
And the ABC approximation to the posterior density of θ, lim n→∞ π ε (·|η n ), may be similarly constrained: using S.ii. And by assumption S.i, S.ii and the continuous mapping theorem we have that giving result (a); result (b) follows immediately from S.iv.
It is noteworthy that this proposition suggests that, at least in the large sample limit, one can use any estimate of the mode which lies within the convex hull of the sampled parameter values. The posterior mean would satisfy this requirement and thus for large enough data sets it is not necessary to employ the nonparametric density estimator at all in order to implement AMLE. This is perhaps an unsurprising result and seems a natural consequence of the usual Bayesian consistency results but it does have implications for implementation of AMLE in settings with large amounts of data for which the summary statistics are with high probability close to their limiting values.

Example 1 (Location-Scale Families and Empirical Quantiles). Consider a simple random sample from a location-scale family, in which we can write the distribution functions in the form:
Allow η 1 n =F −1 (q 1 ) and η 2 n =F −1 (q 2 ) to denote to empirical quantiles. By the Glivenko-Cantelli theorem, these empirical quantiles converge almost-surely to the true quantiles: In the case of the location-scale family, we have that: and we can find explicitly the mapping g −1 : which can be assured if F 0 is strictly increasing and q 1 = q 2 . In this case we even obtain an explicit form for α.

Use of kernel density estimators
In this section we demonstrate that the simple Parzen estimator can be employed within the AMLE context with the support of the results of the previous section.
Definition 1. (Parzen, 1962) Consider the problem of estimating a density with support on R n from m independent random vectors (Z 1 , ..., Z m ). Let K be a kernel, h m be a bandwidth such that h m → 0 when m → ∞, then a kernel density estimator is defined by Under the conditions h m → 0 and mh n m / log(m) → ∞ together with Theorem 1 from Abraham et al. (2003), we have thatθ m a.s.
− − →θ as m → ∞. Therefore, the results presented in the previous section apply to the use of kernel density estimation. This demonstrates that this simple non-parametric estimator is adequate for approximation of the MLE via the AMLE strategy, at least asymptotically. This is, of course, just one of many ways in which the density could be estimated and more sophisticated techniques could very easily be employed and justified in the AMLE context.

Examples
We is conducted using the R command 'kde' together with the bandwidth matrix obtained via the smoothed cross validation approach of Duong and Hazelton (2005) using the command 'Hscv' from the R package 'ks' (Duong, 2011). R source code for these examples is available from the first author upon request.

Binomial Model
Consider a sample of size 30 simulated from a Binomial(10, 0.5) withx = 5.53. Using the prior θ ∼ Unif(0, 1), a tolerance ε = 0.1, a sufficient statistic η(x) =x and the Euclidean metric we simulate an ABC sample of size 10, 000 which, together with Gaussian kernel estimation of the posterior, gives the AMLEθ = 0.552.
There are three quantities affecting the precision in the estimation of θ: D, n and ε. Figure 1
Consider the priors µ ∼ Unif(−0.25, 0.25) and σ ∼ Unif(0.75, 1.25) (crude estimates of location and scale can often be obtained from data, justifying such a choice; using broader prior support here increases computational cost but does not prevent good estimation), a tolerance ε = 0.01, a sufficient statistic η(x) = (x, s), the Euclidean metric, an ABC sample of size 5, 000, and Gaussian kernel estimation of the posterior. Figure 3 illustrates Monte Carlo variability of the AMLE of (µ, σ). Boxplots were obtained using 50 replicates of the algorithm.

Financial application
Logarithmic daily return prices are typically modelled using Lévy processes. For this reason, it is necessary to model the increments (logarithmic returns) using an infinitely divisible distribution. It has been found empirically that these observations have tails heavier than those of the normal distribution, and therefore an attractive option is the use of the 4−parameter (α, β, µ, σ) α−stable family of distributions, which can account for this behaviour. It is well known that maximum likelihood estimation for this family of distributions is difficult. Various numerical approximations of the MLE have been proposed (see e.g. McCulloch, 1986;Nolan, 2001). From a Bayesian perspective, Peters et al. (2010) proposed the use of ABC methods to obtain an approximate posterior sample of the parameters. They propose six summary statistics that can be used for this purpose.
Given the symmetry observed and in the spirit of parsimony, we consider the skewness parameter β to be 0 in order to calculate the AMLE of the parameters (α, µ, σ). Based on the interpretation of these parameters (shape, location and scale) and the data we use the priors

Superposed gamma point processes
The modelling of an unknown number of superposed gamma point processes provides another scenario with intractable likelihoods which is currently attracting some attention (Cox and Kartsonaki, 2012;Mengersen et al., 2012). Intractability of the likelihood in this case is a consequence of the dependency between the observations, which complicates the construction of their joint distribution. Superposed point processes have applications in a variety of areas, for instance Cox and Smith (1954) present an application of this kind of processes in the context of neurophysiology. In this example we consider a simulated sample of size 88 of N = 2 superposed point processes with inter-arrival times identically distributed as a gamma random variable with shape parameter α = 9 and rate parameter β = 1 observed in the interval (0, t 0 ), with t 0 = 420. This choice is inspired by the simulated example presented in Cox and Kartsonaki (2012).
In order to make inference on the parameters (N, α, β) using the AMLE approach, we implement two ABC samplers using the priors N ∼ Unif{1, 2, 3, 4, 5}, α ∼ Unif(5, 15), β ∼ Unif(0.25, 1.5), tolerances ε ∈ {0.5, 0.4, 0.3, 0.2, 0.15} and two sets of summary statistics. The first set of summary statistics, proposed in Cox and Kartsonaki (2012) and subsequently used in Mengersen et al. (2012), consists of the mean rate of occurrence, the coefficient of variation of the intervals between successive points, the sum of the first five autocorrelations of the intervals, the mean of the intervals, and the Poisson indexes of dispersion, variance divided by mean, for intervals of length 1, 5, 10 and 20. Cox and Kartsonaki (2012) mention that summary statistics based on the intervals between successive points are likely to be useful when N is small, therefore we consider a second set of summary statistics by adding a ninth quantity based on the third moment: the sample skewness of the intervals between successive points n j=1 (x j −x) 3 /( n j=1 (x j −x) 2 /n) 3/2 . Note that, unlike Cox and Kartsonaki (2012) and Mengersen et al. (2012), we are taking the discrete nature of the parameter N into account. The AMLE approach is still applicable in this context given that the maximisation of the joint posterior distribution of (N, α, β) can be conducted by conditioning on N . We also considered a continuous prior, uniform over [1, 5] and obtained comparable results (not shown) -although, naturally, by using a discrete prior on N instead of a continuous one, the uncertainty in the estimation of α and β is reduced. Although allowing N to take a continuous range of values leads to an analysis which is arguably more immediately comparable to those presented previously in the literature, we prefer to restrict N to a discrete set as this is consistent with the statistical interpretation of the parameter and the possibility of doing so is a clear advantage of the AMLE methodology. Figure 6 shows the Monte Carlo variability, estimated by using 50 AMLE samples, for each of the two AMLE approaches based on ABC samples of size 5000.
We can notice that the precision in the estimation of (α, β) increases faster, as the tolerance decreases, when using 9 summary statistics. We can observe the same phenomenon from Table 1 in the estimation of N . (Note that the horizontal line shows the parameters used to generate the data not the true value of the MLE). The uncertainty in the estimation of α and β using the AMLE approach with the set of 9 summary statistics seems to be qualitatively comparable with that in Cox and Kartsonaki (2012) for a small tolerance ε. Figure 7 shows scatter plots of the AMLE estimators of β and α for ε = 0.15 and both sets of summary statistics. This scatterplot demonstrates that the mean (α/β) of the gamma distribution is much more tightly constrained by the data than the shape parameter, leading to a nearly-flat ridge in the likelihood surface. The variability in the estimated value of α/β is, in fact, rather small; while the variability in estimation of the shape parameter reflects the lack of information about this quantity in the data and the consequent flatness of the likelihood surface.  To show the variability of the estimator with different data, we also compare the variability of the estimators obtained using 50 different data sets. For each data set we obtain the corresponding AMLE of (N, α, β) by using the priors N ∼ Unif{1, 2, 3, 4, 5}, α ∼ Unif(5, 13) and β ∼ Unif(0.5, 1.5), tolerances ε ∈ {0.5, 0.4, 0.3, 0.2, 0.15} and the two sets of summary statistics mentioned above. Figure   8 shows the boxplots of the AMLEs for (α, β) obtained using ABC samples of size 5000. We can observe that the behaviour of the estimators of (α, β) is fairly similar for both sets of summary statistics. Table 1 also suggests an improvement in the estimation of N produced by the inclusion of the sample skewness.

Discussion
This paper presents a simple algorithm for conducting maximum likelihood estimation via simulation in settings in which the likelihood cannot (readily) be evaluated and provides theoretical and empirical support for that algorithm. This adds another tool to the "approximate computation" toolbox. This allows the (approximate) use of the MLE in most settings in which ABC is possible: desirable both in itself and because it is unsatisfactory for the approach to inference to be dictated by computational considerations. Furthermore, even in settings in which one wishes to adopt a Bayesian approach to inference it may be interesting to obtain also a likelihood-based estimate as agreement or disagreement between the approaches. Naturally, both ABC and AMLE being based upon the same approximation, the difficulties and limitations of ABC are largely inherited by AMLE. Selection of statistics in the case in which sufficient statistics are not available remains a critical question. There has been considerable work on this topic in recent years (see e.g. Fearnhead and Prangle, 2012).
A side-effect of the AMLE algorithm is an approximate characterisation of the likelihood surface, or in Bayesian settings of the posterior surface. We would strongly recommend that this surface be inspected whenever ABC or related techniques are used as even in settings in which the original likelihood contains strong information about the parameters it is possible for a poor choice of summary statistic to lead to the loss of this information. Without explicit consideration of the approximation, perhaps combined with prior sensitivity analysis, this type of issue is difficult to detect.