High-dimensional ABC
D. J. Nott V. M. -H. Ong Y. Fan S. A. Sisson


Discovering Key Nodes in a Temporal Social Network
    Jinshuo LIU1 
Chenghao MOU2
Donghong JI1 ()
    December 30, 2023
==================================================







§ INTRODUCTION





Other chapters in this volume have discussed the curse of dimensionality that is inherent to most standard ABC methods.  For a p-dimensional parameter of interest θ=(θ_1,…,θ_p)^⊤, ABC implementations make use of a summary statistic s=S(y) for data y∈ Y of dimension q, where typically q≥ p. When either θ or s is high dimensional,  standard ABC methods have difficulty in producing  simulated summary data that are acceptably close to the observed summary s_obs=S(y_obs), for observed data y_obs. This means that standard ABC methods have limited applicability in high dimensional problems.

More precisely, write π(θ) for the prior, 
p(y|θ) for the data model, p(y_obs|θ) for the likelihood function and π(θ|y_obs)∝ p(y_obs|θ)π(θ) for the intractable posterior distribution.
Standard ABC methods based on S(y) typically approximate the posterior as
π(θ|y_obs)≈π_ABC,h(θ|s_obs), where  

    π_ABC,h(θ|s_obs)   ∝∫ K_h(s-s_obs)p(s|θ)π(θ) ds,

and where K_h(u) is a kernel weighting function with bandwidth h≥ 0.  
A Monte Carlo 
approximation of (<ref>) involves a kernel density estimation of the 
intractable likelihood based on s-s_obs, the distance between simulated and observed summary statistics. As a result,  the quality of the approximation decreases rapidly as the dimension of the summary statistic q increases, as the distance between s and s_obs necessarily increases with their dimension, even setting aside the approximations involved in the choice of an informative S(y).  

Several authors (e.g. blum10, barber+vw15) have given results which illuminate the way that the dimension of the summary
statistic  q impacts the performance of standard ABC methods.  For example, blum10 obtains the result that the minimal mean squared error
of certain kernel ABC density estimators is of the order of N^-4/(q+5), where N is the number of Monte Carlo samples 
in the kernel approximation.  barber+vw15 consider a simple rejection ABC algorithm where the kernel K_h is uniform,
and obtain a similar result concerned with optimal estimation of posterior expectations.  
biau+cg15 extend the analysis of blum10 using
a nearest neighbour perspective, which accounts for the common ABC practice of choosing h adaptively based on a large pool of samples (e.g. blum+nps13).

Regression adjustments  (e.g. blum10,beaumont+zb02,blum+f10,blum+nps13) are extremely valuable in practice for extending the applicability
of ABC approximations to higher dimensions, since the regression model has some ability to compensate for the mismatch between the simulated
summary statistics s and the observed value s_obs.   Regression adjustment
However, except when the true relationship between θ and s is known precisely (allowing for a perfect adjustment), these approaches may only extend ABC applicability to moderately higher dimensions. For example, nott+fms14 demonstrated a rough doubling of the number of acceptably estimated parameters for a fixed computational cost when using regression adjustment compared to just rejection sampling, for a simple toy model.
Nonparametric regression approaches are also subject
to the curse of dimensionality, and the results of blum10
also apply to certain density estimators which include nonparametric regression adjustments.  Nevertheless, it has been
observed that these theoretical results may be overly pessimistic in practice for some problems.  
See  li+f16 for some recent progress on theoretical aspects of regression adjustment for uncertainty quantification.


This chapter considers the question of whether it may be possible to conduct reliable ABC-based inference for high-dimensional models, or when the number of summary statistics q≥ p is large. 


As a general principle, any methods that improve the efficiency of existing ABC techniques, such as more efficient Monte Carlo sampling algorithms, will as a result help extend ABC methods to higher dimensions, simply because they permit a greater inferential accuracy (measured by an effectively lower kernel bandwidth h) for the same computational overheads. However there is a limit to the extent to which these improvements can produce substantial high-dimensional gains, as ultimately the bottleneck is determined by the s-s_obs term within the kernel K_h embedded as part of the approximate posterior π_ABC,h(θ|s_obs).


Instead, we examine ways in which the reliance on the q-dimensional comparison s-s_obs can be reduced. One technique for achieving this is
by 
estimating low-dimensional marginal posterior distributions for subsets of θ and then reconstructing an estimate
of the joint posterior distribution from these. This approach takes advantage of the fact that the marginal posterior distribution 
π_ABC,h(θ^(1)|s_obs)=∫π_ABC,h(θ|s_obs)dθ^(2) for some partition of the parameter vector θ=(θ^(1)^⊤,θ^(2)^⊤)^⊤
can be much more accurately approximated using ABC directly as π_ABC,h(θ^(1)|s^(1)_obs), since the corresponding necessary set of summary statistics s^(1)⊂ s would be a lower 
dimensional vector compared with the summary statistics s required to estimate the full joint distribution π_ABC,h(θ|s_obs). 
The same idea can also be implemented when
approximating the likelihood function, where it is the sampling distribution of the summary statistics p(s|θ) that is approximated
based on low-dimensional estimates for subsets of s. 

The above techniques are applicable for general ABC inference problems without any particular
exploitable model structure, and are the primary focus of this chapter.  For models with a known exploitable structure it may be possible to achieve better results (e.g. barthelme+c14,white+kp15,bazin+db10,tran+nk16,ong+ntsd16), and we also discuss these briefly. 






§ DIRECT ABC APPROXIMATION OF THE POSTERIOR




In this section we consider direct approximation of the posterior distribution π(θ|s_obs) given the observed summary statistics s_obs.  
We first describe the marginal adjustment approach of nott+fms14, in which the standard ABC approximation of the joint posterior distribution is improved by replacing its univariate margins with more precisely estimated marginal approximations.

These more precise marginal distributions are obtained by implementing standard ABC methods to construct
each univariate marginal posterior separately, for which only low-dimensional summary statistics are required. These univariate marginal posteriors 
then replace the margins in the original approximate joint posterior sample, via an appropriate replacement of order statistics.

While the marginal adjustment can work well
we show an instructive toy example where this strategy fails to adequately estimate 
the posterior dependence structure. We subsequently discuss the Gaussian copula ABC approach of  li+nfs15, which extends the marginal adjustment to improve estimation of all
pairwise dependences of the joint posterior, in combination
with the marginal estimates, by use of a meta-Gaussian distribution fang+fk02. 
These ideas are illustrated by several examples.






 §.§ The marginal adjustment strategy




The marginal adjustment method of nott+fms14 is motivated by the following observation.  Marginal adjustment Suppose we wish to estimate accurately 
the univariate marginal posterior distribution π(θ_j|s_obs) of the parameter θ_j. If we can find a summary statistic, say s^(j)⊂ s, that is nearly marginally sufficient for θ_j in the data model p(y|θ_j), then π(θ_j|s_obs) ≈π(θ_j|s^(j)_obs) and
this summary statistic can be used to obtain marginal ABC posterior inferences about θ_j.    Because θ_j is univariate, the summary statistic s^(j) can be low-dimensional.  

Accordingly, the marginal ABC model takes the form

    [                                         π_ABC, h(θ_j|s^(j)_obs)              ∝∫ K_h(s^(j)-s^(j)_obs)p(s^(j)|θ_j)  π(θ_j) ds^(j);                                                                 =∫∫ K_h(s^(j)-s^(j)_obs)   p(s|θ) π(θ_-j|θ_j)  π(θ_j) dθ_-j  ds ]

where θ_-j denotes the elements of θ excluding θ_j, and π(θ_-j|θ_j) denotes the conditional prior
of θ_-j given θ_j.  


The idea of nott+fms14 is to exploit the observation that marginal posterior inferences are much easier in the ABC framework as they only involve a 
lower dimensional  subset of summary statistics, s^(j)⊂ s.
A sample from the joint ABC posterior π_ABC, h(θ|s_obs) is first obtained, and then this joint sample is adjusted so that it's marginal distributions match
those estimated from the lower-dimensional ABC analyses, π_ABC, h(θ_j|s^(j)_obs).

Write s=(s_1,…,s_q)^⊤ for the summary statistics used to approximate the joint posterior π_ABC, h(θ_j|s_obs), and s^(j)=(s_1^(j),…,s^(j)_q_j)^⊤ for the summary statistics used to approximate the marginal posterior distribution of θ_j, π_ABC, h(θ_j|s^(j)_obs).
The marginal adjustment algorithm is then implemented as follows: 


  * Using standard ABC methods (including regression adjustments) obtain an approximate sample from the joint posterior distribution π(θ|s_obs), θ^J1,…,θ^Jr say, based on the full summary statistic s.  

  * Using standard ABC methods, for each j=1,…,p, obtain an approximate sample from the univariate marginal distribution π(θ_j|s^(j)_obs), θ_j^M1,…,θ_j^Mr' say, based on the lower-dimensional summary statistic s^(j).

  * Write θ_j^M(k) for the k-th order statistic of the (marginally estimated) sample θ_j^M1,…,θ_j^Mr' and 
θ_j^J(k) for the k-th order statistic of the (jointly estimated marginal) sample θ_j^J1,…,θ_j^Jr. Also write 
R(j,k) for the rank of θ_j^Jk within the sample θ_j^J1,…,θ_j^Jr.  Define

    θ^Ak=(θ^M_1(R(1,k)),…,θ^M_p(R(p,k)))^⊤.

Then θ^Ak, k=1,…,r, is a marginally adjusted approximate sample from π(θ|s_obs).  

It is worth stating in words what is achieved by step 3 above.  The samples θ^Ak, k=1,…,r are the same as
θ^Jk, except that componentwise  the order statistics θ_j^J(k) have been replaced by the 
corresponding order statistics θ_j^M(k).  If we were to convert the samples θ^Ak and
θ^Jk to ranks componentwise they would be exactly the same, and so the dependence structure in the original
samples θ^Jk is preserved in θ^Ak in this sense.  However, the estimated marginal distribution
in θ^Ak for θ_j is simply the estimated marginal distribution obtained from the samples θ_j^M1,…,θ_j^Mr', 
so that the adjusted samples θ^Ak give the more precisely estimated marginal distributions from the low-dimensional analyses of step 2, while
preserving the dependence structure from the joint samples of step 1.  

While it is true that the dependence structure obtained at step 1 may not
be well estimated due to standard ABC curse-of-dimensionality arguments, it is also the case that the marginal adjustment improves the estimation of the marginal posterior distributions.

These ideas are illustrated in the following example.





 §.§ A toy example






Following li+nfs15 we let  the data y=(y_1,…,y_p)^⊤, p≥ 2 follow a N(θ,I_p) distribution where θ=(θ_1,…,θ_p)^⊤ is the parameter of interest and I_p denotes the p× p identity matrix. The prior π(θ) is specified as the twisted normal form haario+st99

    π(θ)∝exp(-θ_1^2/200-(θ_2-bθ_1^2+100 b)^2/2-∑_j=3^p θ_j^2)

where we set b=0.1, and if p=2 the ∑_j=3^p θ_j^2 term is omitted.  
A contour plot of π(θ) for p=2 is shown in Figure <ref>.  
This is an interesting example because the likelihood only provides location information
about θ.  The dependence structure in the posterior comes mostly from the prior, and the assocation between θ_1 and
θ_2 changes direction in the left and right tails of the prior (Figure <ref>).  So the posterior dependence changes direction
depending on whether the likelihood locates the posterior in the left or right tail of the prior.  This feature makes it difficult for
standard regression adjustment methods, which merely translate and scale particles (e.g. generated from (s,θ)∼ p(s|θ)π(θ)), to work in high-dimensions.   









Figures <ref> and <ref> show what happens in an analysis of this example with p=5 and p=50 respectively. Four ABC approximation methods are considered with observed data y_obs = (10,0,...,0)^⊤. The contour plots of the bivariate posterior estimates π(θ_1,θ_2|s_obs) are represented by solid lines while the contour plot of the true bivariate margin is represented by grey dashed lines. For both Figures, panel (a) shows the estimates obtained via standard rejection ABC while panels (b), (c) and (d) show estimates obtained after marginal, linear regression and both linear regression and marginal adjustment respectively. Note the regression adjustment step is performed after the rejection sampling stage, and before the marginal adjustment. 

For the case when p=5 (Figure <ref>), rejection sampling alone captures the correlation between θ_1 and θ_2, but the univariate margins are too dispersed.  Performing a marginal adjustment following rejection sampling is not good enough, as it only corrects the margin to the right scale and is not able to recover dependence structure. On the other hand, rejection sampling with linear regression adjustment is able to give a good approximation to the true posterior. Performing a subsequent marginal adjustment (Figure <ref>(d)) shows no further visual improvement. 

The example for p=50 (Figure <ref>) shows both the strengths and limitations of the marginal and regression adjustment strategies. It is very clear that standard rejection ABC estimates do not seem to be using much of the information given by the likelihood, as the posterior estimate follows the shape of the prior distribution π(θ). Performing either regression or marginal adjustment centres the estimates on the right location but the shape of the contour plots for the adjustments are incorrect. Moreover, applying marginal adjustment after regression adjustment corrects the univariate margins well, but does not recover the dependence structure. In this example, all four approaches are not able to recover the dependence structure of the true bivariate posterior. It is worth noting in this example that for the normal case with b=0 the marginal adjustment approach works very well even in high dimensions. 

This example shows some of the limitations of the marginal adjustment strategy. 

One possible approach to improve the estimation of the dependence structure (not discussed by nott+fms14) is to use the marginal adjustment on a reparameterised parameter vector θ^*, where the margins of θ^* account for the dependence structure in θ, while  θ^*_i and θ^*_j, i≠ j remain approximately independent. 

This approach would require some prior knowledge of the dependence structure.

 
Since the key idea of the marginal adjustment approach is to build up a more accurate approximation of the joint posterior from estimates of univariate marginal posterior distributions, it is natural to ask if it is possible to consider estimation of marginal posterior distributions of dimension larger than one and to use these
to help estimate the joint dependence structure of π(θ|s_obs) more accurately.











 §.§ Gaussian copula ABC




One way to implement this idea is the Gaussian copula ABC method of li+nfs15.  Gaussian copula ABC Suppose that
𝒞(u)=P(U_1≤ u_1,…, U_p≤ u_p) is the distribution function of a random vector U=(U_1,…,U_p) where the marginal distribution of each
U_j∼ U(0,1) is uniform.  Then 𝒞(u) is called a copula.  Multivariate distributions can always be written in terms of a copula and their marginal distribution functions, which is an implication of Sklar's theorem <cit.>. This allows for a decoupling of the modelling of marginal distributions and the dependence structure of a multivariate distribution.  One useful type of copula derives from a multivariate Gaussian distribution.  
Suppose that η∼ N(0,C) is a p-dimensional multivariate Gaussian random vector where C is
a correlation matrix.  The distribution of U=(Φ(η_1),…,Φ(η_p))^⊤ where Φ(·) denotes the standard normal distribution function
is then a copula.  This kind of copula, called a Gaussian copula, characterises the dependence structure of a multivariate Gaussian distribution and it is parametrised
by the correlation matrix C.  


Suppose now that we further transform U as
γ=(F_1^-1(U_1),…,F_p^-1(U_p))^⊤ where F_1(·),…,F_p(·) are distribution functions with corresponding density 
functions f_1(·),…,f_p(·). The components of γ then have the marginal densities f_1(·),…,f_p(·) respectively, 
and the dependence structure is being described by the Gaussian copula with correlation matrix C.  Clearly if the densities f_j(·), j=1,…,p are
themselves univariate Gaussian then γ is multivariate Gaussian.  A distribution constructed from a Gaussian copula and given marginal distributions 
is called meta-Gaussian fang+fk02 and its density function is

    h(γ)=|C|^-1/2exp(1/2 z^⊤(I-C^-1)z )∏_j=1^p f_j(γ_j)

where z=(z_1,…,z_p)^⊤ and z_j=Φ^-1(F_j(γ_j)).  

li+nfs15 considered using a meta-Gaussian distribution to approximate the posterior distribution π(θ|s_obs) in ABC.  It is easily seen that a meta-Gaussian distribution is 
determined by its bivariate marginal distributions, so that if we are prepared to accept a meta-Gaussian approximation to the joint posterior distribution
in a Bayesian setting, then it can be constructed based on bivariate posterior marginal estimates.   
Asymptotically the posterior will tend to be Gaussian, but a meta-Gaussian approximation may work well even when we are far from
this situation since it  allows for flexible estimation of the marginal distributions.  As with the marginal adjustment, since the bivariate marginal posterior distributions can 
be estimated using low-dimensional summary statistics, this can help to circumvent the ABC curse of dimensionality in estimation of the joint posterior dependence
structure. 

As before, write s^(j) for the statistics that are informative for ABC estimation of the univariate posterior marginal π(θ_j|s_obs), and now write s^(i,j) for the summary statistics informative
for ABC estimation of the bivariate posterior margin π(θ_i,θ_j|s_obs), i≠ j.  Construction of the Gaussian copula ABC approximation to the posterior π(θ|s_obs) proceeds as follows:


  *  Using standard ABC methods (including regression adjustments), for each j=1,…,p, obtain an approximate sample from the univariate marginal distribution  π(θ_j|s^(j)_obs),  θ_j^U1,…,θ_j^Ur say, based on the lower dimensional summary statistic s^(j).  Use kernel density estimation to construct an approximation ĝ_j(θ_j) to π(θ_j|s^(j)_obs).

  *  Using standard ABC methods, for i=1,…,p-1 and j=i+1,…,p, obtain an approximate sample
from the bivariate marginal distribution π(θ_i,θ_j|s^(i,j)_obs), (θ_i^Bj1,θ_j^Bi1),…,(θ_i^Bjr,θ_j^Bir) say, based on the low-dimensional summary statistics
s^(i,j).


  * Write R(i,j,k) as the rank of θ_i^Bjk within the sample θ_i^Bj1,…,θ_i^Bjr.  With this notation
R(j,i,k), j>i, is the rank of
θ_j^Bik within the sample θ_j^Bi1,…,θ_j^Bir.  Estimate C_ij by Ĉ_ij,  the sample correlation between the vectors

    (Φ^-1(R(i,j,1)/r+1),Φ^-1(R(i,j,2)/r+1),…,Φ^-1(R(i,j,r)/r+1))^⊤

and

    (Φ^-1(R(j,i,1)/r+1),Φ^-1(R(j,i,2)/r+1),…,Φ^-1(R(j,i,r)/r+1))^⊤.


  * Construct the Gaussian copula ABC approximation of π(θ|s_obs) as  the meta-Gaussian distribution with marginal distributions ĝ_j(θ_j), j=1,…,p 
(step <ref>), and Gaussian
copula correlation matrix Ĉ=[Ĉ_ij]_i,j=1,…,p where Ĉ_ij, j>i, is as in step <ref>, Ĉ_ji=Ĉ_ij and Ĉ_ii=1.  

While the estimated correlation matrix Ĉ can fail to be positive definite using this procedure (although this did not occur in our analyses), methods to adjust this can be easily implemented e.g. loland+hhf13. Note that by using
the approximate posterior sample from π(θ_i,θ_j|s^(i,j)_obs) from step 2 and the fitted (bivariate) copula model for the pair, it is possible to
investigate whether the Gaussian copula dependence structure at least represents the true bivariate posterior dependence structure well (though not the full multivariate dependence structure).  This can be supplemented
by  application specific goodness of fit checking of posterior predictive densities based on the joint copula approximation.  





In the twisted normal toy example of Section <ref>, the copula strategy can succeed where the marginal adjustment strategy alone fails. 
Similar to Figure <ref>, Figure <ref> illustrates both the bivariate estimates of π(θ_1,θ_2|s_obs)  based on the Gaussian copula ABC approximation (black solid lines) and the true margins (grey dashed lines) for the p=50 dimensional model. From the contour plots, the ABC copula approximation is able to produce estimates largely similar to the true bivariate margins, in stark contrast to the marginal adjustment alone in Figure <ref>. Thus, in this example where standard ABC sampling with regression and/or marginal adjustment fails, the copula strategy succeeds. 

In order to investigate the performance of each ABC posterior estimation method more precisely, we follow li+nfs15 and vary the dimension of the model,  p, from 2 to 250. Table <ref> shows the mean estimated Kullback-Leibler (KL) divergence between the true bivariate margin π (θ_1,θ_2 | s_obs) and the bivariate margin of the full ABC posterior approximation based on 100 replicate approximations, for all five approaches. 

Observe that when the dimension p increases the performance of the standard rejection ABC approach deteriorates. Adopting any of the adjustment strategies improves the overall performance but the estimated KL divergences still increase with dimension p up to fixed limits). This suggests that if accurate estimation of the posterior dependence structure is important, then regression and marginal adjustment strategies alone may be limited to low dimensional models. From Table <ref> it is clear that Gaussian copula ABC outperforms all other methods in terms of KL divergence and its performance does not deteriorate with increasing dimension, p. This is not surprising as the Gaussian copula ABC approximation  is constructed from bivariate estimates of π (θ_1,θ_2 | s_obs), and is therefore able to capture the dependence structure of all bivariate pairs of the full posterior distribution π(θ|s_obs). 





In the following  sections we implement Gaussian copula ABC for two real data analyses: an analysis of multivariate currency exchange data, and simultaneous estimation of multiple quantile regressions.





 §.§ A multivariate g-and-k model for a foreign currency exchange data set



 g-and-k distribution

The g-and-k distribution (raynerm02) is a flexible model for univariate data. It is typically specified through its quantile function

    Q(p|A,B,g,k)=A+B[1+c1-{-g z(p)}/1+{-g z(p)}](1+z(p)^2)^k z(p),

where A,B>0,g and k>-0.5 are parameters respectively controlling location, scale, skewness and kurtosis of the distribution. The parameter c is conventionally
fixed at 0.8 (resulting in k>-0.5), and z(p) denotes the p-quantile of the standard normal distribution.  Many distributions can be recovered or well approximated for
appropriate values of A,B,g and k (such as the normal when g=k=0).  Despite its attractive properties as a model, inference using the g-and-k distribution is challenging since the density, given by the derivative of the inverse of the quantile function, has no closed form.  However, since simulation from the model is trivial by 
transforming uniform variates on [0,1] through the quantile function, an ABC implementation is one possible inferential approach.  This idea was first explored by 
peters+s06 and allinghamkm09. Here we consider a multivariate extension of the model developed by drovandi+p11.  
This model has a univariate g-and-k distribution for each margin, and the dependence structure is specified through a Gaussian
copula.  Note that this use of a Gaussian copula to describe the dependence structure in the data model (likelihood) is distinct from the use of a Gaussian
copula to approximate the dependence structure of the posterior distribution.  

Suppose that the data are n independent multivariate realisations y=(y^1,…,y^n) where y^i=(y^i_1,…,y^i_q)^⊤.  
We assume that 
marginally each y^i_j i=1,…,n follows a
g-and-k distribution with parameters (A_j,B_j,g_j,k_j), j=1,…,q.  Gaussian copula ABC approximates the joint distribution of  y^i by a meta-Gaussian distribution, 
with  Gaussian copula correlation matrix C.  For a q-dimensional data model, there are 4q marginal parameters, and q(q-1)/2 
distinct parameters in the correlation matrix, 
giving p=q(q +7)/2 parameters in total.  
We consider 
an analysis of log daily returns for q=16 currencies  (resulting in p=184 parameters) versus the Australian dollar for 1757 trading days covering the period 1st January 2007 to 31st 
December 2013 <cit.>. 


We adopt as a prior on C the distribution obtained by sampling V∼(I_q,q),
 and then rescaling V to be a valid correlation matrix with 1's on the diagonal.  The priors on A, B, g and k
for each marginal are independent and uniform over the parameter support, although we adopted uniform distributions with ranges of [-0.1,0.1], [0,0.05], [-1,1] and [-0.2,0.5] for A_j, B_j, g_j and k_j to produce samples (s,θ) proportional to p(s|θ)π(θ) but restricted to a region of high posterior density following an initial pilot analysis (see e.g. fearnhead+p12). 



Following the strategy of li+nfs15, the following summary statistics
were considered informative for each marginal parameter:  writing L_kj, k=1,2,3 for the quantiles and O_kj, k=1,…,7 for the octiles of y_j^1,…,y_j^n,
the marginally informative summary statistics were chosen as L_2j for A_j, (L_3j-L_1j,(E_7j-E_5j+E_3j-E_1j)/(L_3j-L_1j))^⊤ for
B_j, (L_3j+L_1j-2L_2j)/(L_3j-L_1j) for g_j, and (E_7j-E_5j+E_3j-E_1j)/(L_3j-L_1j) for k_j.  These summary statistic
choices were guided by similar summary statistics in drovandi+p11 and preliminary analyses to determine which sets of the distinct summaries 
were marginally informative for individual parameters.  For pairs of parameters the summary statistics for individual parameters were simply combined.  
For the correlation parameters in the Gaussian copula, we follow drovandi+p11 and use the robust normal scores correlation coefficient for the marginal summary statistic.
 

Contour plots of various estimates of the bivariate (B_1,k_1) posterior marginal distribution using ABC rejection sampling are illustrated in Figure <ref>. The top panels
show estimates using the full (p-dimensional) vector of summary statistics with (a)  regression adjustment and (b) marginal adjustment, respectively. The performance of each approach individually is poor as the distributions do not exhibit the more accurately estimated dependence structures observed in the remaining panels.      
These estimates are based on ABC rejection sampling with both marginal and regression adjustments, using (c) the full vector of summary statistics, and (d)
the marginally informative summary statistics for (B_1,k_1).
The similarity between panels (c) and (d) indicates that the marginally informative summary statistics are indeed highly informative for the parameter pair  (B_1,k_1). Finally, panel (e) illustrates the Gaussian copula ABC approximation. 
The similarity between panels (d) and (e) indicates that the copula model provides an excellent approximation of the bivariate posterior marginal distribution. 











 §.§ A non-linear multiple quantile regression analysis






Quantile regression can provide a robust alternative to standard mean regression. Quantile regression Model estimates obtained at multiple quantile levels can also provide a more complete picture of the conditional distribution between predictor and response. For a regression with a single 
covariate x, and response y, the linear model corresponding to the τ-th quantile,  Q_y(τ|x), is given by

    Q_y(τ|x)=α_τ+β_τx

where the coefficients α_τ and β_τ depend on the quantile level, τ∈(0,1).

Standard methods fit quantile regressions independently for each quantile level, which can lead to problems of quantiles crossing and a lack of borrowing of information across the quantile levels <cit.>.



Bayesian approaches to quantile regression require the specification of a likelihood. However, exact and tractable likelihood functions are often not available for these models. 
Quantile regression requires the inversion of
many conditional quantile distributions, which are often not analytically available, although numerical grid search can be used (e.g. tokdar2012; Reich2010). However, in the presence of larger data sets, numerical grid searches can become
computationally prohibitive, see for example Reich2010 who suggests using approximations as an alternative.



We consider a dataset for analysing immunodeficiency in infants. In the search for reference ranges to help diagnose infant immunodeficiency, Isaac1983 measured the serum concentration of immunoglobulin-G (IgG) in 298 preschool children. We are interested in estimating the IgG conditional quantiles at the  levels τ=0.1,0.2,0.3, 0.7, 0.75, 0.8, 0.95. A quadratic model in age (x) is used to fit the data due to the expected smooth change of IgG with age, so that

    Q_y(τ|x)=α_τ+β_τx+η_τx^2.

Figure <ref> illustrates this dataset. The black lines show the separately fitted regression lines for the different quantile levels, based on a frequentist estimator using the quantreg package in R <cit.>. Since these curves are fitted separately, no correlation is assumed between the quantile curves, and for close quantile levels τ the fitted quantile estimates can easily cross each other. In practice, strong correlations can exist between curves close to each other, and the true quantile levels will not cross.




We follow the linearly-interpolated likelihood function approach of fengch2015 as a data model p(s|θ), while extending their quantile function Q_y(r|x) to contain more than one predictor as in (<ref>). For each observed covariate x_obs,i, i=1,…,n, a synthetic data point y_i can be obtained via

    y_i=Q_y(τ_j|x_obs,i) + Q_y(τ_j+1|x_obs,i) - Q_y(τ_j|x_obs,i)  /τ_j+1-τ_j(u_i-τ_j ),

where u_i∼ U(0,1), where j is determined so that
τ_j< u_i < τ_j+1, and where Q_τ(y|x) is the model (<ref>) which depends on parameters α_τ, β_τ and η_τ. If u_i<τ_1, y_i is generated from a normal distribution centred on y̅_obs, with standard deviation 3 times the sample standard deviation of y_obs, and truncated below Q_y(τ_1|x). Similarly, if u_i>τ_m, we simulate from the same distribution except that it is truncated above Q_y(τ_m|x).

The parameters α_τ, β_τ and η_τ are sampled from multivariate Gaussian prior distributions π(θ), with mean vector and covariance matrix based on the estimates obtained using quantreg. 
This prior is constrained to satisfy the quantile monotonicity condition so that the fitted quantile regression lines do not cross.


The full vector of summary statistics is constructed as
 
    [                                                    s=S(y)=( .      α̂_τ_1,…, α̂_τ_m,  β̂_τ_1, …, β̂_τ_m, η̂_τ_1, …, η̂_τ_m,;                                                               pu_τ_1, …, pu_τ_m, pl_τ_1,…, pl_τ_m, q_1(y), …, q_100(y) .)^⊤ ]

where α̂_τ, β̂_τ and η̂_τ are the independent frequentist estimators for α_τ, β_τ  and η_τ at  quantile level τ, pu_τ is the proportion of data points above the τth quantile curve, pl_τ is the proportion of data points below 
the τth quantile curve, and q_1(y), …, q_100(y) are the 100 equally
spaced quantiles of the data y.  The summary statistics for α_τ_1 are α̂_τ_i, pu_τ_i, pl_τ_i and the closest 20 quantiles q_1(y), …, q_100(y) to the level τ_i. Similarly, for β_τ_j, the marginally
informative summary statistics will be β̂_τ_j, pu_τ_j, pl_τ_j and the closest 20 quantiles q_1(y), …, q_100(y) to the level τ_j; and so on. Then for the summaries of the bivariate margin, (α_τ_i, β_τ_j), we concatenate the two sets of summaries.


The following analysis is based on N=1,000,000 samples (s^(ℓ),θ^(ℓ))∼ p(s|θ)π(θ), ℓ=1,…,N.
We specify the smoothing kernel K_h(·) as uniform over the range (-h,h) and determine h as the 0.001 quantile of the Euclidean distances between observed and simulated summary statistics. Our model simultaneously fits the seven quantile levels shown in Figure <ref>, resulting in a p=21 dimensional model with q=135 total summary statistics. Note that with  the application of post-hoc adjustments,
monotonicity of the conditional quantiles may not be preserved. 
If this occurs, the offending samples may simply be discarded, although a preferable solution is the development of adjustments that flexibly respect constraints.



Figure <ref> (left panel) shows the mean predicted conditional quantile estimates for the levels τ=0.1, 0.3, 0.75, 0.95 based on fitting the seven quantile level model. 
Although the true quantile curves are not known here, we might expect the independently fitted frequentist estimates to provide a reasonable guide to the truth in this analysis. When the sample size is reasonably large (here n=298), the frequentist approach can produce estimators with good properties (such as a reduced chance for neighbouring quantiles to overlap as n gets large).  As a result, in the current example, the frequentist estimates should be expected to produce similar results to the Bayesian approaches, particularly in the non-extreme regions where there is more data. However the Bayesian analyses naturally enforce non-crossing of quantiles, and so are preferable for this reason, in spite of the approximate posterior. Results from three different ABC variants are shown in Figure <ref> (left panel).
For most quantile levels there are small differences between the marginal univariate quantile estimates, although quantile non-crossing is enforced in each of the Bayesian estimates. For the lower τ=0.1 quantile where data is more scarce, increasing the quality of the ABC posterior approximation from standard rejection ABC (dashed line) to regression adjusted ABC (dot-dash line) to regression and marginally adjusted ABC (dotted line), produces a marginal quantile that is increasingly close to the frequentist estimate, and which roughly partitions 10% of the data below it. This suggests that there is some ABC approximation error (although this is less obvious in the upper τ=0.95 quantile), but that this is less apparent the better the ABC approximation becomes.



In the case of these marginal quantile estimates, Gaussian copula ABC produces quantile estimates (not shown) that are highly similar to the regression and marginally adjusted estimates (dotted line). However, the real differences here are in the quality of the dependence structure of the ABC posterior.
 Figure <ref> (right panel)  shows the correlation in the estimated posterior bivariate margins of (α_i,α_j), (β_i,β_j) and (η_i,η_j) for i≠ j when using Gaussian copula ABC (x-axis) and standard ABC with regression and marginal adjustment using the full vector of summary statistics (y-axis). Here

it is evident that Gaussian copula ABC is able to capture correlations in the bivariate margins that are missed by regular ABC, even when using the univariate marginal adjustment. 
The quality of the posterior approximation will be vital when considering analyses that critically depend on  full, multiple quantile inference. This lends support to
 the Gaussian copula approach as a viable ABC model approximation able to capture much of the bivariate dependence structure of π(θ|s_obs).











§ ABC APPROXIMATION OF THE SAMPLING DISTRIBUTION OF SUMMARY STATISTICS





An alternative to direct ABC approximation of the posterior distribution π(θ|s_obs) is to instead approximate the sampling distribution of summary statistics p(s|θ)
(leuenberger+w10,fan+ns13), thereby approaching the intractable likelihood problem from the more usual ABC conditional density estimation perspective.
The resulting estimated density is then an analytically tractable approximation of the likelihood function for a Bayesian 
analysis using conventional Bayesian computational tools.  Such approaches may be preferable in problems where inference is required for 
multiple datasets arising from the same model.


One way to achieve this is to first
estimate the joint distribution of (s,θ) flexibly and to then condition on observing s=s_obs in the joint model.  This approach was considered by bonassi+yw11 using multivariate normal mixture 
models for the density estimator on (s,θ).  

Synthetic  likelihood <cit.> is another method that directly approximates  the likelihood via an assumed density such as p(s|θ)≈ N_q(μ(θ),Σ(θ)) where the mean μ(θ) and covariance matrix Σ(θ) are unknown functions of the parameter θ. Various techniques are then needed to estimate θ. For further details on synthetic likelihoods see e.g. wood10, wood18 and drovandi+gkr18.






 §.§ A flexible regression density estimator




We describe the flexible conditional density estimation approach of fan+ns13.  Regression density estimation
As with other ABC density estimators, it is constructed from a sample of N summary statistic and parameter pairs
(s^1,θ^1),…,(s^N,θ^N) drawn from a distribution p(s|θ)h(θ).  Note that while the summary statistics are generated given θ from the sampling distribution for the intractable model of interest, the parameters are not necessarily generated from the prior.  Instead, h(θ) is a distribution
chosen to reflect the region over which the likelihood should be well approximated.  Some rough knowledge of the high likelihood region of the parameter space, 
perhaps based on an initial pilot analysis, is useful for setting h(θ).  The method of
fan+ns13 is based on relating
the summary statistics s 
to θ 
by regression approximations, and so it is useful if the actual relationships between s and θ are as simple as possible. One convenient procedure to achieve this is the
semi-automatic summary statistic
approach of fearnhead+p12 which constructs one summary statistic per parameter where each summary statistic is an estimate of the posterior mean value of the parameter, based on a pilot run.  

That is, s_k is the univariate summary statistic informative for θ_k, k=1,…, p, with 
s^j=(s^j_1,…,s^j_p).  

The first step 
is to build marginal regression models for each component of s conditional on θ.  The training data (s_k^1,θ^1),…,(s_k^N,θ^N) is used to build the marginal model for
s_k resulting in an estimated marginal density f̂_k(s_k|θ) for s_k.  fan+ns13 use a fast variational method for fitting mixture of heteroscedastic regression models (nott+tvk12,tran+nk12)  
for the conditional density estimation.  

Then a conditional density estimate for the joint distribution of s given θ is constructed, using a method closely related to that 
considered in giordani+mtk13 for the unconditional case.  The data (s^j,θ^j) are transformed to (U^j,θ^j), where U_k^j=Φ^-1(F̂_k(s_k^j|θ^j)), 
where F̂_k(s_k|θ) is the distribution function corresponding to the density f̂_k(s_k|θ).  If the marginal densities for each
s_k are well estimated, the transformation to U^j makes each component of U^j approximately standard normal regardless of the value of θ.  
A mixture of normals model is then fitted to the data (U^j,θ^j), j=1,…,N.  Write the fitted normal mixture as 

    ∑_k=1^K w_k N(μ_k,Ψ_k)

where N(μ,Ψ) denotes the multivariate normal distribution with mean μ and covariance matrix Ψ, (μ_k,Ψ_k), k=1,…,K are means
and covariances of K normal mixture components, and w_k, k=1,…,K are mixing weights, w_k≥ 0, ∑_j=1^Kw_j=1.  The mixture model
for the joint distribution of (U,θ) then implies a normal mixture model for the conditional density of U|θ, 

    ∑_k=1^K w_k^c N(μ_k^c,Ψ_k^c)

where

    w_k^c=w_k ϕ(θ;μ_k,Ψ_k)/∑_j=1^K w_j ϕ(θ; μ_j,Ψ_j)

are mixing weights with ϕ(θ;μ,Ψ) denoting the multivariate normal density function in θ with mean μ and covariance matrix Ψ, 
and μ_k^c and Ψ_k^c are the conditional mean and covariance of U given θ in the k-th multivariate normal component
N(μ_k,Ψ_k) in the joint mixture model.  Write ĝ(U|θ) for the resulting estimated conditional density of U given θ.  
Inverting the transformation of s to U then produces an estimate of the conditional density of s given θ, 

    L̂(s|θ)     =ĝ(U|θ)∏_j=1^K f̂_j(s_j|θ)/ϕ(U_j;0,1).

An approximation of the observed data likelihood is then given by L̂(s_obs|θ).

The purpose of the transformation from s to U is to simplify the mixture modelling of the joint distribution (U,θ) compared to what would be required
to estimate the joint distribution of (s,θ).  Note that in L̂(s|θ) the marginal density of s_k is not exactly f̂(s_k|θ) due
to the fact that the estimated marginal distributions in ĝ(U|θ) are not exactly standard normal.  
giordani+mtk13 suggest replacing the ϕ(U^j;0,1)  in (<ref>) by its exact marginal distribution in ĝ(U|θ), but fan+ns13 found that
good approximations to L(s_obs|θ) were obtained without this step.  

The above conditional density estimation method seeks to estimate each univariate marginal conditional distribution s_k|θ arbitrarily well, while approximating the overall joint dependence structure by a mixture of normals model. This approach can work well in relatively high dimensions, in the order of tens to hundreds, provided that the dependence structure is relatively straightforward to capture. This also underlines the importance of techniques that can produce summary statistics with simple relationships to θ, such as the method developed by fearnhead+p12.






 §.§ Analysis of stereological extremes
 Stereological extremes




To illustrate the regression density estimation approach we reanalyse a dataset originally analysed using ABC methods by bortot+cs07, and which was previously considered in sisson+fb18.  The data comprise information about the intensity and size distribution of inclusions
 in a 3 dimensional block of clean steel, with the recorded observations being the inclusion sizes (above a threshold of ν_0=5μm), and their number, observed in a 2-dimensional cross-section.  


bortot+cs07 considered models assuming spherical or ellipsoidal inclusion shapes. For the elliptical model the inclusion size is the length of the major axis of the two-dimensional planar ellipse.
In both models the locations of the inclusions above 5μm in size follow a Poisson process with intensity λ.  
Conditional on having an inclusion larger than ν_0, the distribution of the inclusion size is generalized Pareto, with scale parameter σ>0
and shape parameter ξ.  So in both models there
are 3 parameters, θ=(λ,σ,ξ)^⊤.  For the analysis the priors are logλ∼ N(0,100^2), σ∼(0.01,0.0001) and
ξ∼ N(0,100^2).  For the spherical inclusion model it is possible to directly evaluate the likelihood, but for the ellipsoidal inclusion model this is not possible and so
ABC methods are an attractive option.  Here we focus on the ellipsoidal inclusion model.
An analysis of standard rejection ABC with regression adjustment for the spherical model can be found in erhardt+s16.

The high-dimensionality aspect of this analysis comes from the number of summary statistics, rather than the number of parameters.
The summary statistics used comprise the logarithm of the number of inclusions observed in the two-dimensional cross-section (s_1=111), and s_j+1=log (q_(j+1)-q_(j)), j=1,…,111, 
where the q_(j), j=1,…,112 are 112 equally spaced quantiles of the observed inclusion sizes.  This gives q=112 summary statistics in total, and corresponds to conditional density estimation in 112+3=115 dimensions.


The conditional density estimation method requires the choice of h(θ). This is achieved via a pilot analysis by firstly sampling values (s^i,θ^i), i=1,…,n, 
 where the θ^i are sampled from a uniform distribution over a range wide enough to include the support of the posterior and
the s^i are sampled from p(s|θ^i).  
The sample mean μ̂ and covariance matrix  Σ̂ are then calculated for those θ values for which 
s^i-s_obs≤ 20.
The distribution h(θ) is then specified as the truncated normal distribution

    h(θ)∝ N(μ̂,Σ̂)I((θ-μ̂)^⊤Σ̂^-1(θ-μ̂)<9).

The conditional density estimation method for estimating p(s|θ) is then implemented using N=5000 draws (s^i,θ^i)∼ p(s|θ)h(θ), i=1,…,N.


For comparison with the 115-dimensional regression density estimation approach, an additional analysis is performed in only 6-dimensions, using the 3-dimensional summary statistics obtained using
the semi-automatic method of 
fearnhead+p12.   
Figure <ref> shows pairwise scatterplots of the components of s and θ for the samples generated from h(θ)p(s|θ) (plotting the fearnhead+p12 statistics analysis for clarity).  

The resulting scatterplots
after fitting the flexible models f_k(s|θ) to 
the univariate marginal distributions and transforming to the statistics U are illustrated in Figure <ref>.  
Clearly the dependence structure
has been greatly simplified, which facilitates the accurate mixture modelling of (U,θ).  



The histograms in Figure <ref> show the regression density estimated
marginal posterior distributions obtained by using
the original 112 summary statistics (top panels), and the lower dimensional fearnhead+p12 statistics (bottom panels). The solid line illustrates the density estimates obtained by the `gold standard' ABC-MCMC analysis
of bortot+cs07  using large computational overheads.
It is apparent that even when modelling the original high-dimensional set of summary statistics, reasonable answers are obtained using the regression density approach, although using the same method but with
the fearnhead+p12 summary statistics naturally results in an improved performance.  















§ OTHER APPROACHES TO HIGH-DIMENSIONAL ABC





Beyond the density estimation techniques described above, there are a few alternative approaches for extending ABC analyses to higher dimensions.
ABC methods have been previously developed for  functional parameters, specifically in the case of non-parametric hierarchical density estimation rodrigues+ns16. However, while these `infinite-dimensional' parameters require the development of specialised ABC methods (such as a functional regression adjustment), the dimensionality of these techniques is strictly not high-dimensional in the sense considered in this chapter.

Various possibilities are available when the model of interest has a known and exploitable structure.

The simplest of these is where the model factorises into a hierarchical structure p(s|θ,ϕ)=f(θ|ϕ)∏_i p_i(s^(i)|θ^(i)) (e.g. bazin+db10), where s^(i) and θ^(i) denote mutually exclusive partitions of s and θ. In this case, the ABC approximation to the joint posterior π(θ|s_obs) may be naturally constructed using the lower dimensional comparisons s^(i)-s^(i)_obs only.



When the data model can be written in a conditional factorisation form
 p(s|θ) = p(s_1|θ)∏_i=1^q p(s_i|s_1:(i-1),θ), where s_1:k=(s_1,…,s_k)^⊤, and where conditional simulation from p(s_i|s_1:(i-1),θ) is possible, barthelme+c14 (see also white+kp15) proposed an expectation-propogation ABC scheme.

If s_i is low dimensional then p(θ|s_i,obs,s_1:(i-1),obs) (that is, the posterior obtained by matching s_i-s_i,obs based on simulating conditionally on s_1:(i-1),obs) can be well estimated via regular ABC, which barthelme+c14 then approximate by a Gaussian density. Expectation propagation This leads to a Gaussian approximation of p(θ|s_obs), which may be accurate if the number of summary statistics is large. This may be realistic if
the summary statistics S(y)=y are the observed data.
See barthelme+cc17 (this volume) for further details of this approach.



kousthanas+lhqfw15 consider constructing an MCMC-sampler with univariate updates to sample from each univariate conditional distribution π(θ_i|θ_-i,s_obs). Here they note that if a low-dimensional summary statistic can be identified that is sufficient for the conditional distribution of θ_i|θ_-i, then an ABC-MCMC sampler can be implemented that compares summary statistics of a much lower dimension than than the full vector s at each update step. They demonstrate this approach on a high-dimensional linear model with univariate summary statistics for each parameter update.






Finally, synthetic likelihood methods were discussed in Section <ref> as a method to approximate the likelihood function using an assumed parametric form e.g. p(s|θ)≈ N_q(μ(θ),Σ(θ)) <cit.>. Synthetic likelihood As this technique relies on estimating μ(θ) and Σ(θ) for each θ based on a potentially large number of Monte Carlo samples from p(s|θ), this approach can have high computational overheads. Variational Bayes has only recently been considered as a possible approach for fitting intractable models with synthetic likelihoods, but with greatly reduced computational costs. This then allows higher dimensional analyses to be implemented. See e.g. tran+nk16 and ong+ntsd16 for further details on this technique.









§ DISCUSSION




Given that a direct ABC approximation of the joint posterior distribution π(θ|s_obs) involves a kernel density approximation of the likelihood, where the dimensionality involved is the dimension of the summary statistic s, it might initally seem that development of useful, general purpose methods for high-dimensional ABC may not be possible.  However, if we are prepared to step away from the limiting comparison of s-s_obs within the likelihood approximation of standard ABC methods, and build an approximations to π(θ|s)  or p(s|θ) from approximations of lower dimensional distributions,
then it may be possible to develop useful ABC posterior approximations even in high dimensional settings.  The key idea in these approaches is that instead of matching a single vector of summary statistics in high dimensions, s-s_obs, we instead match many different low dimensional summary statistic vectors in constructing our joint posterior approximation. 

While the methods described here will not always work for posterior distributions with a highly complex dependence structure, or in very high dimensions,
further development 
of related methods using the same “divide and conquer" strategy may be a promising direction for future research in high-dimensional ABC.    

This may be particularly true for those methods that are easily parallelisable in their implementation.

Another area that perhaps has good potential for future research involves those techniques related to pseudo-marginal MCMC methods (see andrieu+lv17), which is currently seeing a surge of research interest beyond ABC. These methods have opened up ways to perform exact estimation and sampling for models with intractable likelihood functions, the ideas of which can be extended to implement various forms of approximation of posterior distributions.  These include synthetic likelihoods (see drovandi+gkr18) and variational Bayes methods, which can both be fast to implement, and for which the latter tends to underestimate uncertainty.

The extension of likelihood-free inference methods to problems of higher dimension is a very active research area and promises to be so for the forseeable future.




§ ACKNOWLEDGEMENTS


SAS is supported by the Australian Research Council under the Discovery Project scheme (DP160102544), and the Australian Centre of Excellence in Mathematical and Statistical Frontiers (CE140100049).

chicago

empty
