How to choose the scaling matrix in ABC (without cheating!)?

Question

I am doing a numerical experiment involving comparing Approximate Bayesian Computation (ABC) with other methods.

I am simulating data $\boldsymbol{y}$ from a model and I'm using ABC to get a sample from the posterior of the parameters $\boldsymbol{\theta}$. The summary statistics $S(\boldsymbol{y})$ are on very different scales, with covariance strongly dependent on the parameters $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$. To take this into account, in the acceptance step I use the quadratic form: $$ (\boldsymbol{s_i}-\boldsymbol{s}_0)\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}}^{-1}(\boldsymbol{s_i}-\boldsymbol{s}_0) < \epsilon $$ where $\boldsymbol{s_i}$ is a simulated vector of stats, $\boldsymbol{s_0}$ the observed stats and $\epsilon$ the tolerance.

So far I have "cheated" and I have used $\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}_0}$ where $\boldsymbol{\theta}_0$ are the true parameters. I thought that the choice of the scaling matrix was secondary (of little importance as $\epsilon \rightarrow$ 0) but it turns out that in my case if I don't estimate $\hat{\boldsymbol{\Sigma}}_{\boldsymbol{\theta}}$ really close to $\boldsymbol{\theta}_0$ the algorithm is simply stuck (extremely low acceptance rate and few statistics dominating everything).

Does anybody know any practical method that can be used to choose the weighting matrix when the real parameters are unknown? Thanks

Answer 1

One solution to your problem, which you may or may not deem as "cheating", is to use dimension reduction before running your sampler (see this review). The most straight-forward approach that comes to mind (and is also discussed in the aforementioned review) is partial least squares regression (PLS). The first use of PLS in ABC (at least that I know of) was by Wegmann et al.

In this approach, one introduces a "calibration" step before running whatever ABC algorithm (rejection sampling, MCMC, SMC, etc.) wherein a large number of summary statistics are simulated based on draws from the prior and PLS is used to calculate a transform of the summary statistics explaining the most variance in the parameter space. The observed summary statistics are then transformed and in the ensuing ABC rejection steps the distance between transformed simulated and observed summary statistics will be used instead of the original summary statistics. In this way, you will avoid having to calculate (or invert) $\Sigma_\theta$.

The calibration step can be outlined as follows-

1. Sample $B$ draws from $\theta' \sim P(\theta)$
2. Conduct $B$ simulations from $S \sim P(S|\theta')$
3. Estimate PLS transform of of $S$ to $S'$

PLS is available for R and python and C++.