<figure><img src="image_1.png" /></figure>

HAL
open science

# Approximate Bayesian Computation with Deep Learning and Conformal prediction

Meïli Baragatti, Bertrand Cloez, David Métivier, Isabelle Sanchez

<figure><img src="image_2.png" /></figure>

## To cite this version:

Meïli Baragatti, Bertrand Cloez, David Métivier, Isabelle Sanchez. Approximate Bayesian Computation with Deep Learning and Conformal prediction. 2024. hal-04595814

HAL Id: hal-04595814

https://hal.inrae.fr/hal-04595814

Preprint submitted on 3 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Approximate Bayesian Computation with Deep
Learning and Conformal prediction

Meïli Baragatti¹, Céline Casenave², Bertrand Cloez², David Métivier², and
Isabelle Sanchez²

¹MISTEA, Université de Montpellier, INRAE, Institut Agro, Montpellier,
France

²MISTEA, Université de Montpellier, INRAE, Montpellier, France

June 19, 2025

Abstract

Approximate Bayesian Computation (ABC) methods are commonly used to approximate posterior distributions in models with unknown or computationally intractable likelihoods. Classical ABC methods are based on nearest neighbor type algorithms and rely on the choice of so-called summary statistics, distances between datasets and a tolerance threshold. Recently, methods combining ABC with more complex machine learning algorithms have been proposed to mitigate the impact of these “user-choices”. In this paper, we propose the first, to our knowledge, ABC method completely free of summary statistics, distance, and tolerance threshold. Moreover, in contrast with usual generalizations of the ABC method, it associates a confidence interval (having a proper frequentist marginal coverage) with the posterior mean estimation (or other moment-type estimates).

Our method, named ABCD-Conformal, uses a neural network with Monte Carlo
Dropout to provide an estimation of the posterior mean (or other moment type function-
als), and conformal theory to obtain associated confidence sets. Efficient for estimating
multidimensional parameters and amortized, we test this new method on four different
applications and compare it with other ABC methods in the literature. Supplementary
materials are available online (appendix, source codes, notebooks, datasets).

Keywords: Likelihood-free inference · Approximate Bayesian computation · Convolu-
tional neural networks · Dropout · Conformal prediction · probability matching crite-
rion.
# Approximate Bayesian Computation with Deep Learning and Conformal prediction

Meïli Baragatti¹, Bertrand Cloez¹, David Métivier¹, and Isabelle Sanchez¹

¹MISTEA, Université de Montpellier, INRAE, Institut Agro, Montpellier, France

## Abstract

Approximate Bayesian Computation (ABC) methods are commonly used to approximate posterior distributions in models with unknown or computationally intractable likelihoods. Classical ABC methods are based on nearest neighbor type algorithms and rely on the choice of so-called summary statistics, distances between datasets and a tolerance threshold. Recently, methods combining ABC with more complex machine learning algorithms have been proposed to mitigate the impact of these "user-choices". In this paper, we propose the first, to our knowledge, ABC method completely free of summary statistics, distance and tolerance threshold. Moreover, in contrast with usual generalizations of the ABC method, it associates a confidence interval (having a proper frequentist marginal coverage) with the posterior mean estimation (or other moment-type estimates).

Our method, ABCD-Conformal, uses a neural network with Monte Carlo Dropout to provide an estimation of the posterior mean (or others moment type functional), and conformal theory to obtain associated confidence sets. Efficient for estimating multidimensional parameters, we test this new method on three different applications and compare it with other ABC methods in the literature.

**Keywords:** Likelihood-free inference · Approximate Bayesian computation · Convolutional neural networks · Dropout · Conformal prediction · probability matching criterion.

## 1 Introduction

In many situations, the likelihood function plays a key role in both frequentist or Bayesian inference. However, it is not unusual that such likelihood is difficult to work with either because no closed form expression exists, e.g., the corresponding statistical model is too complex, or because its computational cost is prohibitive. For example, the first situation occurs with max-stable processes [de Haan, 1984] or Gibbs random fields [Besag, 1974; Grelaud et al., 2009], while the second case frequently occurs with mixture models [Simola et al., 2021].

In this paper, we will restrict our attention to parameter inference in the Bayesian framework and in particular to Approximate Bayesian Computation (ABC), which is one of the most popular likelihood-free parameter inference method. Indeed, ABC can be an efficient alternative to overcome the situation of intractable likelihood functions [Marin et al., 2012]. Within the ABC framework, the need for computing the likelihood function is substituted by the need to be able to simulate from the statistical model under consideration. The first proposed ABC methods enabled to obtain an approximation of the posterior distribution of interest, from simulations [Pritchard et al., 1999]. For these early methods, the degree of approximation depends strongly on the choices of the tolerance threshold, the distance and of the so-called summary statistics used during the ABC estimation. Unfortunately, identifying the most relevant summary statistics is problem dependent and in many situations the statistics are chosen by the practitioner, even if several studies provide guidance in this choice, see for instance Joyce and Marjoram [2008], McKinley et al. [2009], Blum et al. [2012], Fearnhead and Prangle [2012] among others.
# 1 Introduction

In many situations, the likelihood function plays a key role in both frequentist or Bayesian inference. However, it is not unusual that such likelihood is difficult to work with either because no closed form expression exists, e.g., the corresponding statistical model is too complex, or because its computational cost is prohibitive. For example, the first situation occurs with max-stable processes [de Haan, 1984] or Gibbs random fields [Besag, 1974; Greloaud et al., 2009], while the second case frequently occurs with mixture models [Simola et al., 2021].

To overcome the intractability of the likelihood function, several approaches have been explored, among which the use of composite likelihoods [Lindsay, 1988; Varin et al., 2011], that is linear combinations of smaller order log-likelihoods. The use of composite likelihoods in place of the (full) likelihood has been shown to be very efficient, but their application is typically limited to the frequentist framework as Bayes' theorem does not apply with composite likelihoods (although some work exists to obtain what we shall call "relevant composite posterior distributions" [Ribatet et al., 2012; Stoehr and Friel, 2015]). Other likelihood-free approaches, called Simulation Based Inference (SBI), use neural conditional density estimation to learn a parametric approximation of the likelihood [Papamakarios et al., 2019] or of the exact posterior [Papamakarios and Murray, 2016]. Both of these can also be approximated using nonparametric approaches; see Gloeckler et al. [2024] who propose a new amortized Bayesian inference method using neural networks.

In this paper, we will restrict our attention to parameter inference in the Bayesian framework and in particular to Approximate Bayesian Computation (ABC), which is one of the most popular likelihood-free parameter inference method. Indeed, ABC can be an efficient alternative to overcome the situation of intractable likelihood functions [Marin et al., 2012]. Within the ABC framework, the need for computing the likelihood function is substituted by the need to be able to simulate from the model under consideration. The first proposed ABC methods enabled to obtain an approximation of the posterior distribution of interest, from simulations [Pritchard et al., 1999]. The simplest sampler is an accept-reject algorithm, which will be referred to as "standard ABC" in the following. This vanilla sampler can be quite inefficient, as numerous draws may be rejected before one is accepted. Alternative samplers have therefore been proposed in a search for computational efficiency, using MCMC or filtering algorithms (see Marjoram et al. [2003], Sisson et al. [2007], Baragatti and Pudlo [2014] and references therein). In particular, the Sequential Monte Carlo ABC (ABC-SMC, Sisson et al. [2007]; Beaumont et al. [2009]; Del Moral et al. [2012]) has become the gold standard among ABC methods. Even though more efficient than the standard ABC, these early ABC samplers are generally difficult to use when simulating complex systems or inferring high-dimensional parameters. This is why, with the increasing complexity of the data and models used nowadays, ABC samplers remain an active field of research. For instance, in order to improve the ABC-SMC sampler, guided approaches have been proposed with the goal of rapidly targeting the region of the parameter space that is compatible with the observed data [Picchini and Tamborrino, 2024]. However, for most of these ABC methods, the degree of approximation of the posterior distribution depends strongly on the choices of
Recently, several approaches managed to mitigate the need of identifying relevant summary statistics. These approaches combine the ABC framework with machine learning methods, such as random forests or neural networks. For instance, Raynal et al. [2018] propose an automatic selection of relevant summary statistics. They propose to consider a large set of potentially relevant summary statistics and use the random forest framework as a black box to get some point estimates of the posterior distribution given these summary statistics. Another approach of interest is the one of Akesson et al. [2022], which uses the estimated posterior mean as a summary statistic, this posterior mean being estimated through a Convolutional Neural Network (CNN). Consequently, since the emergence of ABC methods, a wide variety of methods have been proposed, from the standard "basic" ABC to methods using random forests or neural networks, and including those using MCMC or filtering algorithms (see Marjoram et al. [2003], Sisson et al. [2007], Baragatti and Pudlo [2014] and references therein).

To quantify the uncertainties, Bayesian methods provide posterior credible intervals from the posterior distribution, through posterior quantiles or highest posterior density regions. This can be justified because, in many simple parametric models, the Bernstein-von Mises theorem ensures that posterior credible intervals are asymptotically as well, reconciling the frequentist and Bayesian approaches. However, these credible intervals are not necessarily relevant in practice. Indeed, note for instance that in nonparametric and high-dimensional models, research is still ongoing to understand better the link between the coverage of credible intervals and the coverage of frequentist, see Rousseau and Szabo [2016]; Datta et al. [2000]; Hoff [2023] among others. Indeed, the posterior coverage probability of a Bayesian credible interval does not always match with the corresponding frequentist coverage probability. Quite often, they can even diverge in nontrivial ways, see for instance Wasserman [2011]. In this case, we say that the probability matching criterion is not verified. This criterion can be hard to check when the exact posterior distribution is known, and even more so when the posterior is approximated [Frazier et al., 2018], which is the cases of most applications requiring an ABC method. Another drawback of using posterior credible intervals in ABC methods is that they only take into account the inherent randomness of the model when estimating the prediction uncertainty. Nevertheless, by definition of the method, the a posteriori distribution is approximated, so there is an error due to the method to take into account.

**Our contribution.** In this paper, we propose a new ABC implementation that does not require summary statistics, tolerance threshold and distance. Our method produce point estimates of the posterior distribution e.g. mean, variance, quantiles, without aiming to approximate the whole posterior distribution. It also provide associated confidence sets of the estimates. The estimates are obtained using neural networks with Monte Carlo Dropout while the confidence sets are obtained through conformal prediction. Our ABC estimation method only requires a 1) sampling function, 2) a neural network (NN) whose architecture is tailored for the problem of interest e.g. CNN for images, GNN for data with a graph structure, RNN for sequential data, and 3) a confidence level. While 1) and 3) are common to all ABC methods, requirement 2) is more original, but still easy to implement because well known architectures can be used with minimal changes, i.e. the inclusion of Dropout (or concrete Dropout) layers. Another originality of this work is that it produces directly confidence sets with proper marginal coverage.

We test the proposed ABCD-conformal algorithm in three different problems. We compare the results against the standard ABC, the ABC Random Forest (ABC-RF) of Raynal et al. [2018] and the ABC-CNN of Akesson et al. [2022]. We measure both the accuracy of the estimation and the frequentist coverage of confidence sets.

The paper is organized as follows. In section 2 we briefly present some existing ABC methods: standard ABC, ABC-CNN and ABC-RF. In section 3, we propose the ABCD-conformal methodology, combining neural networks and conformal prediction in a likelihood-free framework. Then, we illustrate and compare all these methods on three examples in section 4, before a discussion in section 5.
the tolerance threshold, the distance and of the so-called summary statistics used during the ABC estimation. Unfortunately, identifying the most relevant summary statistics is problem dependent and in many situations the statistics are chosen by the practitioner, even if several studies provide guidance in this choice (see for instance Joyce and Marjoram [2008], McKinley et al. [2009], Blum et al. [2012], Fearnhead and Prangle [2012] among others).

Recently, several approaches managed to avoid using the tolerance level or to mitigate the need of identifying relevant summary statistics. Most of these approaches combine the ABC framework with machine learning methods, such as random forests or neural networks. For instance, Raynal et al. [2018] propose an automatic selection of relevant summary statistics. They propose to consider a large set of potentially relevant summary statistics and use the random forest framework as a black box to get some point estimates of the posterior distribution given these summary statistics. They do not need any tolerance level, like the previously cited Simulation Based Inference (SBI) methods [Papamakarios and Murray, 2016; Papamakarios et al., 2019], which are called “ε-free” by their authors and use neural networks. Another approach of interest is the one of Akesson et al. [2022], which uses the estimated posterior mean as a summary statistic, this posterior mean being estimated through a Convolutional Neural Network (CNN). To summarise, since the emergence of ABC methods, a wide variety of methods have been proposed, from the standard “basic” ABC to methods using random forests or neural networks.

To quantify the uncertainties, Bayesian methods provide posterior credible intervals from the posterior distribution, through posterior quantiles or highest posterior density regions. This can be justified because, in many simple parametric models, the Bernstein-von Mises theorem ensures that posterior credible intervals are asymptotically as well, reconciling the frequentist and Bayesian approaches. However, these credible intervals are not necessarily relevant in practice. Indeed, note for instance that in nonparametric and high-dimensional models, research is still ongoing to understand better the link between the coverage of credible intervals and the coverage of frequentist confidence intervals, see Rousseau and Szabo [2016]; Datta et al. [2000]; Hoff [2023] among others. Indeed, the posterior coverage probability of a Bayesian credible interval does not always match with the corresponding frequentist coverage probability. Quite often, they can even diverge in nontrivial ways, see for instance Wasserman [2011]. In this case, we say that the *probability matching criterion* is not verified. This criterion can be hard to check when the exact posterior distribution is known, and even more so when the posterior is approximated [Frazier et al., 2018], which is the case of most applications requiring an ABC method. Another drawback of using posterior credible intervals in ABC methods is that they only take into account the inherent randomness of the model when estimating the prediction uncertainty. Nevertheless, by definition of the method, the a posteriori distribution is approximated, so there is an error due to the method to take into account.

**Our contribution.** In this paper, we propose a new ABC implementation that does not require summary statistics, tolerance threshold and distance. Our method produces point estimates of the posterior distribution e.g., mean, variance, quantiles, without aiming to
# 2 A few approaches of Approximate Bayesian Computation

## 2.1 Standard ABC: assets and drawbacks

Consider a parametric statistical model {$f(x | \theta) : x \in \mathbb{R}^d, \theta \in \Theta$} on which we assume a prior distribution $\pi$ on the parameter $\theta$ leading to the Bayesian parametric model {$f(\cdot | \theta), \pi$}. Based on the sample $\mathbf{x} = (x_1, \dots, x_n)$, the cornerstone of any Bayesian inference is the posterior distribution

$$ \pi(\theta | \mathbf{x}) = \frac{f(\mathbf{x} | \theta)\pi(\theta)}{\int_{\Theta} f(\mathbf{x} | \theta)\pi(\theta)d\theta} \propto f(\mathbf{x} | \pi)\pi(\theta), \quad (1) $$

where $f(\mathbf{x} | \theta)$ denotes the likelihood function for the sample $\mathbf{x}$ and parameter $\theta$. In a setting where we assume that the likelihood is not tractable, the simple, standard ABC methodology was the first powerful algorithm to bypass this hurdle (see Pritchard et al. [1999] or Marin et al. [2012]) and, in one of its simplest form, is described in Algorithm 1.

**Algorithm 1:** Pseudo code for the standard ABC sampler

**Input** : A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, the size of the training set $N_{train}$, summary statistics $\eta: \mathbb{R}^d \to \mathbb{R}^l$, a (pseudo)-distance $d$ on $\mathbb{R}^l$, and a tolerance threshold $\alpha \in (0, 1]$.

**Output:** A set $\theta_1, \dots, \theta_{[N_{train}\alpha]}$ whose distribution is approximately $\pi(\theta | \mathbf{x})$

1 Generation of a reference table :

2 for $j \leftarrow 1$ to $N_{train}$ do

3 Draw $\theta_j \sim \pi$;

4 Draw synthetic sample $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^\top$ from the model $f(\cdot | \theta_j)$

5 end

6 Summary statistics and distances :

7 for $j \leftarrow 1$ to $N_{train}$ do

8 Compute the vector of summary statistics $\eta(\mathbf{x}_j)$;

9 Compute the (pseudo)-distance $d_j = d(\eta(\mathbf{x}), \eta(\mathbf{x}_j))$

10 end

11 Order these distances, i.e., $d_{(1)} < d_{(2)} < \dots < d_{(N_{train})}$;

12 Keep the $\theta_j$ corresponding to the $\lfloor N_{train}\alpha \rfloor$ smallest distances.

In Algorithm 1, $d$ is typically a pseudo distance of the form

$$ d: \mathbb{R}^l \times \mathbb{R}^l \to \mathbb{R}_+ \\ (\eta(\mathbf{x}), \eta(\mathbf{x})) \mapsto \|\eta(\mathbf{x}) - \eta(\mathbf{x})\|, $$

where $\eta$ summarizes the data and is called the vector of summary statistics. $\| \cdot \|$ is typically the usual Euclidean norm.

This algorithm samples $\theta$ and $\mathbf{x}$ from the joint posterior density

$$ \pi_{\alpha,d}(\theta, \mathbf{x} | \eta(\mathbf{x})) = \frac{\pi(\theta) f(\mathbf{x} | \theta) \mathbf{1}_{\alpha,d}(\mathbf{x})}{\int \pi(\theta) f(\mathbf{x} | \theta) \mathbf{1}_{\alpha,d}(\mathbf{x}) d\mathbf{x} d\theta}, $$

where $\mathbf{1}_{\alpha,d}(\mathbf{x}) = \mathbf{1}[d(\eta(\mathbf{x}), \eta(\mathbf{x})) < d_{([\alpha N]+1)}]$. The approximate Bayesian computation posterior density is defined as:

$$ \pi_{\alpha,d}(\theta | \eta(\mathbf{x})) = \int \pi_{\alpha,d}(\theta, \mathbf{x} | \eta(\mathbf{x})) d\mathbf{x} $$

The basic idea behind ABC is that using a representative (enough) summary statistic $\eta$ coupled with a small (enough) tolerance $\alpha$ should produce a good (enough) approximation to the posterior distribution;
approximate the whole posterior distribution. It also provides associated confidence sets of the estimates. The estimates are obtained using neural networks with Monte Carlo Dropout, while the confidence sets are obtained through conformal prediction. Our ABC estimation method only requires 1) a sampling function, 2) a neural network (NN) whose architecture is tailored for the problem of interest e.g., CNN for images, GNN for data with a graph structure, RNN for sequential data. Furthermore, to obtain confidence intervals we need 3) a confidence level $\delta$. While 1) is common to all ABC methods, requirement 2) is more original, but still easy to implement because well-known architectures can be used with minimal changes, i.e., the inclusion of Dropout (or concrete Dropout) layers. Note that the choice of the architecture in 2) is not directly comparable to the choice of a summary statistics because the same architecture (with only minor changes) can be reused for very different applications. Another originality of this work is that it directly produces confidence sets with proper marginal coverage, whereas usual ABC methods enable one to obtain credible sets using $\delta$-quantiles of the estimated posterior distribution.

We test the proposed ABCD-Conformal algorithm in four different problems. We compare the results against the standard ABC, the ABC-SMC, the ABC Random Forest (ABC-RF) of Raynal et al. [2018] and the ABC-CNN of Akesson et al. [2022]. We measure both the accuracy of the estimation and the frequentist coverage of confidence sets. We provide a link to the Git repository https://forge.inrae.fr/mistea/codes_articles/abcdconformal where are stored all the codes, instructions and notebooks to reproduce our results in Julia and R. The Julia package ABCMethods.jl [Métivier, 2025] was developed for this paper.

The paper is organized as follows. In section 2 we briefly present some existing ABC methods: standard ABC, ABC-SMC, ABC-CNN and ABC-RF. In section 3, we propose the ABCD-Conformal methodology, combining neural networks and conformal prediction in a likelihood-free framework. Then, we illustrate and compare all these methods on three examples in section 4, before a discussion in section 5.

# 2 A few approaches of Approximate Bayesian Computation

## 2.1 Standard ABC, vanilla sampler

Consider a parametric statistical model {$f(x | \theta) : x \in \mathbb{R}^d, \theta \in \Theta$} on which we assume a prior distribution $\pi$ on the parameter $\theta$ leading to the Bayesian parametric model {$f(\cdot | \theta), \pi$}. Based on the sample $\mathbf{x} = (x_1, \dots, x_n)$, the cornerstone of any Bayesian inference is the posterior distribution,

$$ \pi(\theta | \mathbf{x}) = \frac{f(\mathbf{x} | \theta)\pi(\theta)}{\int_{\Theta} f(\mathbf{x} | \theta)\pi(\theta)d\theta} \propto f(\mathbf{x} | \pi)\pi(\theta), \quad (1) $$

where $f(\mathbf{x} | \theta)$ denotes the likelihood function for the sample $\mathbf{x}$ and parameter $\theta$. In a setting where we assume that the likelihood is not tractable, the simple, standard ABC methodology
that is $\pi_{\alpha,d}(\theta | \eta(\mathbf{x})) \approx \pi(\theta | \mathbf{x})$. This approximation has been validated by different theoretical results. Let us cite the work of Frazier et al. [2018] which prove that the posterior distribution concentrates on sets containing the true parameter under general assumptions and with some rates of convergence. Cite also Fearnhead and Prangle [2012] which describe the best summary statistics to perform this approximation.

Nevertheless, for the algorithm to give a valid approximation of the true posterior distribution, not only $\pi_{\alpha,d}$ has to be close to $\pi$, but also the Monte Carlo approximation error :

$$ \frac{1}{[N\alpha]} \sum_{j,d_j \le d_{(N\alpha)}} \delta_{\tilde{\theta}_j} \approx \pi_{\alpha,d} \quad (2) $$

needs to be considered. In Eq. (2), $d_j$ and $d_{(\cdot)}$ are defined in Algorithm 1. Here it only corresponds to the rate of convergence of the empirical measure in the classical law of large number. This question has recently attracted a great deal of attention Fournier and Guillin [2015]; Dereich et al. [2013]; Weed and Bach [2019]. There is a unavoidable curse of dimensionality : the empirical measure over a $n$-sample (here $n = [N\alpha]$) is at distance $n^{-1/d}$ to the target measure. This explain why the standard ABC lacks efficiency when the dimension of $\mathbf{x}$ is large or when the dimension of the summaries is large. This also induces high computation times. Using a small number of relevant summaries is then crucial, but, when the size of the data tends to infinity, these summaries should converge toward an injective mapping from parameters to summaries. This injectivity depends on both the true structural model and the particular choice of $\eta(\cdot)$, hence there is no general method for finding such statistics $\eta(\cdot)$. This can restrict the (computational) use or the validity of the standard ABC algorithm, because without this condition posterior concentration is not guaranteed.

## 2.2 Towards a preprocessing summary statistics free ABC

In this section, we briefly review two attempts in the ABC literature to obtain methods depending less on the summary statistics choice (and distance) than the original ABC method.

### 2.2.1 ABC-Convolutional Neural Network (Akesson et al. [2022])

The earliest work to bypass the choice of summary statistics in the ABC method is maybe the one of Fearnhead and Prangle [2012]. Indeed, they showed that the optimal summary statistic is the true posterior mean of the parameter of interest. Although it cannot be calculated analytically, they proposed an extra stage of simulation to estimate it. This estimates is then used as a summary statistic within a standard ABC algorithm. In practice, they used linear regression to estimate the posterior mean, with appropriate functions of the data as predictors. This approach has inspired Jiang et al. [2017] who proposed to estimate the posterior means using deep Neural Networks (NN) instead of linear regression. Akesson et al. [2022] goes a bit further by proposing to estimate them using Convolutional Neural Networks (CNN). Indeed, CNNs are often more effective than classical dense NNs for tasks involving structured grid data, such as images, audio, video and time-series data (see for example LeCun et al. [1999] or Boston et al. [2022]). The approach of Akesson et al. [2022] is then a slight modification of the standard ABC: CNN are used to estimate the posterior mean of a parameter, which is then used as a summary statistic. Hence to have asymptotic consistency and valid frequentist coverages, the same kind of conditions are needed than for Algorithm 1. For sake of completeness, this algorithm is briefly described and commented in appendix, see Algorithm 3.

### 2.2.2 ABC-Random Forest [Raynal et al., 2018]

Another interesting approach is the one of Raynal et al. [2018]. Their goal is slightly different from the goal of standard ABC, as they do not seek to approximate the full posterior distribution. Instead, they focused on functional of the posterior distribution, like posterior mean, posterior variance or posterior quantiles. Indeed, they are easier to interpret and report than the whole posterior distribution, and are the main interest of practitioners. Formally, they predict a scalar transform $T(\theta)$ of a multidimensional parameter $\theta$. Without
was the first powerful algorithm to bypass this hurdle (see Pritchard et al. [1999] or Marin et al. [2012]) and, in one of its simplest form, is described in Algorithm 1. In Algorithm 1,

<caption>**Algorithm 1:** Pseudocode for the standard ABC sampler</caption>

**Input** : A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, the size of the training set $N_{\text{train}}$, summary statistics $\eta: \mathbb{R}^d \to \mathbb{R}^l$, a (pseudo)-distance $d$ on $\mathbb{R}^l$, and a tolerance threshold $\alpha \in (0, 1]$.

**Output:** A set $\theta_1, \dots, \theta_{[N_{\text{train}}\alpha]}$ whose distribution is approximately $\pi(\theta | \mathbf{x})$

1 **Generation of a reference table :**

2 **for** $j \leftarrow 1$ **to** $N_{\text{train}}$ **do**

3 **Draw** $\theta_j \sim \pi$;

4 **Draw synthetic sample** $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^\top$ from the model $f(\cdot | \theta_j)$

5 **end**

6 **Summary statistics and distances :**

7 **for** $j \leftarrow 1$ **to** $N_{\text{train}}$ **do**

8 **Compute the vector of summary statistics** $\eta(\mathbf{x}_j)$;

9 **Compute the (pseudo)-distance** $d_j = d(\eta(\mathbf{x}), \eta(\mathbf{x}_j))$

10 **end**

11 **Order these distances, i.e.,** $d_{(1)} < d_{(2)} < \dots < d_{(N_{\text{train}})}$;

12 **Keep the** $\theta_j$ **corresponding to the** $\lfloor N_{\text{train}}\alpha \rfloor$ **smallest distances.**

$d$ is typically a pseudo distance of the form:

$$d: \mathbb{R}^l \times \mathbb{R}^l \to \mathbb{R}_+ \\ (\eta(\mathbf{x}), \eta(\mathbf{x})) \mapsto \|\eta(\mathbf{x}) - \eta(\mathbf{x})\|,$$

where $\eta$ summarizes the data and is called the vector of summary statistics. $\| \cdot \|$ is typically the usual Euclidean norm. This algorithm samples $\theta$ and $\mathbf{x}$ from the joint posterior density

$$\pi_{\alpha,d}(\theta, \mathbf{x} | \eta(\mathbf{x})) = \frac{\pi(\theta) f(\mathbf{x} | \theta) \mathbf{1}_{\alpha,d}(\mathbf{x})}{\int \pi(\theta) f(\mathbf{x} | \theta) \mathbf{1}_{\alpha,d}(\mathbf{x}) d\mathbf{x} d\theta},$$

where $\mathbf{1}_{\alpha,d}(\mathbf{x}) = \mathbf{1}[d(\eta(\mathbf{x}), \eta(\mathbf{x})) < d_{([\alpha N]+1)}]$. The approximate Bayesian computation posterior density is defined as:

$$\pi_{\alpha,d}(\theta | \eta(\mathbf{x})) = \int \pi_{\alpha,d}(\theta, \mathbf{x} | \eta(\mathbf{x})) d\mathbf{x}$$

The basic idea behind ABC is that using representative (enough) summary statistics $\eta$ coupled with a small (enough) tolerance level $\alpha$ should produce a good (enough) approximation to the posterior distribution; that is $\pi_{\alpha,d}(\theta | \eta(\mathbf{x})) \approx \pi(\theta | \mathbf{x})$. Different theoretical results have validated this approximation. Frazier et al. [2018] prove that the posterior distribution concentrates on sets containing the true parameter under general assumptions and
loss of generality, their interest is in

$$ \psi(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \mathbf{x}] = \int_{\Theta} T(\theta)\pi(\theta | \mathbf{x})d\theta = c(\mathbf{x}) \int_{\Theta} T(\theta)f(\mathbf{x} | \theta)\pi(\theta)d\theta, \quad (3) $$

where $c(\mathbf{x})$ is a normalizing constant. From this functional, we easily recover posterior mean, posterior variance or posterior quantiles using respectively $T(\theta) = \theta$, $T(\theta) = \theta^2$ and $T(\theta) = \theta$, and $T(\theta) = \mathbf{1}_{\theta \le q}$.

To mitigate the curse of dimensionality, the dimension of (3) can be reduced by using summary statistics, i.e., the focus is now on

$$ \psi_{\eta}(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \eta(\mathbf{x})] = \int_{\Theta} T(\theta)\pi\{\theta | \eta(\mathbf{x})\}d\theta = c_{\eta}(\mathbf{x}) \int_{\Theta} T(\theta)f\{\eta(\mathbf{x}) | \theta\}\pi(\theta)d\theta, \quad (4) $$

which is hoped to be a good approximation of $\psi(\mathbf{x})$ but much simpler to estimate. Although many studies focus on how to calculate efficiently an average from a sample see for example Gobet et al. [2022] and references therein, the idea of directly approximating the functional of interest in Eq. (3) (or similarly in Eq. (4)) bypasses the costly step of approximating the empirical measure Eq. (2), especially when the dimension is large.

Raynal et al. [2018] then proposed approximations of $\psi_{\eta}(\mathbf{x})$, of the posterior variance $\mathbb{V}_{\pi}[T(\theta) | \eta(\mathbf{x})]$ and of posterior quantiles, using (regression) random forests [Breiman, 2001]. Their algorithm is briefly described and commented in Algorithm 4. A drawback of this method (as presented in Raynal et al. [2018]) is that it only works on unidimensional parameter inference. When the parameter $\theta$ of interest is multidimensional, Raynal et al. [2018] recommend to construct one random forest (RF) for each unidimensional component of $\theta$, and covariance between components is not taken into account. Additional RFs might be constructed if one is interested by estimating posterior covariance between pairs of components. Piccioni et al. [2022] used for instance this approach with an additional sensitivity analysis. To quantify the uncertainty of their estimation, Raynal et al. [2018] proposes, as usual, to use their approximated estimation of quantiles and variances (in Raynal et al. [2018, Section 2.3.3]) but also to use the so-called out-of-bag simulations inherent from RF algorithms (in Raynal et al. [2018, Section 2.3.4]). This latter gives a first hint to quantify uncertainty without giving nevertheless precise confidence intervals.

# 3 The ABCD-conformal algorithm

In this section we propose a new ABC method based on Deep learning with Dropout and conformal prediction, entitled ABCD-conformal.

## 3.1 Motivation and principle

Like Raynal et al. [2018], instead of approximating the whole posterior $\pi(\theta | \mathbf{x})$, we are interested in estimating functional of the form $\psi(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \mathbf{x}]$ (see equation Eq. (3)). The ABCD-conformal method will both output a prediction and an associated exact confidence interval.

Our method uses Neural Networks (NN), similarly to Akesson et al. [2022]. However, in contrast to their approach, the NN directly outputs the functional of interest and an associated uncertainty, enabling us to bypass the nearest neighbors step detailed in Section 2.1 to overcome the shortcomings of standard ABC.

To use conformal prediction, we need both a prediction and an associated measure of uncertainty that we will name heuristic uncertainty and denote $\hat{\mathbb{V}}(\mathbf{x})$. To this end, we will use neural networks with Dropout layers. Dropout layers are classically used to prevent overfitting for neural networks, see Hinton et al. [2012] and Srivastava et al. [2014a]. However, they can also be used to associate uncertainties e.g. variance with predictions from neural networks, see the Bayesian Neural Network literature [Gal, 2016]. We can then obtain valid confidence sets using these uncertainties through conformal prediction [Angelopoulos and Bates, 2023]. These two steps are described in Sections 3.2 and 3.3.
with some rates of convergence. Fearnhead and Prangle [2012] obtain the best theoretical summary statistics to perform this approximation, under the hypothesis of a quadratic loss.

Nevertheless, for the algorithm to give a valid approximation of the true posterior distribution, not only $\pi_{\alpha,d}$ has to be close to $\pi$, but also the Monte Carlo approximation error:

$$ \frac{1}{[N\alpha]} \sum_{j,d_j \le d_{(N\alpha)}} \delta_{\tilde{\theta}_j} \approx \pi_{\alpha,d} \quad (2) $$

needs to be considered. In Eq. (2), $\delta$ denotes the Dirac measure and $d_j$ and $d_{(\cdot)}$ are defined in Algorithm 1. Here it only corresponds to the rate of convergence of the empirical measure in the classical law of large number. This question has recently attracted a great deal of attention [Fournier and Guillin, 2015; Dereich et al., 2013; Weed and Bach, 2019]. There is an unavoidable curse of dimensionality: the empirical measure over a $n$-sample (here $n = \lfloor N\alpha \rfloor$) is at distance $n^{-1/d}$ to the target measure. This explains why the standard ABC lacks efficiency when the dimension of $\mathbf{x}$ is large or when the dimension of the summaries is large. This also induces high computation times. Using a small number of relevant summaries is then crucial, but, when the size of the data tends to infinity, these summaries should converge toward an injective mapping from parameters to summaries. This injectivity depends on both the true structural model and the particular choice of $\eta(\cdot)$, hence there is no general method for finding such statistics $\eta(\cdot)$. This can restrict the (computational) use or the validity of the standard ABC algorithm, because without this condition, posterior concentration is not guaranteed.

## 2.2 Sequential Monte Carlo ABC (ABC-SMC)

The standard ABC can be quite inefficient, as numerous draws may be rejected before accepting one. One of the most widely used alternative, now considered the “gold standard” among ABC methods, is the Sequential Monte Carlo ABC (ABC-SMC, Sisson et al. [2007]; Beaumont et al. [2009]; Del Moral et al. [2012]). In this approach, several weighted populations of samples are generated sequentially, these populations being associated with a decreasing sequence of tolerance levels $\epsilon_1 > \epsilon_2 > \dots > \epsilon_T$. At each iteration, as sample of population $t$ is obtained by randomly drawing a sample from the weighted population at iteration $t-1$ and perturbing it. As the iterations increase and the tolerance levels decrease, these weighted populations approximate the target posterior distribution more and more closely. As in the standard ABC approach, summary statistics are used to summarize the samples. This algorithm is briefly described and discussed in appendix E.1, see Algorithm 3. But it is important to note that compared to standard ABC, this approach is sequential and no longer based on a reference table. As a result, we cannot simulate all the required samples in advance and reuse them for any observation of interest. Another consequence of this sequential accept-reject scheme is that the total number of samples needed to reach a given tolerance level is not known beforehand. Moreover, this approach requires a sequence of decreasing tolerance levels, with the final level defining the approximation of the posterior, or at least an initial level $\epsilon_1$ and a method for decreasing the level between populations.
## 3.2 Prediction and heuristic uncertainty using Monte Carlo Dropout

For an introduction to Deep Learning, we refer the reader to Courville et al. [2016]; Guéron [2019] and references therein. Neural networks have a strong predictive power however they are often regarded as black boxes making their predictions not explainable and without uncertainty. This is an issue in many areas like medicine or autonomous vehicles where false predictions can have big consequences Filos et al. [2020].

A lot of efforts have been made recently either to explain predictions, see the Explainable AI [Xu et al., 2019; Dantas et al., 2023] literature or to associate an uncertainty with the prediction, see Izmailov et al. [2021] or Kompa et al. [2021] among many others approaches.

In this paper we focus on the Dropout method proposed in Gal [2016] and Gal and Ghahramani [2016] even if more recent uncertainty methods might be more powerful as for example Variational BNN methods Folgoc et al. [2021]; Jospin et al. [2022]. We choose the classical Dropout method as it is generic, simple to implement and to integrate into an already existing NN architecture. As explained in Section 3.3, using conformal prediction will enable us to compensate the shortcomings of Dropout.

At every training step, Dropout layers randomly drop (set to zero) some neurons (elements of the weights matrices) [Hinton et al., 2012; Srivastava et al., 2014b] according to some user-choice parameters (such as the Dropout probability). This mechanism prevents complex co-adaptations on the training data: an input or hidden unit cannot rely on other hidden units being present, this can also be seen as a regularization method during the training, to prevent overfitting. Standard Dropout layers are deactivated during test/prediction phase. In their seminal paper Gal and Ghahramani [2016] proposed to activate Dropout during test, introducing randomness in the NN, hence in the output. Repeating the same prediction task several times produces a distribution and is sometimes called Monte Carlo Dropout (MC Dropout).

We briefly describe the MC Dropout procedure. Let $D = \{(x_j, \theta_j), j = 1, \dots, N\}$ be the training set, containing inputs $x$ and outputs $\theta$, where $x$ corresponds in the ABC framework to the sampled data and $\theta$ (or a transformation $T(\theta)$) to the unknown parameter of interest. The vector of parameters of the network (weights and bias) is denoted $\omega$. During training, the goal is to find parameters of the network $\omega$ that are likely to have generated the outputs $(\tilde{\theta}_j, j = 1, \dots, N)$, given the inputs $(x_j, j = 1, \dots, N)$. For this goal, we use another Bayesian inference model to infer such $\omega$ and we then aim to estimate the distribution $\pi(\omega | D)$. After training, the goal is then to predict the parameter of interest $\theta$ associated to a new $x$.

The posterior distribution of interest is then:

$$ \pi(\theta | \mathbf{x}, D) = \int \pi(\theta | \mathbf{x}, \omega)\pi(\omega | D)d\omega. \quad (5) $$

The density $\pi(\omega | D)$ is quite complex and cannot be evaluated analytically. It is approximated using a variational approach by $q(\omega)$. The approximate posterior distribution of interest is then given by

$$ q(\theta | \mathbf{x}) = \int \pi(\theta | \mathbf{x}, \omega)q(\omega)d\omega. \quad (6) $$

The first two moments of this distribution can be estimated empirically following Monte Carlo integration with $K$ samples. For each Monte Carlo sample, different units of the network are randomly dropped out Hinton et al. [2012], we note $\tilde{\omega}_k$ the estimated weights of the associated network with dropped units.

Hence, these weights $\tilde{\omega}_k$ are different for each Monte Carlo iteration. Denote by $\mathbf{f}^\omega(\mathbf{x})$ the model's stochastic output for input $\mathbf{x}$ and parameters $\omega$. Assuming that $\theta | \mathbf{x}, \omega \sim \mathcal{N}(\mathbf{f}^\omega(\mathbf{x}), \tau^{-1}\mathbf{I})$ and $\tilde{\omega}_k \sim q(\omega)$, $t = 1, \dots, K$; an unbiased (consistent) estimator for $\mathbb{E}_{q(\theta|\mathbf{x})}[\theta | \mathbf{x}]$ is given by:

$$ \tilde{\theta}(\mathbf{x}) = \frac{1}{K} \sum_{k=1}^{K} \mathbf{f}^{\tilde{\omega}_k}(\mathbf{x}), \quad (7) $$

which corresponds to the average of $K$ stochastic forward passes through the network with Dropout. An unbiased (consistent) estimator for the second moment $\mathbb{E}_{q(\theta|\mathbf{x})}[\theta^T\theta]$ is given by:

$$ \hat{\mathbb{E}}[\theta^T\theta | \mathbf{x}] = \tau^{-1}\mathbf{I} + \frac{1}{K} \sum_{k=1}^{K} \mathbf{f}^{\tilde{\omega}_k}(\mathbf{x})^T \mathbf{f}^{\tilde{\omega}_k}(\mathbf{x}) $$
It also requires a transition kernel to perturb samples randomly drawn from the previous population. The construction of such kernels, also referred to as proposal samplers, has been studied in particular by Filippi et al. [2013].

## 2.3 Towards a preprocessing summary statistics free ABC

In this section, we briefly review two attempts in the ABC literature to obtain methods depending less on the summary statistics choice (and distance) than the original ABC method.

### 2.3.1 ABC-Convolutional Neural Network [Akesson et al., 2022]

The earliest work to bypass the choice of summary statistics in the ABC method is maybe the one of Fearnhead and Prangle [2012]. Indeed, they showed that the optimal summary statistic is the true posterior mean of the parameter of interest (for quadratic losses). Although it cannot be calculated analytically, they proposed an extra stage of simulation to estimate it. This estimate is then used as a summary statistic within a standard ABC algorithm. In practice, they used linear regression to estimate the posterior mean, with appropriate functions of the data as predictors. This approach has inspired Jiang et al. [2017] who proposed to estimate the posterior means using deep Neural Networks (NN) instead of linear regression. Akesson et al. [2022] go a bit further by proposing to estimate them using Convolutional Neural Networks (CNN). Indeed, CNNs are often more effective than classical dense NNs for tasks involving structured grid data, such as images, audio, video and time-series data (see for example LeCun et al. [1999] or Boston et al. [2022]). The approach of Akesson et al. [2022] is then a slight modification of the standard ABC: CNN are used to estimate the posterior mean of a parameter, which is then used as a summary statistic. Hence, to have asymptotic consistency and valid frequentist coverages, the same kind of conditions are needed as for Algorithm 1. For sake of completeness, this algorithm is briefly described and commented in appendix E.2, see Algorithm 4.

### 2.3.2 ABC-Random Forest [Raynal et al., 2018]

Another interesting approach is the one of Raynal et al. [2018]. Their goal is slightly different from the goal of standard ABC, as they do not seek to approximate the full posterior distribution. Instead, they focused on functionals of the posterior distribution, like posterior mean, posterior variance or posterior quantiles. Indeed, they are easier to interpret and report than the whole posterior distribution, and are the main interest of practitioners. Formally, they predict a scalar transform $T(\theta)$ of a multidimensional parameter $\theta$. Without loss of generality, their interest is in

$$ \psi(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \mathbf{x}] = \int_{\Theta} T(\theta)\pi(\theta | \mathbf{x})d\theta = c(\mathbf{x}) \int_{\Theta} T(\theta)f(\mathbf{x} | \theta)\pi(\theta)d\theta, \quad (3) $$

where $c(\mathbf{x})$ is a normalizing constant. From this functional, we easily recover posterior mean, posterior variance or posterior quantiles using respectively $T(\theta) = \theta$, $T(\theta) = \theta^2$, $T(\theta) = \theta$, and $T(\theta) = \mathbf{1}_{\theta \le q}$.
Then to obtain a variance associated with the prediction of $\theta$ we can use the following unbiased (consistent) estimator:

$$ \hat{V}[\theta | \mathbf{x}] = \underbrace{\tau^{-1}\mathbf{I}}_{\hat{V}_a[\theta|\mathbf{x}]} + \underbrace{\frac{1}{K} \sum_{k=1}^{K} \mathbf{f}^{\hat{\omega}_k}(\mathbf{x})^T \mathbf{f}^{\hat{\omega}_k}(\mathbf{x}) - \hat{\theta}(\mathbf{x})^T \hat{\theta}(\mathbf{x})}_{\hat{V}_e[\theta|\mathbf{x}]} \quad (8) $$

which corresponds to the inverse model precision plus the sample variance of K stochastic forward passes through the network with Dropout. Hence, using Dropout and its associated randomness, we can obtain an estimate of $E_{q(\theta|\mathbf{x})}[\theta]$, as well as associated uncertainty. This estimate can be considered as an approximation of $E_{\pi}[\theta]$ because $q(\theta | \mathbf{x})$ is an approximation of the posterior of interest $\pi(\theta | \mathbf{x}, D)$. Moreover, Gal [2016] interprets the uncertainty Eq. (8) as the sum of an aleatoric uncertainty $\hat{V}_a[\theta | \mathbf{x}]$ and an epistemic uncertainty $\hat{V}_e[\theta | \mathbf{x}]$. The aleatoric uncertainty $\hat{V}_a[\theta | \mathbf{x}]$ is interpreted as the noise in the data, it is the result of measurement imprecision (it is often modelled as part of the likelihood, and this is often a Gaussian corrupting noise). The epistemic uncertainty $\hat{V}_e[\theta | \mathbf{x}]$ is interpreted as the model uncertainty: uncertainty in the model parameters or in the structure of the model. The total predictive uncertainty $\hat{V}[\theta | \mathbf{x}]$ combines these both types of uncertainties and will be used as our heuristic uncertainty $\hat{V}(\mathbf{x})$.

In practice, for each of the K Monte Carlo iterations, the weigths $\hat{\omega}_k$ are different, and the NN with input $\mathbf{x}$ gives as outputs $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$ and $\tau^{-1}_{\hat{\omega}_k}$. The aleatoric uncertainty $\tau^{-1}$ is estimated by the mean of the $\tau^{-1}_{\hat{\omega}_k}$, and the epistemic uncertainty by the sample variance of the $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$.

The main drawback of Dropout layers is that the Dropout probability is a new model hyperparameter. To circumvent this issue, Gal et al. [2017] propose a Dropout variant, called Concrete Dropout, that allows for automatic tuning of this probability in large models, using gradient methods. It improves performance and tunes automatically the Dropout rate, producing better uncertainties compared to classical Dropout. This is the Dropout version we will be using in the following.

**Remark 1.** Gal [2016] and Gal and Ghahramani [2016] showed that optimising any neural network with Dropout is equivalent to a form of approximate inference in a probabilistic interpretation of the model. In other words, the optimal parameters found through the optimisation of a Dropout neural network are the same as the optimal variational parameters in a probabilistic Bayesian neural network with the same structure (the parameters are found through a variational inference approach). It means that a network trained with Dropout is equivalent to a Bayesian Neural Network and possesses all the properties of such a network (see MacKay [1992] or Neal [2012] concerning Bayesian Neural Networks). For this equivalence to be true, only one condition should be verified in the variational inference approach, which is about the Kullback-Leibler divergence between the prior distribution of the parameters $\omega$ and an approximating distribution for $\omega$. We will not need to verify such assumption because the theoretical guaranty of our confidence sets will be assured through the conformal prediction.

## 3.3 Conformal prediction

For a test input $\mathbf{x}$, the NN with concrete Dropout method enables us to obtain an approximation $\hat{\theta}(\mathbf{x})$ of $E_{\pi}[\theta | \mathbf{x}]$, as well as associated heuristic uncertainties $\hat{V}(\mathbf{x})$. However, it does not give valid confidence sets. In order to obtain such confidence sets, we propose to apply a conformal procedure, as explained in Angelopoulos and Bates [2023]. Conformal prediction is indeed a straightforward way to generate valid confidence sets for wide variety of models from a proxy of the uncertainty. The better this proxy will be, the smaller confidence interval will be.

It requires a small amount of additional calibration data $((\theta_j, \mathbf{x}_j), j = 1, \dots, N_{cal})$, hundreds of calibration data are theoretically enough (see Eq. (12) below). We propose to use as heuristic uncertainties $\hat{V}(\mathbf{x})$ one of the variances obtained from the NN with Dropout in Section 3.2, see Eq. (8).

This heuristic uncertainty can be transformed into a rigorous confidence interval $C(\mathbf{x})$ through the following conformal procedure.
To mitigate the curse of dimensionality, the dimension of (3) can be reduced by using summary statistics, i.e., the focus is now on

$$ \psi_{\eta}(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \eta(\mathbf{x})] = \int_{\Theta} T(\theta) \pi\{\theta | \eta(\mathbf{x})\} d\theta = c_{\eta}(\mathbf{x}) \int_{\Theta} T(\theta) f\{\eta(\mathbf{x}) | \theta\} \pi(\theta) d\theta, \quad (4) $$

which is expected to be a good approximation of $\psi(\mathbf{x})$ but much simpler to estimate. Although many studies focus on how to calculate efficiently an average from a sample (see for example Gobet et al. [2022] and references therein), the idea of directly approximating the functional of interest in Eq. (3) (or similarly in Eq. (4)) bypasses the costly step of approximating the empirical measure Eq. (2), especially when the dimension is large.

Raynal et al. [2018] propose approximations of $\psi_{\eta}(\mathbf{x})$, of the posterior variance $\mathbb{V}_{\pi}[T(\theta) | \eta(\mathbf{x})]$ and of posterior quantiles, using (regression) random forests [Breiman, 2001]. Their algorithm is briefly described and commented in appendix E.3, Algorithm 5. A drawback of this method (as presented in Raynal et al. [2018]) is that it only works on unidimensional parameter inference. When the parameter $\theta$ of interest is multidimensional, Raynal et al. [2018] recommend constructing one random forest (RF) for each unidimensional component of $\theta$, and covariance between components is not taken into account. Additional RFs might be constructed if one is interested in estimating posterior covariance between pairs of components. Piccioni et al. [2022] used, for instance, this approach with an additional sensitivity analysis. To quantify the uncertainty of their estimation, Raynal et al. [2018] propose to use their approximated estimation of quantiles and variances (in Raynal et al. [2018, Section 2.3.3]) but also to use the so-called out-of-bag simulations inherent from RF algorithms (in Raynal et al. [2018, Section 2.3.4]). This latter gives a first hint to quantify uncertainty without giving nevertheless precise confidence intervals.

# 3 The ABCD-Conformal algorithm

In this section, we propose a new ABC method based on Deep learning with Dropout and conformal prediction, entitled ABCD-Conformal.

## 3.1 Motivation and principle

Like Raynal et al. [2018], instead of approximating the whole posterior $\pi(\theta | \mathbf{x})$, we are interested in estimating functional of the form $\psi(\mathbf{x}) = \mathbb{E}_{\pi}[T(\theta) | \mathbf{x}]$ (see Eq. (3)). The ABCD-Conformal method will both output a prediction and an associated exact confidence set. Our method uses Neural Networks (NN), similarly to Akesson et al. [2022]. However, in contrast to their approach, the NN directly outputs the functional of interest and an associated uncertainty, enabling us to bypass the nearest neighbors step detailed in Section 2.1 to overcome the shortcomings of standard ABC.

To use conformal prediction, we need both a prediction and an associated measure of uncertainty that we will name heuristic uncertainty and denote $\hat{\mathbb{V}}(\mathbf{x})$. To this end, we will use neural networks with Dropout layers. Dropout layers are classically used to prevent
Using the heuristic uncertainty $\hat{V}(\mathbf{x})$, we compute a score $s(\theta, \mathbf{x})$ for each data sample in the calibration set, from these scores we can define $\hat{q}$ the conformal quantile and then a confidence set. More precisely, as in Messoudi et al. [2022], the score function will be

$$s(\theta, \mathbf{x}) = \sqrt{(\theta - \hat{\theta}(\mathbf{x}))^T \hat{V}(\mathbf{x})^{-1} (\theta - \hat{\theta}(\mathbf{x}))}. \qquad (9)$$

We can then define $\hat{q}$ as the $\frac{[(N_{cal}+1)(1-\alpha)]}{N_{cal}}$ quantile of $\{s(\theta_j, \mathbf{x}_j), j = 1, \dots, N_{cal}\}$ to have as the confidence set:

$$C(\mathbf{x}) = \{\theta \mid s(\theta, \mathbf{x}) \le \hat{q}\}. \qquad (10)$$

This corresponds to an ellipsoid which center is $\hat{\theta}(\mathbf{x})$ and covariance matrix is $\hat{V}(\mathbf{x})^{-1}/\hat{q}^2$. In one dimension, we have $s(\theta, \mathbf{x}) = |\theta - \hat{\theta}(\mathbf{x})|/\sqrt{\hat{V}(\mathbf{x})}$, and

$$C(\mathbf{x}) = [\hat{\theta}(\mathbf{x}) - \hat{q}\sqrt{\hat{V}(\mathbf{x})}; \hat{\theta}(\mathbf{x}) + \hat{q}\sqrt{\hat{V}(\mathbf{x})}]. \qquad (11)$$

This confidence set satisfies the marginal coverage property, i.e. for a chosen $\delta \in [0, 1]$:

$$1 - \delta \le \mathbb{P}[\theta \in C(\mathbf{x})] \le 1 - \delta + \frac{1}{N_{cal} + 1}. \qquad (12)$$

Indeed, see [Angelopoulos and Bates, 2023, Appendix D] (in particular the right hand side inequality needs that random variables $s(\theta_j, \mathbf{x}_j)$ to be continuous). This probability is marginal over the randomness in the calibration and test points [Angelopoulos and Bates, 2023], see Algorithm 2. One great advantage of this procedure is that it is easy to implement, fast and generic. Moreover, the obtained prediction intervals are non-asymptotic, with no distributional or model assumptions. The only assumption needed is that the calibration data and the test data are i.i.d [Vovk et al., 1999] which is naturally the case in our setting. Finally, this procedure provides a frequentist perspective on a Bayesian approach, reconciling the two points of view.

The performance of this approach in terms of intervals length depends only on the quality of the uncertainty measure used. A "good" uncertainty measure should reflect the magnitude of model error: it should be smaller for easy input, and larger for hard ones. If the uncertainty measure is "bad", then the intervals obtained by the conformal procedure will be quite large and their lengths will not differ a lot between easy and hard inputs. From the previous MC Dropout method, three types of uncertainties were obtained: epistemic, aleatoric and overall (sum of aleatoric and epistemic). In section 4, we show that depending on the examples, one or the other of these uncertainties will produce the best results.

**Remark 2.** Note that a conditional coverage would be preferable, that is:

$$\mathbb{P}[\theta \in C(\mathbf{x}) \mid \mathbf{x}] \ge 1 - \delta. \qquad (13)$$

Indeed, marginal coverage does not guarantee the coverage of $C(\mathbf{x})$ for any specific $\mathbf{x}$, only the average coverage over the whole domain. However, this conditional property is not guaranteed by the conformal procedure in general, but some metrics can be used to check how close we are to this stronger property, see Angelopoulos and Bates [2023].

## 3.4 Implementation of ABCD-conformal

The training step of the proposed method is as follows. a) Generate the training dataset i.e. reference table in ABC vocabulary and the calibration dataset, of sizes $N_{train}$ and $N_{cal}$ respectively. b) The NN with Dropout is trained using the reference table. In principle, one could select the network architecture using an extra validation dataset. c) Monte Carlo Dropout prediction is performed, that is each sample of the calibration dataset is passed through the trained network with Dropout, K times. For each of these samples,
overfitting for neural networks, see Hinton et al. [2012] and Srivastava et al. [2014a]. However, they can also be used to associate uncertainties e.g., variance with predictions from neural networks, see the Bayesian Neural Network literature [Gal, 2016]. We can then obtain valid confidence sets using these uncertainties through conformal prediction [Angelopoulos and Bates, 2023]. These two steps are described in Sections 3.2 and 3.3.

## 3.2 Prediction and heuristic uncertainty using Monte Carlo Dropout

For an introduction to Deep Learning, we refer the reader to Courville et al. [2016]; Guéron [2019] and references therein. Neural networks have a strong predictive power, however they are often regarded as black boxes making their predictions not explainable and without associated uncertainty. This is an issue in many areas like medicine or autonomous vehicles, where false predictions can have big consequences [Filos et al., 2020]. A lot of efforts have been made recently either to explain predictions, see the Explainable AI [Xu et al., 2019; Dantas et al., 2023] literature or to associate an uncertainty with the prediction, see Izmailov et al. [2021] or Kompa et al. [2021] among many others approaches.

In this paper we focus on the Dropout method proposed in Gal [2016] and Gal and Ghahramani [2016] even if more recent uncertainty methods might be more powerful as for example Variational BNN methods [Folgoc et al., 2021; Jospin et al., 2022]. We choose the classical Dropout method as it is generic, simple to implement and to integrate into an already existing NN architecture. As explained in Section 3.3, using conformal prediction will enable us to compensate the shortcomings of Dropout. At every training step, in each Dropout layer, some neurons are randomly dropped (their corresponding weights are set to zero) according to some user-choice parameters (such as the Dropout probability). This mechanism prevents complex co-adaptations on the training data. Co-adaptation means that different neurons have highly correlated behavior. Dropout avoid this, because a neuron cannot rely on other neurons being present, hence neurons are encouraged to detect features independently of each other. This can also be seen as a regularization method during the training, to prevent overfitting. Standard Dropout layers are deactivated during test/prediction phase. Details concerning dropout procedures can be found in [Hinton et al., 2012; Srivastava et al., 2014b]. In their seminal paper Gal and Ghahramani [2016] proposed to activate Dropout during test, introducing randomness in the NN, hence in the output. Repeating the same prediction task several times produces a distribution and is sometimes called Monte Carlo Dropout (MC Dropout).

We briefly describe the MC Dropout procedure. Let $D = \{(\mathbf{x}_j, \theta_j), j = 1, \dots, N\}$ be the training set, where $\mathbf{x}_j$ corresponds in the ABC framework to the sampled data and $\theta_j$ (or a transformation $T(\theta_j)$) to the unknown parameter of interest. The vector of parameters of the network (weights and bias) is denoted $\omega$. During training, the goal is to find parameters of the network $\omega$ that are likely to have generated the outputs ($\theta_j, j = 1, \dots, N$), given the inputs ($\mathbf{x}_j, j = 1, \dots, N$). For this goal, we use another Bayesian inference model to infer such $\omega$ and we then aim to estimate the distribution $\pi(\omega | D)$. After training, the goal is then to predict the parameter of interest $\theta$ associated to a new $\mathbf{x}$.
an approximation of $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ is obtained as the average of the K predictions (see Eq. (7)), associated with an heuristic uncertainty which is here the variance, see Eq. (8). d) Using this approximation and associated uncertainty, the score of each sample of the calibration dataset is calculated. The conformal quantile is then the $\frac{\frac{(N_{cal}+1)(1-\delta)}{N_{cal}}}{N_{cal}}$ quantile of the calibration scores.

Then, for a new data sample $\mathbf{x}$, to approximate $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ and obtain an associated confidence interval: e) Monte Carlo Dropout prediction is performed to obtain an approximation of $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ (average of the predictions) and the associated heuristic uncertainty (variance of the predictions). A confidence interval is then calculated using the conformal quantile calculated in d) Eq. (11).

The pseudo-code of the algorithm proposed is detailed in Algorithm 2.

# 4 Applications

In this section, we will illustrate our model on three examples: The Moving Average 2 model (Section 4.1), a two-dimensional Gaussian random field (Section 4.2) and a Lokta-Volterra model (Section 4.3).

For all these examples, we compare the ABCD-Conformal algorithm with the standard ABC, the ABC-CNN of Akesson et al. [2022] and the ABC-RF of Raynal et al. [2018] (only for the two-dimensional Gaussian random field of Section 4.2 because it is not designed for multidimensional output).

We simulate a reference table (training set) of size $N_{train}$. For the ABC-CNN and the ABCD-conformal we need a validation set of size $N_{val}$, plus a calibration set of size $N_{cal}$ for the ABCD-conformal. We will test these methods on a test set of size $N_{test}$. The sizes used in the different examples are given in Table 1). These datasets are the same for the different methods studied. To approximate the posterior, the standard ABC and the ABC-CNN are used with a $\alpha$ acceptance ratio, hence $\alpha \times N_{train}$ generated samples are kept to approximate the posterior distribution.

Also, the type of results obtained with the different methods are not exactly the same. Using the standard ABC approach and the ABC-CNN, we obtain approximations of the whole posterior distributions. Using ABC-RF or ABCD-conformal approach, we obtain directly estimates of some posterior functional estimate of the parameter of interest $\theta$, in particular the posterior expected values $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ which can be used to estimate the true parameters $\theta$. To compare the methods, we then compute estimates of the posterior expected values from the approximated posterior distributions obtained with standard ABC and ABC-CNN. Performing that for the $N_{test}$ datasets of the test set, we can compare these estimates with the true values of the parameters used to generate these datasets. We compute NMAE as well as standard deviations of the absolute differences between the estimates and the true values of the parameters. More precisely, for a parameter of interest $\theta$, we compute

$$NMAE = \frac{\sum_{i=1}^{N_{test}} | \theta_i - \hat{\theta}_i |}{\sum_{i=1}^{N_{test}} \theta_i},$$

where $\theta_i$ is the true value of $\theta$ for the $i^{th}$ dataset, and $\hat{\theta}_i$ the estimate of $\theta$ for the $i^{th}$ dataset. Note that this NMAE is defined for unidimensional parameters. Hence, in case of parameters of interest that are multidimensional, we calculate it for each marginal of $\theta$ separately. For standard ABC and ABC-CNN, $\hat{\theta}_i$ corresponds to the empirical mean of the $\theta$ corresponding to the $\alpha N_{train}$ samples kept to approximate the posterior distribution of $\theta$ given the $i^{th}$ dataset. For ABCD-conformal, $\hat{\theta}_{1,i}$ is the empirical mean of $N_{val}$ stochastic forward passes through the network with Dropout, given the $i^{th}$ dataset as input (see Eq. (7)). To have an idea of the variability of the absolute errors, we compute $sd(|\theta - \hat{\theta}|)$ where

$$sd(|\theta - \hat{\theta}|) = \sqrt{\frac{1}{N_{test}} \sum_{i=1}^{N_{test}} (\theta - \hat{\theta}_i)^2}.$$

Finally, as we are interested by the uncertainty associated with an estimate, we also want to compare credible and confidence sets for the three methods. Concerning ABC-RF and ABCD-conformal, confidence sets (for
The posterior distribution of interest is then:

$$ \pi(\theta | \mathbf{x}, D) = \int \pi(\theta | \mathbf{x}, \omega)\pi(\omega | D) \mathrm{d}\omega. \qquad (5) $$

The density $\pi(\omega | D)$ is quite complex and cannot be evaluated analytically. It is approximated using a variational approach by $q(\omega)$. The approximate posterior distribution of interest is then given by

$$ q(\theta | \mathbf{x}) = \int \pi(\theta | \mathbf{x}, \omega)q(\omega) \mathrm{d}\omega. \qquad (6) $$

The first two moments of this distribution can be estimated empirically following Monte Carlo integration with K i.i.d. samples with law $q(\omega)\mathrm{d}\omega$. For each Monte Carlo sample $k$, different units of the network are randomly dropped out [Hinton et al., 2012]; we note $\hat{\omega}_k$ the estimated weights of the associated network with dropped units. Hence, these weights $\hat{\omega}_k$ are different for each Monte Carlo iteration. Denote by $\mathbf{f}^\omega(\mathbf{x})$ the neural network's stochastic output for input $\mathbf{x}$ and parameters $\omega$. If we assume that $\theta | \mathbf{x}, \omega \sim \mathcal{N}(\mathbf{f}^\omega(\mathbf{x}), \tau^{-1}\mathbf{I})$ and $\hat{\omega}_k \sim q(\omega)$, $k = 1, \dots, K$, an unbiased (consistent) estimator for $\mathbb{E}_{q(\theta|\mathbf{x})}[\theta | \mathbf{x}]$ is given by:

$$ \hat{\theta}(\mathbf{x}) = \frac{1}{K} \sum_{k=1}^{K} \mathbf{f}^{\hat{\omega}_k}(\mathbf{x}), \qquad (7) $$

which corresponds to the average of K stochastic forward passes through the network with Dropout. We refer the reader to the proposition 2 on page 47 of Yarin Gal's thesis [Gal, 2016] for more details about how to obtain this expression and the followings. An unbiased (consistent) estimator for the second moment $\mathbb{E}_{q(\theta|\mathbf{x})}[\theta^T\theta]$ is given by:

$$ \hat{\mathbb{E}}[\theta^T \theta | \mathbf{x}] = \tau^{-1} \mathbf{I} + \frac{1}{K} \sum_{k=1}^{K} \mathbf{f}^{\hat{\omega}_k}(\mathbf{x})^T \mathbf{f}^{\hat{\omega}_k}(\mathbf{x}) $$

Then to obtain a variance associated with the prediction of $\theta$ we can use the following unbiased (consistent) estimator:

$$ \hat{\mathbb{V}}[\theta | \mathbf{x}] = \underbrace{\tau^{-1} \mathbf{I}}_{\hat{\mathbb{V}}_a[\theta|\mathbf{x}]} + \frac{1}{K} \underbrace{\sum_{k=1}^{K} \mathbf{f}^{\hat{\omega}_k}(\mathbf{x})^T \mathbf{f}^{\hat{\omega}_k}(\mathbf{x})}_{\hat{\mathbb{V}}_e[\theta|\mathbf{x}]} - \hat{\theta}(\mathbf{x})^T \hat{\theta}(\mathbf{x}), \qquad (8) $$

which corresponds to the inverse model precision plus the sample variance of K stochastic forward passes through the network with Dropout. Hence, using Dropout and its associated randomness, we can obtain an estimate of $\mathbb{E}_{q(\theta|\mathbf{x})}[\theta]$, as well as associated uncertainty. This estimate can be considered as an approximation of $\mathbb{E}_\pi[\theta]$ because $q(\theta | \mathbf{x})$ is an approximation of the posterior of interest $\pi(\theta | \mathbf{x}, D)$. Moreover, Gal [2016] interprets the uncertainty Eq. (8) as the sum of an aleatoric uncertainty $\hat{\mathbb{V}}_a[\theta | \mathbf{x}]$ and an epistemic uncertainty $\hat{\mathbb{V}}_e[\theta | \mathbf{x}]$. The
**Algorithm 2:** Pseudo code for ABCD-conformal

**Input** : A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, integers $N_{\text{train}}$, $N_{\text{cal}}$ and $K$ representing sizes of training and calibration sets, and the number of stochastic passes for the Dropout. $\delta$ between 0 and 1 to obtain $(1 - \delta)\%$ confidence sets.

**Output:** An approximation of the posterior expected value $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ and a confidence interval for $\theta$.

1 a) Generation of a reference table (training dataset) and a calibration dataset :

2 These datasets are of sizes $N_{\text{train}}$, $N_{\text{val}}$ and $N_{\text{cal}}$ respectively. For the reference table for instance:

3 for $j \leftarrow 1$ to $N_{\text{train}}$ do

4 Draw $\theta_j \sim \pi$;

5 Draw synthetic sample $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^{\top}$ from the model $f(\cdot | \theta_j)$;

6 end

7 b) Train a CNN with concrete Dropout on the reference table: inputs are the $\{\mathbf{x}_j, j = 1, \dots, N_{\text{train}}\}$, and the outputs are the $\{\theta_j, j = 1, \dots, N_{\text{train}}\}$.

8 c) Monte Carlo Dropout prediction on the calibration set:

9 for $j \leftarrow 1$ to $N_{\text{cal}}$ do

10 for $k \leftarrow 1$ to $K$ do

11 $\mathbf{x}_j$ is given as input to the trained network with Dropout to obtain outputs $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x}_j)$ and $\tau_{\hat{\omega}_k}^{-1}$;

12 end

13 Obtain $\hat{\theta}_j$ by averaging the $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x}_j)$ outputs (see eq. (7))), $\tau^{-1}$ by averaging the $\tau_{\hat{\omega}_k}^{-1}$;

14 and an associated uncertainty $\hat{\mathbf{V}}(\mathbf{x}_j)$ that can be $\hat{\mathbf{V}}_a[\theta_j | \mathbf{x}_j], \hat{\mathbf{V}}_e[\theta_j | \mathbf{x}_j]$ or $\hat{\mathbf{V}}[\theta_j | \mathbf{x}_j]$ (see eq. (8)).

15 end

16 d) Computation of the conformal quantile on the calibration set :

17 for $j \leftarrow 1$ to $N_{\text{cal}}$ do

18 Compute the calibration score $s_j$ using the score function, see eq. (9):

$$s_j = \sqrt{(\theta_j - \hat{\theta}_j)^T \hat{\mathbf{V}}(\mathbf{x}_j)^{-1} (\theta_j - \hat{\theta}_j)},$$

19 end

20 Compute the conformal quantile $\hat{q}$ as the $\frac{[(N_{\text{cal}}+1)(1-\delta)]}{N_{\text{cal}}}$ quantile of the calibration scores $s_1, \dots, s_{N_{\text{cal}}}$.

21 e) For the new data sample $\mathbf{x}$, approximation of $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ and confidence set for $\theta$ :

22 for $k \leftarrow 1$ to $K$ do

23 $\mathbf{x}$ is given as input to the trained network, to obtain an output $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$;

24 end

25 Obtain $\hat{\theta}(\mathbf{x})$ an approximation of $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ by averaging these outputs (see eq. (7)), and an associated uncertainty $\hat{\mathbf{V}}(\mathbf{x})$ (using eq. (8));

26 The confidence set for $\theta$ is an ellipsoid which center is $\hat{\theta}(\mathbf{x})$ and covariance matrix is $\hat{\mathbf{V}}(\mathbf{x})^{-1}/\hat{q}^2$, see eq. (10).
aleatoric uncertainty $\hat{V}_a[\theta | \mathbf{x}]$ is interpreted as the noise in the data, it is the result of measurement imprecision (it is often modelled as part of the likelihood, and this is often a Gaussian corrupting noise). The epistemic uncertainty $\hat{V}_e[\theta | \mathbf{x}]$ is interpreted as the model uncertainty: uncertainty in the model parameters or in the structure of the model. The total predictive uncertainty $\hat{V}[\theta | \mathbf{x}]$ combines these both types of uncertainties and will be used as our heuristic uncertainty $\hat{V}(\mathbf{x})$. In practice, for each of the K Monte Carlo iterations, the weights $\hat{\omega}_k$ are different, and the NN with input $\mathbf{x}$ gives as outputs $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$ and $\tau_{\hat{\omega}_k}^{-1}$. The aleatoric uncertainty $\tau^{-1}$ is estimated by the mean of the $\tau_{\hat{\omega}_k}^{-1}$, and the epistemic uncertainty by the sample variance of the $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$.

The main drawback of Dropout layers is that the Dropout probability is a new model hyperparameter. To circumvent this issue, Gal et al. [2017] propose a Dropout variant, called Concrete Dropout, that allows for automatic tuning of this probability in large models, using gradient methods. It improves performance and tunes automatically the Dropout rate, producing better calibrated uncertainties

**Remark 1.** Gal [2016] and Gal and Ghahramani [2016] showed that optimizing any neural network with Dropout is equivalent to a form of approximate inference in a probabilistic interpretation of the model. In other words, the optimal parameters found through the optimization of a Dropout neural network are the same as the optimal variational parameters in a probabilistic Bayesian neural network with the same structure (the parameters are found through a variational inference approach). It means that a network trained with Dropout is equivalent to a Bayesian Neural Network and possesses all the properties of such a network (see MacKay [1992] or Neal [2012] concerning Bayesian Neural Networks). For this equivalence to be true, only one condition should be verified in the variational inference approach, which is about the Kullback-Leibler divergence between the prior distribution of the parameters $\omega$ and an approximating distribution for $\omega$. We will not need to verify such assumption because the theoretical guaranty of our confidence sets will be assured through the conformal prediction.

## 3.3 Conformal prediction

For a test input $\mathbf{x}$, the NN with concrete Dropout method enables us to obtain an approximation $\hat{\theta}(\mathbf{x})$ of $E_\pi[\theta | \mathbf{x}]$, as well as associated heuristic uncertainties $\hat{V}(\mathbf{x})$. However, it does not give valid confidence sets. In order to obtain such confidence sets, we propose to apply a conformal procedure, as explained in Angelopoulos and Bates [2023]. Conformal prediction is indeed a straightforward way to generate valid confidence sets for a wide variety of models from a proxy of the uncertainty. The better this proxy will be, the smaller the confidence interval will be. This procedure requires a level $\delta \in [0, 1]$ and a small amount of additional calibration data $((\theta_j, \mathbf{x}_j), j = 1, \dots, N_{\text{cal}})$, hundreds of calibration data are theoretically enough (see Eq. (12) below). We propose to use as heuristic uncertainty $\hat{V}(\mathbf{x})$ one of the variances obtained from the NN with Dropout in Section 3.2, see Eq. (8). This heuristic uncertainty can be transformed into a rigorous confidence interval $C(\mathbf{x})$ (with confidence level $1 - \delta$) through the following conformal procedure.
marginals of $\theta$) are directly given by the proposed method and we used $\delta = 0.05$. Two uncertainty measures were used for comparison in ABCD-conformal: the overall uncertainty and the epistemic uncertainty returned by the Dropout procedure. Concerning standard ABC and ABC-CNN, we used the $\delta/2$ and $1 - \delta/2$ quantiles of the approximated posterior distributions to obtain credible intervals. We then compare the coverage of these credible and confidence sets for these methods, as well as the mean lengths of these sets. The coverage is estimated by the number of true values of marginals of $\theta$ in the obtained sets, for the $N_{test}$ simulations from the test set.

All the parameters used in the different examples are synthesized in Table 1.

<table><caption>Table 1: Sizes of datasets used for the different methods: $N_{train}$ for the training set, $N_{val}$ for the validation set, $N_{cal}$ for the calibration set and $N_{test}$ for the test set, and $\alpha$ the tolerance threshold used for the ABC standard and ABC-CNN.</caption><thead><tr><th>Example</th><th>N<sub>train</sub></th><th>N<sub>val</sub></th><th>N<sub>cal</sub></th><th>N<sub>test</sub></th><th>α</th></tr></thead><tbody><tr><td>Section 4.1</td><td>10<sup>4</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0.01</td></tr><tr><td>Section 4.2</td><td>7 × 10<sup>3</sup></td><td>1900</td><td>10<sup>3</sup></td><td>10<sup>2</sup></td><td>0.01</td></tr><tr><td>Section 4.3</td><td>10<sup>5</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0.005</td></tr></tbody></table>

# 4.1 The Moving Average 2 model

## 4.1.1 The model

The Moving Average 2 model (MA(2)) is a simple benchmark example used in Bayesian and ABC literature (Marin and Robert [2007], Marin et al. [2012], Jiang et al. [2017], Wiqvist et al. [2019] or Akesson et al. [2022] for instance). It is defined for observations $X_1, \dots, X_p$;

$$X_j = Z_j + \theta_1 Z_{j-1} + \theta_2 Z_{j-2}, j = 1, \dots, p,$$

where $(Z_j)_{-2<j\le p}$ is an i.i.d. sequence of standard Gaussians $N(0, 1)$. This model is identifiable in the following triangular region $\mathcal{D}$:

$$-2 < \theta_1 < 2, \quad \theta_1 + \theta_2 > -1, \quad \theta_1 - \theta_2 < 1;$$

See for instance Marin and Robert [2007][Chapter 5] or [Marin et al., 2012, Section 2]. We assume that the prior distribution on $(\theta_1, \theta_2)$ is the uniform distribution over this triangular region. This model allows for exact calculation of the posterior $(\theta_1, \theta_2) | X$.

The goal is to obtain posterior estimates of the bi-dimensional vector of parameters $(\theta_1, \theta_2)$ from an observed dataset of length $p$, $(X_1, \dots, X_p)$. In the following, we will use $p = 100$.

## 4.1.2 Algorithm parametrization

Number of used samples are resumed in Table 1. For the standard ABC, the summary statistics chosen are the first two autocovariances:

$$\tau_1 = \sum_{j=2}^{p} x_j x_{j-1} \quad \text{and} \quad \tau_2 = \sum_{j=3}^{p} x_j x_{j-2}.$$

To compare two samples, we then use the $L^2$ distance between the two associated vector of summary statistics.

In ABC-CNN method, the distance function between a sample from the test set and a sample from the training set is the quadratic distances $L^2$ distance between the parameters predicted by the CNN for the test sample, and true parameters used to simulate the training sample.

Concerning ABC-CNN and ABCD-conformal, architectures of neural networks are the same for the two methods: 3 convolutional 1D layers with 64 neurons and a kernel size of 3 followed by max-pooling for the
Using the heuristic uncertainty $\hat{V}(\mathbf{x})$, we compute a score $s(\theta, \mathbf{x})$ for each data sample in the calibration set, from these scores we can define $\hat{q}$ the conformal quantile and then a confidence set. More precisely, as in Messoudi et al. [2022], the score function will be

$$s(\theta, \mathbf{x}) = \sqrt{(\theta - \hat{\theta}(\mathbf{x}))^t \hat{V}(\mathbf{x})^{-1} (\theta - \hat{\theta}(\mathbf{x}))}. \quad (9)$$

We can then define $\hat{q}$ as the $\frac{[(N_{cal}+1)(1-\delta)]}{N_{cal}}$ quantile of $\{s(\theta_j, \mathbf{x}_j), j = 1, \dots, N_{cal}\}$ to have as the confidence set:

$$C(\mathbf{x}) = \{\theta \mid s(\theta, \mathbf{x}) \le \hat{q}\}. \quad (10)$$

This corresponds to an ellipsoid whose center is $\hat{\theta}(\mathbf{x})$ and whose covariance matrix is $\hat{V}(\mathbf{x})^{-1}/\hat{q}^2$. In one dimension, we have $s(\theta, \mathbf{x}) = |\theta - \hat{\theta}(\mathbf{x})|/\sqrt{\hat{V}(\mathbf{x})}$, and

$$C(\mathbf{x}) = \left[ \hat{\theta}(\mathbf{x}) - \hat{q}\sqrt{\hat{V}(\mathbf{x})}; \hat{\theta}(\mathbf{x}) + \hat{q}\sqrt{\hat{V}(\mathbf{x})} \right]. \quad (11)$$

This confidence set satisfies the marginal coverage property:

$$1 - \delta \le \mathbb{P}[\theta \in C(\mathbf{x})] \le 1 - \delta + \frac{1}{N_{cal} + 1}, \quad (12)$$

as shown in Angelopoulos and Bates [2023, Appendix D] (in particular the right-hand side inequality needs that random variables $s(\theta_j, \mathbf{x}_j)$ to be continuous). This probability is marginal over the randomness in the calibration and test points [Angelopoulos and Bates, 2023], see Algorithm 2. One great advantage of this procedure is that it is easy to implement, fast and generic. Moreover, the obtained prediction intervals are non-asymptotic, with no distributional or model assumptions. The only assumption needed is that the calibration data and the test data are i.i.d [Vovk et al., 1999] which is naturally the case in our setting. Finally, this procedure provides a frequentist perspective on a Bayesian approach, reconciling the two points of view.

The performance of this approach in terms of intervals length depends only on the quality of the uncertainty measure used. A “good” uncertainty measure should reflect the magnitude of model error: it should be smaller for easy input, and larger for hard ones. Hard inputs here can refer to extreme points in the training set e.g. see example in section 4.2 where some $\theta_3$ can be very large, and/or with high aleatoric uncertainty. If the uncertainty measure is “bad”, then the intervals obtained by the conformal procedure will be quite large and their lengths will not differ a lot between easy and hard inputs. From the previous MC Dropout method, three types of uncertainties were obtained: epistemic, aleatoric and overall (sum of aleatoric and epistemic).

**Remark 2.** Note that a conditional coverage would be preferable, that is:

$$\mathbb{P}[\theta \in C(\mathbf{x}) | \mathbf{x}] \ge 1 - \delta. \quad (13)$$

Indeed, marginal coverage does not guarantee the coverage of $C(\mathbf{x})$ for any specific $\mathbf{x}$, only the average coverage over the whole domain. However, this conditional property is not guaranteed by the conformal procedure in general, but some metrics can be used to check how close we are to this stronger property, see Angelopoulos and Bates [2023].
first two layers and by a flatten layer for the third one. Then, 3 dense layers of 100 neurons. Each of the layers use the `relu` activation function. The raw samples of length $p = 100$ are the inputs of the neural networks, and the outputs are the bi-dimensional associated vectors of parameters $(\theta_1, \theta_2)$.

### 4.1.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 2. Figure 1 shows the predicted parameters $\hat{\theta}_1$ and $\hat{\theta}_2$ estimated by the ABCD-conformal method against the true values for the 1000 test samples, and confidence sets using gray intervals.

To highlight the sharpness of the marginal coverage from the ABCD-Conformal algorithm, we have repeated the experiment with 10 different calibration and test sets, and we obtain (in mean) coverages of 94.76% for $\theta_1$ and 94.68% for $\theta_2$.

<table><thead><tr><th></th><th>Standard ABC</th><th>ABC-CNN</th><th>ABCD-conformal<br>overall</th><th>ABCD-conformal<br>epistemic</th></tr></thead><tbody><tr><td>NMAE<sub>1</sub></td><td>0.1852</td><td>0.1963</td><td>0.2008</td><td>0.2008</td></tr><tr><td>NMAE<sub>2</sub></td><td>0.2644</td><td>0.2709</td><td>0.2627</td><td>0.2627</td></tr><tr><td>sd(|θ<sub>1</sub> - θ̂<sub>1</sub>|)</td><td>0.0953</td><td>0.1001</td><td>0.0994</td><td>0.0994</td></tr><tr><td>sd(|θ<sub>2</sub> - θ̂<sub>2</sub>|)</td><td>0.1061</td><td>0.1098</td><td>0.1046</td><td>0.1046</td></tr><tr><td>mean length confidence sets θ<sub>1</sub></td><td>0.6003</td><td>0.2889</td><td>0.7001</td><td>0.6368</td></tr><tr><td>mean length confidence sets θ<sub>2</sub></td><td>0.6385</td><td>0.2891</td><td>0.7771</td><td>0.6621</td></tr><tr><td>coverage confidence sets θ<sub>1</sub></td><td>94.8%</td><td>62.1%</td><td>94.8%</td><td>93.0%</td></tr><tr><td>coverage confidence sets θ<sub>2</sub></td><td>93.7%</td><td>61.2%</td><td>94.4%</td><td>94.6%</td></tr></tbody></table>

Table 2: Moving Average 2 example: comparison of standard ABC, ABC-CNN and ABCD-conformal (using overall or epistemic uncertainty as uncertainty heuristic for the conformal procedure). These indicators are computed on a test set of size $10^3$.

We can see in Table 2 that the three methods obtain quite similar results concerning NMAE and standard deviation of the absolute error. Some differences can be noted for the coverages, and the mean lengths of the confidence sets. The coverage is sharp for standard ABC and ABCD-conformal. However, it is too small for ABC-CNN. This is because the approximated posteriors were too peaked using this method, the dispersions around the mean values were too small, leading to too narrow (hence useless) confidence sets. Concerning the mean lengths, ABC-CNN has the smallest ones, but as the confidence sets are too narrow we can not say that this method outperforms the others. The standard ABC gives the best results in terms of mean lengths, followed closely by ABCD-conformal using the epistemic uncertainty, then by ABCD-conformal using the overall uncertainty. Concerning the particularly good results of the standard ABC, it is not surprising, as we used good summary statistics. Moreover, as showed by Frazier et al. [2018], conditions for good asymptotic properties are verified on this example, which is quite rare.

To sum up, ABC-CNN does not provide sufficient guarantees for these confidence sets, and ABCD-Conformal is performing as well as the standard ABC, which is here in favourable estimation conditions.

## 4.2 Two-dimensional Gaussian random field

### 4.2.1 The model

We study stationary isotropic Gaussian random fields on the domain $[0, 5] \times [0, 5]$ with a regular grid of size $100 \times 100$, with exponential covariance functions [Wood and Chan, 1994]. This covariance between two points $z_1$ and $z_2$ is given by:

$$C(z_1, z_2) = \exp \left( \left| \frac{z_1 - z_2}{\theta} \right|^2 \right),$$

with $\theta$ the range (or scale) parameter, which is the unknown parameter. We assume that the prior distribution on $\theta$ is the uniform distribution between 0 and 1.
## 3.4 Implementation of ABCD-Conformal

The training step of the proposed method is as follows. a) Generate the training dataset i.e., reference table in ABC vocabulary and the calibration dataset, of sizes $N_{train}$ and $N_{cal}$ respectively. b) The NN with Dropout is trained using the reference table. In practice the aleatoric uncertainty $\hat{V}_a$ is learned while training the neural network with an heteroscedastic loss as in [Kendall and Gal, 2017, section 3]. One could select the network architecture using an extra validation dataset. We refer the reader to the many tutorials and books available on the implementation of neural networks. c) Monte Carlo Dropout prediction is performed, that is, each sample of the calibration dataset is passed through the trained network with Dropout, K times. For each of these samples, an approximation of $E_\pi[\theta | \mathbf{x}]$ is obtained as the average of the K predictions (see Eq. (7)), associated with an *heuristic uncertainty* which is here the variance, see Eq. (8). d) Using this approximation and associated uncertainty, the score of each sample of the calibration dataset is calculated. The conformal quantile is then the $\frac{\left[ (N_{cal}+1)(1-\delta) \right]}{N_{cal}}$ quantile of the calibration scores. Then, for a new data sample $\mathbf{x}$, to approximate $E_\pi[\theta | \mathbf{x}]$ and obtain an associated confidence interval: e) Monte Carlo Dropout prediction is performed to obtain an approximation of $E_\pi[\theta | \mathbf{x}]$ (average of the predictions) and the associated heuristic uncertainty (variance of the predictions). A confidence interval is then calculated using the conformal quantile calculated in d) Eq. (11).

The pseudocode of the algorithm proposed is detailed in Algorithm 2.

# 4 Applications

In this section, we will illustrate our method on three examples: The Moving Average 2 (MA2) model (Section 4.1), a Lotka-Volterra model (Section 4.2) and a model for lake ecosystem dynamics (Section 4.3). An additional example about Gaussian random field is given in the appendix D. For all these examples, we compare the performances of ABCD-Conformal with those of standard ABC, ABC-SMC, ABC-CNN (Akesson et al. [2022]) and ABC-RF (Raynal et al. [2018]) when possible.

**Datasets.** For each example, we simulate a reference table (training set) of size $N_{train}$. For ABC-CNN and ABCD-Conformal we need a validation set of size $N_{val}$, plus a calibration set of size $N_{cal}$ for ABCD-Conformal. We will test these methods on a test set of size $N_{test}$. The sizes used in the different examples are given in Table 1). These datasets are the same for the different methods studied. To approximate the posterior, standard ABC and ABC-CNN are used with an $\alpha$ acceptance ratio, hence $\alpha \times N_{train}$ generated samples are kept to approximate the posterior distribution. Concerning ABC-SMC, the approach is not based on a reference table, hence we can not use the same reference table as for the other methods. Moreover, several choices must be made to implement the algorithm, such as how to decrease the sequence of tolerance levels, or how to select the transition kernel used to perturb a random draw from a previous population. Our choices can be found at the end in Appendix E.1.
<figure><img src="image_3.png" /><figcaption>Figure 1: Moving Average 2 example: predicted parameters $\hat{\theta}_1$ and $\hat{\theta}_2$ estimated by the ABCD-conformal against the true values for the 1000 test samples, and confidence sets obtained by the conformal procedure (in gray intervals, for the epistemic uncertainty).</figcaption></figure>

### 4.2.2 Algorithm parametrization

Number of used samples are resumed in Table 1. Concerning the standard ABC and the ABC-RF, we use two summary statistics: the Moran's *I* statistics from lag 1 to 5 (that is the Moran's correlogram from lags 1 to 5), and the semi-variogram up to a distance of 20 (15 values kept per variogram) (see Cliff and Ord [1981] for these notions). We assume the spatial weights matrix to be row-standardised, and the neighbors of a pixel being the 4-nearest pixels. As all Gaussian fields are simulated on the same grid, the variograms of the different fields are calculated at exactly the same distances, so they are comparable. For the standard ABC, the distance used to compare two Gaussian random fields is the sum of the quadratic distance between their Moran's correlograms and of the quadratic distance between their semi-variograms.

Concerning the architectures of neural networks for ABC-CNN and ABC-Dropout, we used 3 convolutional layers having 32, 64 and 64 neurons and the relu activation function, followed by 2 dense layers having both 64 neurons and the relu activation function. The raw samples given as inputs of the neural networks are the 2D Gaussian random fields, and the associated outputs are the estimations of the parameter $\theta$.

### 4.2.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 3. Figure 2 shows the predicted parameters $\hat{\theta}$ estimated by the ABCD-conformal method against the true values for the 100 test samples, and using gray intervals.

In Table 3, we see that on this example, ABC-RF outperforms the others methods, with the smallest NMAE and $sd(|\theta - \hat{\theta}|)$, while having a very good coverage, 99%, which is better than the expected coverage of 95%. The standard ABC and ABC-CNN do not give satisfactory results. Indeed, while having a coverage of 100%, the standard ABC has the largest NMAE and $sd(|\theta - \hat{\theta}|)$, and the good coverage is due to too large coverage. Concerning the ABC-CNN, like in the previous example, the NMAE, $sd(|\theta - \hat{\theta}|)$ and the mean lengths of the confidence sets are quite small, but it is counterbalanced by a too small coverage (88%). Finally, the
**Algorithm 2:** Pseudocode for ABCD-Conformal

**Input** : A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, integers $N_{\text{train}}$, $N_{\text{cal}}$ and $K$ representing sizes of training and calibration sets, and the number of stochastic passes for the Dropout. $\delta$ between 0 and 1 to obtain $(1 - \delta)$ confidence sets.

**Output:** An approximation of the posterior expected value $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ and a confidence interval for $\theta$.

1. **a) Generation of a reference table (training dataset) and a calibration dataset:**

2. For the reference table for instance:

3. **for** $j \leftarrow 1$ **to** $N_{\text{train}}$ **do**

4. &nbsp;&nbsp;&nbsp;&nbsp;Draw $\theta_j \sim \pi$;

5. &nbsp;&nbsp;&nbsp;&nbsp;Draw synthetic sample $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^{\top}$ from the model $f(\cdot | \theta_j)$;

6. **end**

7. **b) Train a Neural Network with concrete Dropout on the reference table:** the training pairs are the {($\mathbf{x}_j, \theta_j$), $j = 1, \dots, N_{\text{train}}$} and the loss is the heteroscedastic loss as in [Kendall and Gal, 2017, Eq. (8)].

8. **c) Monte Carlo Dropout prediction on the calibration set:**

9. **for** $j \leftarrow 1$ **to** $N_{\text{cal}}$ **do**

10. &nbsp;&nbsp;&nbsp;&nbsp;**for** $k \leftarrow 1$ **to** $K$ **do**

11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**x**<sub>j</sub> is given as input to the trained network to obtain outputs $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x}_j)$ and $\tau_{\hat{\omega}_k}^{-1}$;

12. &nbsp;&nbsp;&nbsp;&nbsp;**end**

13. &nbsp;&nbsp;&nbsp;&nbsp;Obtain $\hat{\theta}_j$ and $\tau^{-1}$ by averaging the $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x}_j)$ outputs (see eq. (7))) and the $\tau_{\hat{\omega}_k}^{-1}$;

14. &nbsp;&nbsp;&nbsp;&nbsp;and an associated uncertainty $\hat{V}(\mathbf{x}_j)$ that can be $\hat{V}[\theta_j | \mathbf{x}_j]$ (see eq. (8)).

15. **end**

16. **d) Computation of the conformal quantile on the calibration set:**

17. **for** $j \leftarrow 1$ **to** $N_{\text{cal}}$ **do**

18. &nbsp;&nbsp;&nbsp;&nbsp;Compute the calibration score $s_j$ using the score function, see eq. (9)

$$s_j = \sqrt{(\theta_j - \hat{\theta}_j)^t \hat{V}(\mathbf{x}_j)^{-1} (\theta_j - \hat{\theta}_j)},$$

19. **end**

20. Compute the conformal quantile $\hat{q}$ as the $\frac{[(N_{\text{cal}}+1)(1-\delta)]}{N_{\text{cal}}}$ quantile of the calibration scores $s_1, \dots, s_{N_{\text{cal}}}$.

21. **e) For the new data sample x, approximation of E<sub>π</sub>[θ | x] and confidence set for θ:**

22. **for** $k \leftarrow 1$ **to** $K$ **do**

23. &nbsp;&nbsp;&nbsp;&nbsp;**x** is given as input to the trained network, to obtain an output $\mathbf{f}^{\hat{\omega}_k}(\mathbf{x})$;

24. **end**

25. Obtain $\hat{\theta}(\mathbf{x})$ an approximation of $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ by averaging these outputs (see Eq. (7)), and an associated uncertainty $\hat{V}(\mathbf{x})$ (using Eq. (8));

26. The confidence set for $\theta$ is an ellipsoid whose center is $\hat{\theta}(\mathbf{x})$ and covariance matrix is $\hat{V}(\mathbf{x})^{-1}/\hat{q}^2$, see Eq. (10).
<table><caption>Table 3: Comparison of standard ABC, ABC-CNN, ABC-RF and ABCD-conformal on the two-dimensional Gaussian random field 2 example. These indicators are computed on a test set of size 10².</caption><thead><tr><th></th><th>Standard<br>ABC</th><th>ABC-CNN</th><th>ABC-RF</th><th>ABCD-conformal<br>overall</th><th>ABCD-conformal<br>epistemic</th></tr></thead><tbody><tr><td>NMAE</td><td>0.0436</td><td>0.0223</td><td>0.0154</td><td>0.0303</td><td>0.0303</td></tr><tr><td>sd(|θ - θ̂|)</td><td>0.0229</td><td>0.0106</td><td>0.0077</td><td>0.0141</td><td>0.0141</td></tr><tr><td>coverage confidence sets θ</td><td>100%</td><td>88%</td><td>99%</td><td>93%</td><td>95%</td></tr><tr><td>mean length confidence sets θ</td><td>0.1313</td><td>0.0307</td><td>0.0414</td><td>0.0684</td><td>0.0718</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: Predicted parameters $\hat{\theta}$ estimated by the ABCD-conformal against the true values for the 100 test samples, and (in gray intervals, for the overall uncertainty).</figcaption></figure>

results given by the ABCD-conformal approach are sharp, with coverages of exactly 95% (for the epistemic uncertainty). The overall and the epistemic uncertainty give quite similar results in the conformal procedure.

To sum up, the ABC-RF is efficient here on all criteria. ABCD-Conformal and ABC-CNN give good predictions, but ABCD-Conformal is much better on . Standard ABC is unsatisfactory on all criteria.

## 4.3 Lokta-Volterra model

### 4.3.1 The model

The Lokta-Volterra model describes the dynamics of biological systems in which two species interact, one as a predator and the other as a prey. Here we consider a stochastic Markov jump process version of this model with state $(X_1, X_2) \in \mathbb{Z}^2$ representing prey and predator population sizes, and three parameters $c_1, c_2$ and
**Evaluation.** The type of results obtained with the different methods are not exactly the same. Using the standard ABC approach, the ABC-SMC and the ABC-CNN, we obtain approximations of the whole posterior distributions. Using ABC-RF or ABCD-Conformal approaches, we obtain directly estimates of some posterior functional estimate of the parameter of interest $\theta$, in particular the posterior expected value $\mathbb{E}_{\pi}[\theta | \mathbf{x}]$ which can be used to estimate the true parameter $\theta$. To compare the methods, we then compute an estimate of the posterior expected value from the approximated posterior distribution obtained with standard ABC, ABC-CNN and ABC-SMC. Performing that for the $N_{test}$ datasets of the test set, we can compare the obtained estimates with the true values of the parameters used to generate these datasets. We compute NMAE as well as standard deviations of the absolute differences between the estimates and the true values of the parameters. More precisely, for a parameter of interest $\theta$, we compute

$$NMAE = \frac{\sum_{i=1}^{N_{test}} | \theta_i - \hat{\theta}_i |}{\sum_{i=1}^{N_{test}} | \theta_i |},$$

where $\theta_i$ is the true value of $\theta$ for the $i^{th}$ dataset, and $\hat{\theta}_i$ the estimate of $\theta$ for the $i^{th}$ dataset.

Let us give few comments on this error estimation. First, since the ABC method aims to provide a posterior mean estimator, it may seem odd to compare its output with the true parameter value rather than with the posterior mean, as these two quantities can differ substantially. However, in the (non-trivial) examples we consider, the exact posterior means are unknown and cannot be used as a reference for comparison. Moreover, the likelihoods associated to our examples are based on a large number of observations ($p = 100$ for the MA2 model, a grid of size $100 \times 100$ for the Gaussian random field, $36 \times 2$ observation points for the Lotka-Volterra model and 365 observation points for the Lake model), and so results such as the Bernstein-von Mises Theorem [Van der Vaart, 1998], or its generalizations, ensure that in our examples the true parameter value and the posterior mean are indeed close. Note, however, that for the MA2 example only, the exact posterior is explicit (although costly to evaluate). We therefore numerically compared the posterior mean with the true value of the parameter in our Julia Notebook (see supplementary materials) to confirm this claim. Figure A.2 in Appendix A illustrates this comparison for the 1000 samples in the test set.

Note that this NMAE is defined for unidimensional parameters. Hence, in case of parameters of interest that are multidimensional, we calculate it for each marginal of $\theta$ separately. For standard ABC and ABC-CNN, $\hat{\theta}_i$ corresponds to the empirical mean of the $\theta$ corresponding to the $\alpha N_{train}$ samples kept to approximate the posterior distribution of $\theta$ given the $i^{th}$ dataset. For ABC-SMC, it corresponds to the weighted mean of the last population. For ABCD-Conformal, $\hat{\theta}_i$ is the empirical mean of $N_{val}$ stochastic forward passes through the network with Dropout, given the $i^{th}$ dataset as input (see Eq. (7)). To have an idea of the variability of the absolute errors, we compute $sd(|\theta - \hat{\theta}|)$ where

$$sd(|\theta - \hat{\theta}|) = \sqrt{\frac{1}{N_{test}} \sum_{i=1}^{N_{test}} (\theta_i - \hat{\theta}_i)^2}.$$
$c_3$. Three transitions are possible, with hazard rates $c_1X_1$, $c_2X_1X_2$ and $c_3X_2$ respectively:

$$ (X_1, X_2) \xrightarrow{c_1 X_1} (X_1 + 1, X_2) \quad (\text{prey growth}) $$

$$ (X_1, X_2) \xrightarrow{c_2 X_1 X_2} (X_1 - 1, X_2 + 1) \quad (\text{predation interaction}) $$

$$ (X_1, X_2) \xrightarrow{c_3 X_2} (X_1, X_2 - 1) \quad (\text{predator death}) $$

The estimation of the parameters of this Lokta Volterra model has been studied by several authors in an ABC framework, see for instance Prangle [2017]. The initial conditions are $X_1(0) = 50$ and $X_2(0) = 100$ and a dataset corresponds to observations of state $(X_1, X_2)$ at times 0, 2, 4, ..., 36. As usual, all simulations with an extinction of either the preys or the predators were discarded: we are interested in the conditional law of survival. The prior distributions are independent uniforms $U[-6, 2]$ for the transformed parameters $\log(c_1)$, $\log(c_2)$ and $\log(c_3)$.

We are interested here in the estimation of $\theta = (c_1, c_2, c_3)$.

### 4.3.2 Algorithm parametrization

Number of used samples are resumed in Table 1. The Gillespie's stochastic simulation algorithm is used to generate all the samples, using the Explicit tau-leap method from the package GillespieSSA2 [Cannoodt et al., 2021]. To improve performances of the algorithm, once simulations have been done, we used standardised versions of the parameters. The goal is then to obtain posterior estimates of the tri-dimensional vector of normalized parameters $(c_1, c_2, c_3)$, from an observed sample consisting of two time series of sizes 19.

Concerning the standard ABC, there is no obvious summary statistics. Instead, the distance function $d((x_1, x_2), (x_{1d}, x_{2d}))$ between a sample from the test set $\{(x_{1d}[i], x_{2d}[i]), i = 1, \dots, 19\}$ and a sample from the training set $\{(x_1[i], x_2[i]), i = 1, \dots, 19\}$, is given by the sum of squared differences:

$$ d((x_1, x_2), (x_{1d}, x_{2d})) = \sum_{i=1}^{19} ((x_1[i] - x_{1d}[i])^2 + (x_2[i] - x_{2d}[i])^2) $$

Concerning ABC-CNN, the distance function between a sample from the test set and a sample from the training set is the quadratic distances between the parameters predicted by the CNN for the test sample, and the true parameters used to simulate the training sample.

The architectures used for the neural networks are the same for the ABC-CNN and the ABCD-conformal: 3 convolutional 1D layers with 128 neurons and a kernel size of 2, followed by max-pooling for the first two layers, and by a flatten layer for the third one. Then, 3 dense layers of 100 neurons. Each of the layers use the tanh activation function. The raw samples consisting of two time series of lengths 19 are the inputs of the neural networks, and the outputs are the tri-dimensional associated vectors of parameters $(c_1, c_2, c_3)$.

### 4.3.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 4. Figure 3 shows the predicted normalized parameters and associated confidence sets obtained with standard ABC, ABC-CNN and ABCD-conformal methods respectively for the 1000 test samples, and associated confidence sets.

We can see in Table 4 that ABCD-conformal outperforms standard ABC and obtains quite better results than ABC-CNN, when looking at NMAE and standard deviations of the absolute errors. Predictions of ABC-CNN and ABCD-conformal are better than those of the standard ABC, thanks to the power of CNN, and because we do not have relevant summary statistics for the standard ABC. Thanks to Dropout, the predictions are even better for ABCD-conformal than for ABC-CNN.

The confidence sets coverages appear similar for the three methods. But the mean lengths of confidence sets are different depending on the method and on the parameter. Using this criterion, ABCD-conformal is the best, followed by ABC-CNN then by standard ABC, for the three parameters. In this example, there is also an impact of the heuristic uncertainty measure used for the conformal procedure: the overall variance gives
**Confidence and credible Intervals/Regions.** Finally, as we are interested by the uncertainty associated with an estimate, we compare credible and confidence sets for the five methods. First, we compare confidence and credible sets for each component of $\theta$, hence in dimension one. For standard ABC and ABC-CNN, we used the $\delta/2$ and $1 - \delta/2$ quantiles of the marginals of the approximated posterior distributions to obtain credible intervals, for $\delta = 0.05$. For ABC-SMC we did the same thing, but using weighted quantiles. ABC-RF directly gives confidence intervals, for unidimensional parameters. Concerning ABCD-Conformal, we computed confidence intervals for each component of $\theta$, by using the diagonal terms of the uncertainty measure $\hat{V}(\mathbf{x})$ (by analogy with a variance-covariance matrix, this corresponds to using only the variance terms). Then, we compare confidence and credible sets in multidimensional terms. Confidence sets (ellipsoids) are directly obtained from the ABCD-Conformal approach. For ABC methods approximating the full multidimensional posterior with neighboring samples, i.e. ABC-Standard, ABC-SMC and ABC-CNN, one could define a 95% credible set using the Highest Posterior Density Region [Hyndman, 1996]. Here for simplicity and computational tractability, we assume the posterior to be a Multivariate Normal distribution. The associated credible set is the ellipsoid generated by Eqs. (9, 10) where $\hat{V}(\mathbf{x})$ is replaced by the covariance of the neighboring samples and $\hat{q}$ is replaced by the 0.95 quantile of a Chi distribution whose degrees of freedom are the dimension of $\theta$.

The default uncertainty measure used by ABCD-conformal in all the examples is the overall uncertainty returned by the Dropout procedure. For comparison, we also present results using the epistemic uncertainty for the Gaussian random field example (Appendix D), and results using a constant uncertainty for the Lotka Voltera example (Section 4.2).

For all methods, the coverage is estimated by the number of true values of $\theta$ or marginals of $\theta$ in the obtained sets, for the $N_{\text{test}}$ simulations from the test set. We then compare the coverage of these credible and confidence sets for these methods, the mean and median lengths of credible and confidence intervals for marginals of $\theta$, and volumes for multidimensional credible and confidence sets when possible.

**Code and Software** Source codes and instructions to run the examples, in R and/or Julia, are available as supplementary materials. For some examples, datasets and pretrained models are provided to reproduce exactly our results. For more details, see the dedicated section at the end of the article. For this paper, the authors developed a Julia package ABCMethods.jl [Métivier, 2025] where the ABC-Standard, ABC-SMC, ABC-CNN and ABC-Conformal methods and associated tools are coded. All the parameters, methods and languages used in the different examples are synthesized in Table 1.

# 4.1 The Moving Average 2 model

## 4.1.1 The model

The Moving Average 2 model (MA2) is a simple benchmark example used in Bayesian and ABC literature (see Marin and Robert [2007], Marin et al. [2012], Jiang et al. [2017], Wiqvist
<figure><img src="image_5.png" /><figcaption>Figure 3: Lokta Volterra example: predicted normalized parameters $\hat{c}_1$, $\hat{c}_2$ and $\hat{c}_3$ estimated by the standard ABC, the ABC-CNN and the ABCD-conformal methods, against the true values for the 1000 test samples. In gray are represented, for the overall uncertainty.</figcaption></figure>

<table><caption>Table 4: Lokta Volterra example: comparison of standard ABC, ABC-CNN and ABCD-conformal (using overall or epistemic uncertainty as uncertainty heuristic for the conformal procedure). These indicators are computed on a test set of size $10^3$.</caption><thead><tr><th></th><th>Standard ABC</th><th>ABC-CNN</th><th>ABCD-conformal<br>overall</th><th>ABCD-conformal<br>epistemic</th></tr></thead><tbody><tr><td>NMAE<sub>1</sub></td><td>0.2647</td><td>0.1516</td><td>0.1294</td><td>0.1294</td></tr><tr><td>NMAE<sub>2</sub></td><td>0.3974</td><td>0.2421</td><td>0.1966</td><td>0.1966</td></tr><tr><td>NMAE<sub>2</sub></td><td>0.2770</td><td>0.1400</td><td>0.1206</td><td>0.1206</td></tr><tr><td>sd(|c<sub>1</sub> - ĉ<sub>1</sub>|)</td><td>0.3735</td><td>0.1319</td><td>0.1106</td><td>0.1106</td></tr><tr><td>sd(|c<sub>2</sub> - ĉ<sub>2</sub>|)</td><td>0.4865</td><td>0.1940</td><td>0.1714</td><td>0.1714</td></tr><tr><td>sd(|c<sub>3</sub> - ĉ<sub>3</sub>|)</td><td>0.4368</td><td>0.1485</td><td>0.1534</td><td>0.1534</td></tr><tr><td>mean length confidence sets c<sub>1</sub></td><td>0.999</td><td>0.6847</td><td>0.4699</td><td>0.4956</td></tr><tr><td>mean length confidence sets c<sub>2</sub></td><td>1.447</td><td>0.8447</td><td>0.7671</td><td>1.0184</td></tr><tr><td>mean length confidence sets c<sub>3</sub></td><td>1.012</td><td>0.6218</td><td>0.3807</td><td>0.4176</td></tr><tr><td>coverage confidence sets c<sub>1</sub></td><td>91.7%</td><td>91.7%</td><td>92.1%</td><td>91.0%</td></tr><tr><td>coverage confidence sets c<sub>2</sub></td><td>94.2%</td><td>93.3%</td><td>95.6%</td><td>96.3%</td></tr><tr><td>coverage confidence sets c<sub>3</sub></td><td>95.9%</td><td>97.2%</td><td>91.9%</td><td>90.1%</td></tr></tbody></table>

slightly better results than the epistemic variance for $c_1$ and $c_3$, and this impact is bigger for $c_2$. We also note that the lengths of the confidence sets are quite different depending on the regions of the normalized parameters, for all methods. To understand in more detail what's happening, we compute the mean lengths of confidence sets, as well as coverages, for different regions of the normalized parameter $c_3$. The results are summarized in Table 5, and similar results were obtained for the two other parameters $c_1$.
<table><caption>Table 1: Sizes of datasets used for the different methods: N<sub>train</sub> for the training set, N<sub>val</sub> for the validation set, N<sub>cal</sub> for the calibration set and N<sub>test</sub> for the test set, and α the tolerance threshold used for the ABC standard and ABC-CNN. Methods used on the examples: stand. denotes the standard vanilla ABC, SMC the ABC-SMC, CNN the ABC-CNN and D-conf. our ABCD-conformal method. The language column is to specify which language has been used to obtain the presented results, even if code has been developped in both R and Julia.</caption><thead><tr><th rowspan="2">Example</th><th rowspan="2">N<sub>train</sub></th><th rowspan="2">N<sub>val</sub></th><th rowspan="2">N<sub>cal</sub></th><th rowspan="2">N<sub>test</sub></th><th rowspan="2">α</th><th colspan="5">ABC methods</th><th rowspan="2">Language</th></tr><tr><th>stand</th><th>SMC</th><th>CNN</th><th>RF</th><th>D-conf</th></tr></thead><tbody><tr><td>Section 4.1<br>MA2</td><td>10<sup>5</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0.001</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>Julia</td></tr><tr><td>Section 4.2<br>Lotka-Volterra</td><td>10<sup>5</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0.001</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>Julia</td></tr><tr><td>Section 4.3<br>Lake model</td><td>3 × 10<sup>4</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0.001</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>Julia</td></tr><tr><td>Appendix D<br>Gaussian field</td><td>7 × 10<sup>3</sup></td><td>1900</td><td>10<sup>3</sup></td><td>10<sup>2</sup></td><td>0.01</td><td>✓</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>R</td></tr></tbody></table>

et al. [2019] or Akesson et al. [2022] for instance). It is defined for observations X<sub>1</sub>, ..., X<sub>p</sub>;

$$X_j = Z_j + \theta_1 Z_{j-1} + \theta_2 Z_{j-2}, j = 1, \dots, p,$$

where (Z<sub>j</sub>)<sub>-2<j&le;p</sub> is an i.i.d. sequence of standard Gaussian N(0, 1). This model is identifiable in the following triangular region D:

$$-2 < \theta_1 < 2, \quad \theta_1 + \theta_2 > -1, \quad \theta_1 - \theta_2 < 1;$$

See for instance Marin and Robert [2007][Chapter 5] or [Marin et al., 2012, Section 2]. We assume that the prior distribution on (θ<sub>1</sub>, θ<sub>2</sub>) is the uniform distribution over this triangular region. This model allows for exact calculation of the posterior (θ<sub>1</sub>, θ<sub>2</sub>) | X. The goal is to obtain posterior estimates of the bi-dimensional vector of parameters (θ<sub>1</sub>, θ<sub>2</sub>) from an observed dataset of length p, (X<sub>1</sub>, ..., X<sub>p</sub>). In the following, we will use p = 100.

### 4.1.2 Algorithm parametrization

Number of used samples are given in Table 1. For the standard ABC, the summary statistics chosen are the first two autocovariances:

$$\tau_1 = \frac{1}{p-1} \sum_{j=2}^{p} x_j x_{j-1} - \left( \frac{1}{p} \sum_{j=1}^{p} x_j \right)^2 \quad \text{and} \quad \tau_2 = \frac{1}{p-2} \sum_{j=3}^{p} x_j x_{j-2} - \left( \frac{1}{p} \sum_{j=1}^{p} x_j \right)^2.$$

To compare two samples, we then use the L<sup>2</sup> distance between the two associated vectors of summary statistics. For ABC-SMC, the same summary statistics and distance are used, the algorithm consists of 10 iterations, in each iteration 10<sup>4</sup> simulations are made and the best
and $c_2$.

We can understand that the performances of the different methods can vary a lot depending on the region of the parameter. The coverages on the whole domain are quite different from the coverages on some specific regions. For instance, the coverage of standard ABC (resp. ABCD conformal) on the whole domain is 95.9% (resp 91.9%), while it is only 53.3% for $c_3 > 3$ (resp 80.0%). Hence, even if in theory ABCD-conformal gives marginal coverage and standard ABC gives conditional coverages, in this example standard ABC also does not give similar coverages depending on the region.

<table><caption>Table 5: Lokta Volterra example: comparison of standard ABC, ABC-CNN and ABCD-conformal (using overall uncertainty as uncertainty heuristic for the conformal procedure), in different regions for the normalized parameter $c_3$. These indicators are computed on a test set of size $10^3$.</caption><thead><tr><th></th><th></th><th>Standard ABC</th><th>ABC-CNN</th><th>ABCD-conformal overall</th></tr></thead><tbody><tr><td rowspan="2">$c_3 < 1$</td><td>mean length confidence sets</td><td>0.7271</td><td>0.4970</td><td>0.2808</td></tr><tr><td>coverage</td><td>96.3%</td><td>96.7%</td><td>91.9%</td></tr><tr><td rowspan="2">$1 < c_3 < 3$</td><td>mean length confidence sets</td><td>2.5856</td><td>1.2039</td><td>0.8366</td></tr><tr><td>coverage</td><td>98.4%</td><td>100%</td><td>93.6%</td></tr><tr><td rowspan="2">$3 < c_3$</td><td>mean length confidence sets</td><td>4.2436</td><td>2.9231</td><td>2.3087</td></tr><tr><td>coverage</td><td>53.3%</td><td>100%</td><td>80.0%</td></tr></tbody></table>

To sum up, ABCD-Conformal is better than the other methods in this example, with better predictons and narrower confidences sets, for similar coverages. We noted that the coverages are quite different in different regions of the parameters, illustrating the marginality of the coverage, but this is the case for ALL methods.

# 5 Discussion

In this article we propose a new ABC method that combines several approaches: the ABC framework, Neural Networks with Monte Carlo Dropout and a conformal procedure. This method is free of any summary statistic, distance or tolerance threshold. It is able to deal with multidimensional parameters, and gives exact non-asymptotic confidences sets.

In practice, this method is computationally efficient, and obtains good results. We test the method on three examples and compare its performances with other approaches: standard ABC, ABC-RF and ABC-CNN. We observe that depending on the problem at hand, standard ABC or ABC-RF can have slightly better accuracy. However, we also see in the Lokta-Volterra example that ABCD-conformal can also outperforms them. Regarding ABC-CNN, on all examples it never outperfoms the ABCD-conformal, either because of bad coverages of confidence sets, or because of larger mean lengths of these intervals. In contrast with the other methods, ABCD-conformal is the only one with guaranteed coverages for confidence sets. A big advantage of ABCD-conformal in our practice, is that it always gives both a good estimation accuracy and a good marginal frequentist coverage, which is not always the case for the other methods. It is an alternative to other methods when there is no obvious summary statistic. Moreover, the computing time is comparable with ABC-CNN and ABC-RF. The choice of the summary statistics is replaced by the choice of a network architecture. This choice can be guided by common Deep Learning architectures (imaging, time series, ...) and the massive associated literature. A drawback of our proposed method is that the coverage of the confidence sets is valid marginally, and not conditionally to some values of the parameters of interest. Indeed, on the Lokta-Volterra example we have seen that the performances of the methods can vary a lot depending on the region of the parameter. However, note that this was the case for ALL methods studied. For a more detailed comparison between the four ABC methods used, the reader can refer to Appendix B.

The ABCD-conformal algorithm proposed is promising. However, several improvements and modifications could be considered mostly at the level of the uncertainty proxy outputted by the neural network and the conformal procedure.
$10^3$ are kept to be the population of this iteration. Hence a total of $10^5$ simulations is made for each test sample, to be comparable with standard ABC and ABC-CNN.

In ABC-CNN method, the distance function between a sample from the test set and a sample from the training set is the quadratic $L^2$ distance between the parameters predicted by the CNN for the test sample, and true parameters used to simulate the training sample.

Regarding ABC-CNN and ABCD-Conformal, architectures of neural networks are the same for the two methods: 3 convolutional 1D layers with 128 neurons and a kernel size of 2 followed by max-pooling for the first two layers and by a flatten layer for the third one. Then, 3 dense layers of 100 neurons. Each of the previous layers uses the tanh activation function. The last layer to obtain the output of dimension 2 is linear. The raw samples of length $p = 100$ are the inputs of the neural networks, and the outputs are the bi-dimensional associated vectors of parameters $(\theta_1, \theta_2)$.

### 4.1.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 2. To compare ABCD-Conformal with the other methods, we computed confidence intervals for each component of $\theta$. Figure 1 shows the predicted parameters $\hat{\theta}_1$ and $\hat{\theta}_2$ estimated by the ABCD-Conformal method against the true values for the 1000 test samples, and confidence intervals in gray. For illustration, see Figure A.1 in Appendix A, which compares confidence ellipses on a single test sample. On the 1000 test samples the multidimensional coverage obtained using ABCD-conformal is of 94.6%.

<table><caption>Table 2: Moving Average 2 example: for each method and each component of $\theta$, we give the NMAE, $sd(|\theta - \hat{\theta}|)$, the mean and median lengths of the confidence or credible intervals and the coverage. For each method the mean and median areas of multidimensional confidence or credible ellipses are also given, as well as the associated coverage. These indicators are computed on a test set of size $10^3$. The best values for each line are highlighted in bold.</caption><thead><tr><th></th><th>Standard ABC</th><th>ABC-SMC</th><th>ABC-CNN</th><th>ABC-Conformal</th></tr></thead><tbody><tr><td>NMAE(θ<sub>1</sub>)</td><td>0.174</td><td>0.176</td><td>0.177</td><td><strong>0.166</strong></td></tr><tr><td>sd(|θ<sub>1</sub> - θ̂<sub>1</sub>|)</td><td>0.087</td><td>0.089</td><td>0.091</td><td><strong>0.083</strong></td></tr><tr><td>mean length CI(θ<sub>1</sub>)</td><td>0.544</td><td>0.626</td><td><strong>0.377</strong></td><td>0.560</td></tr><tr><td>median length CI(θ<sub>1</sub>)</td><td>0.535</td><td>0.609</td><td><strong>0.375</strong></td><td>0.559</td></tr><tr><td>coverage CI(θ<sub>1</sub>)</td><td>94.0%</td><td><strong>96.5%</strong></td><td>81.2%</td><td>95.7%</td></tr><tr><td>NMAE(θ<sub>2</sub>)</td><td>0.273</td><td>0.272</td><td>0.264</td><td><strong>0.234</strong></td></tr><tr><td>sd(|θ<sub>2</sub> - θ̂<sub>2</sub>|)</td><td>0.101</td><td>0.101</td><td>0.099</td><td><strong>0.089</strong></td></tr><tr><td>mean length CI(θ<sub>2</sub>)</td><td>0.615</td><td>0.667</td><td><strong>0.391</strong></td><td>0.583</td></tr><tr><td>median length CI(θ<sub>2</sub>)</td><td>0.596</td><td>0.653</td><td><strong>0.394</strong></td><td>0.585</td></tr><tr><td>coverage CI(θ<sub>2</sub>)</td><td>93.2%</td><td>95.8%</td><td>78.1%</td><td><strong>94.8%</strong></td></tr><tr><td>mean area CE(θ) (2D)</td><td>0.437</td><td>0.507</td><td><strong>0.187</strong></td><td>0.409</td></tr><tr><td>median area CE(θ) (2D)</td><td>0.412</td><td>0.479</td><td><strong>0.187</strong></td><td>0.411</td></tr><tr><td>coverage CE(θ) (2D)</td><td><strong>95.1%</strong></td><td>96.7%</td><td>75.6%</td><td>94.5%</td></tr></tbody></table>
Concerning the neural network, we can imagine to use the Dropconnect technique instead of Dropout. These two techniques prevent "co-adaptation" of units in a neural network. Dropout randomly drops hidden nodes, and Dropconnect drops connections (but all nodes can remain partially active). Dropconnect is a generalization of Dropout since there are more possible connections than nodes in a neural network. Another possibility to be explored, could be to use ensembles of neural networks (or deep ensembles) to obtain a random estimation of the parameter of interest and associated uncertainties, see Lakshminarayanan et al. [2017]. As explained by Srivastava et al. [2014c], Dropout can even be interpreted as ensemble model combination.

Regarding the conformal procedure, in this article we focused on conformalizing a scalar uncertainty estimate, because the parameter of interest $\theta$ is a vector of scalars. But we can also use conformalized quantile regression (see Angelopoulos and Bates [2023] for a presentation of this procedure). Finally, in the method proposed, we have considered a *split conformal procedure*. This is computationally attractive, as the model needs to be fitted only one time. But it requires to have a calibration set, in addition to the training and validation ones (even if in general this set is quite small compared to the training set). *Full conformal procedure* could avoid these extra simulations, at the cost of many more model fits, see Angelopoulos and Bates [2023]. Hence, choosing between split or full conformal procedure could depend on the problem at hand.

## References

Akesson, M., Singh, P., Wrede, F., and Hellander, A. (2022). Convolutional neural networks as summary statistics for approximate bayesian computation. *IEEE/ACM Transactions on Computational Biology and Bioinformatics*, 19(06):1–1.

Angelopoulos, A. N. and Bates, S. (2023). Conformal prediction: A gentle introduction. *Foundations and Trends® in Machine Learning*, 16(4):494–591.

Baragatti, M. and Pudlo, P. (2014). An overview on approximate bayesian computation*. *ESAIM: Proc.*, 44:291–299.

Besag, J. (1974). Spatial interactions and the statistical analysis of lattice systems. *Journal of the Royal Statistical Society. Series B. Statistical methodology*, 148:1–36.

Blum, M., Nunes, M., Prangle, D., and Sisson, S. (2012). A comparative review of dimension reduction methods in approximate bayesian computation. *stat sci 28: 189-208*. *Statistical Science*, 28.

Boston, T., Dijk, A., Larraondo, P., and Thackway, R. (2022). Comparing cnns and random forests for landsat image segmentation trained on a large proxy land cover dataset. *Remote Sensing*, 14:3396.

Breiman, L. (2001). Random forests. *Machine Learning*, 45:5–32.

Cannoodt, R., Saelens, W., Deconinck, L., and Saeys, Y. (2021). Spearheading future omics analyses using dyngen, a multi-modal simulator of single cells. *Nature Communications*.

Cliff, A. D. and Ord, J. K. (1981). *Spatial processes: models and applications*. Pion Ltd, London, UK.

Courville, A., Goodfellow, I., and Bengio, Y. (2016). *Deep Learning*. Adaptive computation and machine learning series. MIT Press.

Dantas, C. F., Drumond, T. F., Marcos, D., and Ienco, D. (2023). Counterfactual explanations for remote sensing time series data: An application to land cover classification. In *Joint European Conference on Machine Learning and Knowledge Discovery in Databases*, pages 20–36. Springer.

Datta, G. S., Ghosh, M., Mukerjee, R., and Sweeting, T. J. (2000). Bayesian prediction with approximate frequentist validity. *The Annals of Statistics*, 28(5):1414 – 1426.
We can see in Table 2 that the four methods obtain quite similar results concerning NMAE and standard deviation of the absolute error, with a slight advantage for ABCD-conformal. Some differences can be noted for the coverages, and the mean lengths of the confidence intervals. The coverage is sharp for standard ABC, ABC-SMC and ABCD-Conformal. However, it is too small for ABC-CNN. This is because the approximated posteriors are too peaked around the mean values, leading confidence intervals with a small coverage (and hence useless). Concerning the mean and median CI lengths, ABCD-conformal and standard ABC give the best results which are quite close (discarding ABC-CNN which does not have good coverage). Concerning the confidence ellipses, ABCD-conformal gives the smallest one (excluding again ABC-CNN). Note that it was expected that standard ABC and ABC-SMC performed quite well, as we used good summary statistics. Moreover, as showed by Frazier et al. [2018], conditions for good asymptotic properties are verified for standard ABC on this example, which is quite rare.

To sum up, ABC-CNN does not provide reliable credible intervals, ABCD-Conformal performs slightly better than all other methods in terms of NMAE and coverage (especially for the estimation of $\theta_2$) including the standard ABC and ABC-SMC, which are here in favorable estimation conditions.

**Remark 3.** We also evaluated the performances of the different methods by comparing the estimated $\theta$ values with the exact posterior means of $\theta$ rather than with the true $\theta$ used to generate the samples. In this example the exact posterior means can indeed be computed. As expected, and as explained at the beginning of section 4, the exact posterior means are very close to the true values, see Figure A.2 in Appendix A. Consequently, the results and performances obtained with the different methods are similar to those presented in Table 2, and the same conclusions can be drawn. See supplementary materials for more details.

## 4.2 Lotka-Volterra model

### 4.2.1 The model

The Lotka-Volterra model describes the dynamics of biological systems in which two species interact, one as a predator and the other as a prey. Here we consider a stochastic Markov jump process version of this model with state $(X_1, X_2) \in \mathbb{Z}^2$ representing prey and predator population sizes, and three parameters $\theta_1, \theta_2$ and $\theta_3$. Three transitions are possible, with hazard rates $\theta_1 X_1, \theta_2 X_1 X_2$ and $\theta_3 X_2$ respectively:

$$ (X_1, X_2) \xrightarrow{\theta_1, X_1} (X_1 + 1, X_2) \quad (\text{prey growth}) $$

$$ (X_1, X_2) \xrightarrow{\theta_2, X_1 X_2} (X_1 - 1, X_2 + 1) \quad (\text{predation interaction}) $$

$$ (X_1, X_2) \xrightarrow{\theta_3, X_2} (X_1, X_2 - 1) \quad (\text{predator death}) $$

The estimation of the parameters of this Lotka-Volterra model has been studied by several authors in an ABC framework, see for instance Prangle [2017]. The initial conditions are
de Haan, L. (1984). A spectral representation for max-stable processes. *The Annals of Probability*, 12(4):1194-1204.

Dereich, S., Scheutzow, M., and Schottstedt, R. (2013). Constructive quantization: Approximation by empirical measures. In *Annales de l'IHP Probabilités et statistiques*, volume 49, pages 1183-1203.

Fearnhead, P. and Prangle, D. (2012). Constructing summary statistics for approximate bayesian computation: semi-automatic approximate bayesian computation. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 74(3):419-474.

Filos, A., Tigkas, P., McAllister, R., Rhinehart, N., Levine, S., and Gal, Y. (2020). Can autonomous vehicles identify, recover from, and adapt to distribution shifts? In *International Conference on Machine Learning*, pages 3145-3153. PMLR.

Folgoc, L. L., Baltatzis, V., Desai, S., Devaraj, A., Ellis, S., Manzanera, O. E. M., Nair, A., Qiu, H., Schnabel, J., and Glocker, B. (2021). Is mc dropout bayesian? *arXiv preprint arXiv:2110.04286*.

Fournier, N. and Guillin, A. (2015). On the rate of convergence in wasserstein distance of the empirical measure. *Probability theory and related fields*, 162(3):707-738.

Frazier, D. T., Martin, G. M., Robert, C. P., and Rousseau, J. (2018). Asymptotic properties of approximate Bayesian computation. *Biometrika*, 105(3):593-607.

Gal, Y. (2016). *Uncertainty in deep learning*. PhD thesis, University of Cambridge.

Gal, Y. and Ghahramani, Z. (2016). Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In *international conference on machine learning*, pages 1050-1059.

Gal, Y., Hron, J., and Kendall, A. (2017). Concrete dropout. *31st Conference on Neural Information Processing System*.

Gobet, E., Lerasle, M., and Métivier, D. (2022). Mean estimation for randomized quasi monte carlo method. *Hal preprint hal-03631879v2*.

Grelaud, A., Marin, J.-M., Robert, C. P., Rodolphe, F., and Taly, J.-F. (2009). ABC likelihood-free methods for model choice in gibbs random fields. *Bayesian Analysis*, 4.

Guéron, A. (2019). *Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow, 2nd Edition*. O'Reilly Media, Inc.

Hinton, G. E., Srivastava, N., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2012). Improving neural networks by preventing co-adaptation of feature detectors. *CoRR*, abs/1207.0580.

Hoff, P. (2023). Bayes-optimal prediction with frequentist coverage control. *Bernoulli*, 29(2):901–928.

Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G. G. (2021). What are bayesian neural network posteriors really like? In *International conference on machine learning*, pages 4629-4640. PMLR.

Jiang, B., Wu, T.-Y., Zheng, C., and Wong, W. H. (2017). Learning summary statistic for approximate bayesian computation via deep neural network. *Statistica Sinica*, 27:1595-1618.

Jospin, L. V., Laga, H., Boussaid, F., Buntine, W., and Bennamoun, M. (2022). Hands-On Bayesian Neural Networks—A Tutorial for Deep Learning Users. *IEEE Computational Intelligence Magazine*, 17(2):29-48.

Joyce, P. and Marjoram, P. (2008). Approximately sufficient statistics and bayesian computation. *Statistical applications in genetics and molecular biology*, 7:Article26.

Kompa, B., Snoek, J., and Beam, A. L. (2021). Empirical frequentist coverage of deep learning uncertainty quantification procedures. *Entropy*, 23(12).
<figure><img src="image_1.png" /><figcaption>Figure 1: Moving Average 2 example: estimated parameters $\hat{\theta}_1$ and $\hat{\theta}_2$ using the ABCD-Conformal
against the true values for the 1000 test samples, and confidence sets obtained by the conformal
procedure (in gray intervals, for the epistemic uncertainty).</figcaption></figure>
Lakshminarayanan, B., Pritzel, A., and Blundell, C. (2017). Simple and scalable predictive uncertainty estimation using deep ensembles. In *Proceedings of the 31st International Conference on Neural Information Processing Systems*, NIPS’17, page 6405–6416, Red Hook, NY, USA. Curran Associates Inc.

LeCun, Y., Haffner, P., Bottou, L., and Bengio, Y. (1999). Object Recognition with Gradient-Based Learning, pages 319–345. Springer Berlin Heidelberg, Berlin, Heidelberg.

MacKay, D. J. (1992). A practical bayesian framework for backpropagation networks. *Neural computation*, 4(3):448–472.

Marin, J.-M., Pudlo, P., Robert, C. P., and Ryder, R. J. (2012). Approximate bayesian computational methods. *Statistics and computing*, 22(6):1167–1180.

Marin, J.-M. and Robert, C. (2007). *Bayesian Core: A Practical Approach to Computational Bayesian Statistics*. Springer New York, NY.

Marjoram, P., Molitor, J., Plagnol, V., and Tavaré, S. (2003). Markov chain Monte Carlo without likelihoods. *Proceedings of the National Academy of Sciences of the USA*, 100(26):15324–15328.

McKinley, T., Cook, A. R., and Deardon, R. (2009). Inference in epidemic models without likelihoods. *The International Journal of Biostatistics*, 5(1).

Meinshausen, N. (2006). Quantile regression forests. *Journal of Machine Learning Research*, 7(35):983–999.

Messoudi, S., Destercke, S., and Rousseau, S. (2022). Ellipsoidal conformal inference for multi-target regression. In *Conformal and Probabilistic Prediction with Applications*, pages 294–306. PMLR.

Neal, R. M. (2012). *Bayesian learning for neural networks*, volume 118. Springer Science & Business Media.

Piccioni, F., Casenave, C., Baragatti, M., Cloez, B., and Vinçon-Leite, B. (2022). Calibration of a complex hydro-ecological model through approximate bayesian computation and random forest combined with sensitivity analysis. *Ecological Informatics*, 71:101764.

Prangle, D. (2017). Adapting the ABC Distance Function. *Bayesian Analysis*, 12(1):289 – 309.

Pritchard, J., Seielstad, M., Perez-Lezaun, A., and Feldman, M. (1999). Population growth of human Y chromosomes: a study of Y chromosome microsatellites. *Molecular Biology and Evolution*, 16:1791–1798.

Raynal, L., Marin, J.-M., Pudlo, P., Ribatet, M., Robert, C. P., and Estoup, A. (2018). ABC random forests for Bayesian parameter inference. *Bioinformatics*, 35(10):1720–1728.

Rousseau, J. and Szabo, B. (2016). Asymptotic frequentist coverage properties of bayesian credible sets for sieve priors in general settings. *Annals of Statistics*, 48.

Simola, U., Cisewski-Kehe, J., and Wolpert, R. L. (2021). Approximate bayesian computation for finite mixture models. *Journal of Statistical Computation and Simulation*, 91(6):1155–1174.

Sisson, S. A., Fan, Y., and Tanaka, M. M. (2007). Sequential monte carlo without likelihoods. *Proceedings of the National Academy of Sciences*, 104(6):1760–1765.

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014a). Dropout: a simple way to prevent neural networks from overfitting. *The journal of machine learning research*, 15(1):1929–1958.

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014b). Dropout: A simple way to prevent neural networks from overfitting. *Journal of Machine Learning Research*, 15(56):1929–1958.

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014c). Dropout: a simple way to prevent neural networks from overfitting. *J. Mach. Learn. Res.*, 15(1):1929–1958.
$X_1(0) = 50$ and $X_2(0) = 100$ and a dataset corresponds to observations of state $(X_1, X_2)$ at times 0, 2, 4, ..., 36. As usual, all simulations with an extinction of either the preys or the predators were discarded: we are interested in the conditional law of survival. The prior distributions are independent uniforms $U[-6, 2]$ for the transformed parameters $\log(\theta_1), \log(\theta_2)$ and $\log(\theta_3)$. We are interested here in the estimation of $\theta = (\theta_1, \theta_2, \theta_3)$.

### 4.2.2 Algorithm parametrization

Number of used samples are given in Table 1. The Gillespie's stochastic simulation algorithm is used to generate all the samples. To improve performances of the algorithm, once simulations have been done, we used standardized versions of the parameters. The goal is then to obtain posterior estimates of the tri-dimensional vector of normalized parameters $(\theta_1, \theta_2, \theta_3)$, from an observed sample consisting of two time series of size 19.

Concerning the standard ABC, there is no obvious summary statistic. Instead, the distance function $d((x_1, x_2), (x_{1d}, x_{2d}))$ between a sample from the test set $\{(x_{1d}[i], x_{2d}[i])\}, i = 1, \dots, 19\}$ and a sample from the training set $\{(x_1[i], x_2[i]), i = 1, \dots, 19\}$, is given by the sum of squared differences:

$$d((x_1, x_2), (x_{1d}, x_{2d})) = \sum_{i=1}^{19} ((x_1[i] - x_{1d}[i])^2 + (x_2[i] - x_{2d}[i])^2)$$

For ABC-SMC, the same distance is used, the algorithm consists of 10 iterations, in each iteration 100 × 100 simulations are made and the best 100 are kept to be the population of this iteration. Hence a total of $10^5$ simulations is made for each test sample, to be comparable with standard ABC and ABC-CNN. See Appendix E.1 for the implementation details. Concerning ABC-CNN, the distance function between a test sample and a training sample is the quadratic distance between the parameters predicted by the CNN for the test sample, and the true parameters used to simulate the training sample. The architectures used for the neural networks are the same for the ABC-CNN and the ABCD-Conformal: 3 convolutional 1D layers with 128 neurons and a kernel size of 2, followed by max-pooling for the first two layers, and by a flatten layer for the third one. Then, 3 dense layers of 128 neurons. Each of the previous layers uses the tanh activation function. The last layer to obtain the output of dimension 3 is linear. The raw samples consisting of two time series of lengths 19 are the inputs of the neural networks, and the outputs are the tri-dimensional associated vectors of parameters $(\theta_1, \theta_2, \theta_3)$. Concerning ABCD-conformal, the default uncertainty measure is the overall uncertainty returned by the Dropout procedure. To show the benefits of using a "relevant" measure of uncertainty, we also present results using a constant uncertainty (only for confidence ellipsoids).

Note that this example is very challenging as $\theta$ is conditioned by the non-extinction condition (some prior values of $\theta$ can lead almost systematically to extinction). Moreover, $\theta_1$ and $\theta_3$ both have heavy tailed distributions. It means that in some rare cases their values can be very large. In general Machine Learning algorithms are not very robust to these
Vovk, V., Gammerman, A., and Saunders, C. (1999). Machine-learning applications of algorithmic randomness. In *Sixteenth International Conference on Machine Learning* (ICML-1999) (01/01/99), pages 444–453.

Wasserman, L. (2011). Frasian inference. *Statistical Science*, 26(3):322–325.

Weed, J. and Bach, F. (2019). Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance. *Bernoulli*, 25(4A):2620–2648.

Wiqvist, S., Mattei, P.-A., Picchini, U., and Frellsen, J. (2019). Partially exchangeable networks and architectures for learning summary statistics in approximate Bayesian computation. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 6798–6807. PMLR.

Wood, A. T. A. and Chan, G. (1994). Simulation of stationary gaussian processes in $[0, 1]$ d. *Journal of Computational and Graphical Statistics*, 3:409–432.

Xu, F., Uszkoreit, H., Du, Y., Fan, W., Zhao, D., and Zhu, J. (2019). Explainable AI: A brief survey on history, research areas, approaches and challenges. In Tang, J., Kan, M.-Y., Zhao, D., Li, S., and Zan, H., editors, *Natural Language Processing and Chinese Computing*, pages 563–574, Cham. Springer International Publishing.

# A Pseudo-codes for ABC-CNN and ABC-RF and some comments

## A.1 ABC-CNN

The approach of Akesson et al. [2022] is simply a small modification of the standard ABC, where CNNs are used to estimate the posterior mean of a parameter, which is then used as a summary statistic. The line 8 of the standard ABC algorithm (Algorithm 1) to compute summary statistics, is replaced by the prediction of a CNN trained on the reference table, see lines 11, 14 and 15 of Algorithm 3. The output of the algorithm is similar to the one of the standard ABC: for a given data sample $\mathbf{x}$ and a tolerance threshold, we hope that the parameters $\theta$ accepted during the process are approximately distributed from the posterior distribution $\pi(\theta | \mathbf{x})$. Hence to have asymptotic consistency and valid frequentist coverages, the same kind of conditions are needed than for Algorithm 1. But it's harder to check these conditions for ABC-CNN than for standard ABC, because the summary statistic used in the ABC-CNN is the output of a CNN which is a "black-box". It is therefore practically impossible to verify anything about this summary statistic.

Although, there are some important differences that are worth highlighting between the two methods. The main and most important one is that it is not necessary to give summary statistics, thanks to the use of a CNN which automatically generates a relevant one. However, a distance $d$ for comparing two data samples and a tolerance threshold $\alpha$ are still required. In practice, this distance is often simpler to determine in the case of ABC-CNN than in the case of standard ABC, because comparing two parameters is often simpler than comparing two sets of summary statistics, which can be numerous and very diverse. Theoretically, as explained by Frazier et al. [2018], for a fixed choice of summaries, the two-stage procedure advocated by Fearnhead and Prangle [2012] and used by Akesson et al. [2022], will not reduce the asymptotic variance over a point estimate of a parameter produced via Algorithm 1. However, this two-stage procedure does reduce the Monte Carlo error inherent in estimating the approximate posterior distribution $\pi_{\alpha,d}(\cdot | \eta(\mathbf{x}))$ by reducing the dimension of the statistics on which the matching in approximate Bayesian computation is based. Using this approach, we can hope to have a smaller global error of approximation of the true posterior, compared to standard ABC.

Moreover, a CNN can deal with different types of data, in particular high dimensional and complex data, like temporal or spatial data. Parameters can be easily estimated from these complex data by a CNN, while it can be difficult to find appropriate summary statistics to be used in a classical ABC.
extremes [Drenkow et al., 2022] which can lead to difficulties in the training. In Section 4.2.3, we hence study the results conditional on the value of $\theta$.

### 4.2.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 3. Looking

<table><caption>Table 3: Lotka-Volterra example: for each method and each component of $\theta$, we give the NMAE, $sd(|c - \hat{c}|)$, the mean and median lengths of the confidence or credible intervals and the coverage. For each method the mean and median volumes of multidimensional confidence or credible ellipses are also given, as well as the associated coverage. For the ABCD-conformal method, the overall uncertainty is used by default. A constant uncertainty has also been used for confidence ellipses. All the indicators are computed on a test set of size $10^3$. The best values for each line are highlighted in bold.</caption><thead><tr><th></th><th>Stand. ABC</th><th>ABC-SMC</th><th>ABC-CNN</th><th>ABCD-Conf. overall</th><th>ABCD-conf. constant</th></tr></thead><tbody><tr><td colspan="6" style="text-align:center; font-weight:bold;">θ<sub>1</sub> parameter</td></tr><tr><td>NMAE(θ<sub>1</sub>)</td><td>0.186</td><td>0.161</td><td>0.090</td><td><strong>0.082</strong></td><td></td></tr><tr><td>sd(|θ<sub>1</sub> - θ̂<sub>1</sub>|)</td><td>0.503</td><td>0.527</td><td>0.276</td><td><strong>0.260</strong></td><td></td></tr><tr><td>mean length CI(θ<sub>1</sub>)</td><td>0.764</td><td>0.396</td><td><strong>0.326</strong></td><td>0.381</td><td></td></tr><tr><td>median length CI(θ<sub>1</sub>)</td><td>0.316</td><td>0.208</td><td>0.231</td><td><strong>0.178</strong></td><td></td></tr><tr><td>coverage CI(θ<sub>1</sub>)</td><td><strong>96.7%</strong></td><td>93%</td><td>96.1%</td><td>94.5%</td><td></td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">θ<sub>2</sub> parameter</td></tr><tr><td>NMAE(θ<sub>2</sub>)</td><td>0.1882</td><td>0.1529</td><td>0.1543</td><td><strong>0.1403</strong></td><td></td></tr><tr><td>sd(|θ<sub>2</sub> - θ̂<sub>2</sub>|)</td><td>0.0029</td><td>0.0029</td><td><strong>0.0018</strong></td><td>0.0020</td><td></td></tr><tr><td>mean length CI(θ<sub>2</sub>)</td><td>0.0062</td><td>0.0036</td><td><strong>0.0028</strong></td><td>0.0050</td><td></td></tr><tr><td>median length CI(θ<sub>2</sub>)</td><td>0.0037</td><td><strong>0.0024</strong></td><td><strong>0.0024</strong></td><td><strong>0.0024</strong></td><td></td></tr><tr><td>coverage CI(θ<sub>2</sub>)</td><td><strong>95.6%</strong></td><td>94.2%</td><td>80.7%</td><td>94.9%</td><td></td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">θ<sub>3</sub> parameter</td></tr><tr><td>NMAE(θ<sub>3</sub>)</td><td>0.195</td><td>0.173</td><td>0.091</td><td><strong>0.081</strong></td><td></td></tr><tr><td>sd(|θ<sub>3</sub> - θ̂<sub>3</sub>|)</td><td>0.438</td><td>0.461</td><td><strong>0.218</strong></td><td>0.222</td><td></td></tr><tr><td>mean length CI(θ<sub>3</sub>)</td><td>0.697</td><td>0.365</td><td><strong>0.278</strong></td><td>0.318</td><td></td></tr><tr><td>median length CI(θ<sub>3</sub>)</td><td>0.268</td><td>0.175</td><td>0.177</td><td><strong>0.113</strong></td><td></td></tr><tr><td>coverage CI(θ<sub>3</sub>)</td><td><strong>97.2%</strong></td><td>92.4%</td><td>95.7%</td><td>94.6%</td><td></td></tr><tr><td colspan="6" style="text-align:center; font-weight:bold;">θ parameter (3 dimensional)</td></tr><tr><td>mean volume CE(θ)</td><td>0.0197</td><td>0.00269</td><td><strong>0.00207</strong></td><td>0.0336</td><td>0.0230</td></tr><tr><td>median volume CE(θ)</td><td>0.00028</td><td><strong>0.00008</strong></td><td>0.0004</td><td>0.00015</td><td>0.0230</td></tr><tr><td>coverage CE(θ)</td><td>95.9%</td><td>91.5%</td><td>88.8%</td><td>94.6%</td><td><strong>96.7%</strong></td></tr></tbody></table>

marginally at each component of $\theta$, we see that overall ABC-SMC is better than standard ABC, on all criteria except coverage, for which the standard ABC is slightly better. We also note that ABC-CNN and ABCD-Conformal outperform standard ABC and ABC-SMC and
**Algorithm 3:** Pseudo code for CNN as Summary Statistics for ABC [Akesson et al., 2022]

**Input** : A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, integers $N_{\text{train}}$ and $N_{\text{val}}$ representing sizes of training, and validation sets, a (pseudo)-distance $d$ between two samples, and a tolerance threshold $\alpha \in (0, 1]$.

**Output:** A set $\theta_1, \dots, \theta_{[N_{\text{train}}\alpha]}$ whose distribution is approximately $\pi(\theta | \mathbf{x})$

1 **Generation of a reference table (training dataset):**

2 **for** $j \leftarrow 1$ **to** $N_{\text{train}}$ **do**

3 &nbsp;&nbsp;&nbsp;Draw $\theta_j \sim \pi$;

4 &nbsp;&nbsp;&nbsp;Draw synthetic sample $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^{\top}$ from the model $f(\cdot | \theta_j)$;

5 **end**

6 **Generation of a validation dataset:**

7 **for** $j \leftarrow 1$ **to** $N_{\text{val}}$ **do**

8 &nbsp;&nbsp;&nbsp;Draw $\theta_j \sim \pi$;

9 &nbsp;&nbsp;&nbsp;Draw synthetic sample $\mathbf{x}_j = (x_{1,j}, \dots, x_{d,j})^{\top}$ from the model $f(\cdot | \theta_j)$;

10 **end**

11 **Train a CNN on the reference table:** inputs of the CNN are the {$\mathbf{x}_j, j = 1, \dots, N_{\text{train}}$}, and the outputs are the {$\theta_j, j = 1, \dots, N_{\text{train}}$}. The validation set is used to choose the architecture of the CNN.

12 **Summary statistics and distances:**

13 **for** $j \leftarrow 1$ **to** $N_{\text{train}}$ **do**

14 &nbsp;&nbsp;&nbsp;**$\mathbf{x}_j$ is given as input to the trained CNN $\eta_{\text{CNN}}$ to obtain an estimate of the posterior expected value for $\theta_j = \eta_{\text{CNN}}(\mathbf{x}_j)$;**;

15 &nbsp;&nbsp;&nbsp;Compute the distance $d_j = d(\eta_{\text{CNN}}(\mathbf{x}), \eta_{\text{CNN}}(\mathbf{x}_j))$;

16 **end**

17 Order these distances, i.e., $d_{(1)} < d_{(2)} < \dots < d_{(N_{\text{train}})}$;

18 Keep the $\theta_j$ corresponding to the $\lfloor N_{\text{train}}\alpha \rfloor$ smallest distances.
obtain quite better results than ABC-CNN, when looking at NMAE and standard deviations of the absolute errors. ABCD-Conformal shows similar results to ABC-CNN regarding these indicators, which shows the benefit of using neural networks and no summary statistics to obtain good predictions. Concerning the mean and median lengths of CI, sometimes it is ABC-CNN which is better, other times it is ABCD-conformal. Concerning the coverages of CI, when comparing ABC-CNN and ABCD-conformal we can note that ABC-CNN does not always have a good coverage, unlike ABCD-conformal. Hence, marginally and globally on all three dimensions, ABCD-conformal gives the most acceptable results. Looking at the muldimensional confidence and credible ellipses is quite interesting. Indeed, we can note that ABCD-Conformal give ellipses with the largest volume in mean, but not in median, and that mean and median volumes are very different. ABC-CNN and ABC-SMC are discarded because of their too small coverages.

In this example, to highlight the benefit of using a relevant heuristic uncertainty measure $\hat{V}(\mathbf{x})$ for the conformal procedure, we also show the results when $\hat{V}(\mathbf{x})$ is constant, see Table 3. We notice that in this case, the confidence sets also have a good coverage, thanks to the conformal procedure, but they do not adapt to the input as they are all of the same size (hence same mean and median). We therefore see the benefit of using a relevant uncertainty like the overall variance, which reflects the magnitude of model error: the volumes of the credible sets obtained are then small for most (easy) inputs, and large for the few hard ones. This is demonstrated by the median volume of the confidence sets being much smaller than the mean when using the ABC-Conformal method with non constant uncertainty $\hat{V}(\mathbf{x})$.

As noted before, the distributions of $\theta_1$ and $\theta_3$ have heavy tails. To show that this behaviour yields big differences between mean and median lengths of confidence and credible intervals and ellipses, we investigate the results for different regions of values of $\theta_3$ for illustration. The results are summarized in Table 4, and similar results were obtained for $\theta_1$. We can understand that the performances of the different methods can vary a lot depending on the region of the parameter. For a given component of $\theta$, the coverages on the whole domain are quite different from the coverages on some specific regions. For instance, the coverage of standard ABC (resp. ABCD-Conformal) on the whole domain is 97.2% (resp 94.6%), while it is only 77.5% (resp 87.5%) for extreme parameter region, $\theta_3 > 3$. Hence, even if in theory ABCD-Conformal gives marginal coverage and standard ABC gives conditional coverages, in this example standard ABC also does not give similar coverages depending on the region. Thanks to its good uncertainty heuristic, ABCD-Conformal is able to perform the least bad on harder parameters regions while being very good in easier regions $\theta < 1$ and $1 < \theta < 3$. In these easier regions, it has the best results in terms of size of the confidence sets, compared to standard ABC which is the only other method having a good coverage. Concerning ABC-SMC and ABC-CNN, they should be discarded most of the time because of too small coverages.

To sum up, ABCD-Conformal is most often better than the other methods in this example, with better predictions and smaller confidences intervals or ellipsoids, associated to good coverages. We note that the coverages are quite different in different regions of the parameters, illustrating the marginality of the coverage, but this is the case for all methods.
Another difference with the standard ABC is that a validation dataset should be generated, to be used to choose the network architecture, as usual for neural networks. This dataset is generally much smaller in size than the reference table used as a training set. Therefore, it generally adds little computation time compared with the effort required to generate the reference table. This lost computation time can be recovered, and time can even be saved, in the classical ABC step. Indeed, in this step, the sample of interest x for which we want to approximate the posterior π(θ | x) should be compared with all the samples in the reference table (lines 9 of Algorithm 1 and 15 of Algorithm 3). This comparison can be time consuming, depending on the distance used between summary statistics and the size of the reference table. If there are many summary statistics to be used by the standard ABC, the ABC-CNN which uses only one summary statistic, will be much faster. The larger the reference table, the greater the difference. This is why, in many cases, ABC-CNN will often be faster than standard ABC, despite the generation of a validation set. However, we can not generalise and it clearly depends on the problem at hand.

## A.2 ABC-RF

Raynal et al. [2018] who propose to conduct likelihood-free Bayesian inferences about parameters with no prior selection of the relevant components of summary statistics and bypassing the derivation of the associated tolerance threshold. The approach relies on the random forest (RF) methodology of Breiman [2001] applied in a regression setting.

<table><caption>**Algorithm 4:** Pseudo code for an Approximate Bayesian Computation - Random Forest sampler</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input</td><td>A Bayesian parametric model {f(· | θ), π}, a data sample x, N<sub>train</sub> the size of the training set, summary statistics η: ℝ<sup>d</sup> → ℝ<sup>l</sup>.</td></tr><tr><td>Output</td><td>For T(θ) a scalar transform of θ, approximations of the posterior expected value ψ<sub>η</sub>(x) = E<sub>π</sub>[T(θ) | η(x)], the posterior variance V<sub>π</sub>[T(θ) | η(x)], and posterior quantiles.</td></tr><tr><td>1</td><td>Generation of a reference table:</td></tr><tr><td>2</td><td>for j ← 1 to N<sub>train</sub> do</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Draw θ<sub>j</sub> ~ π;</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Draw synthetic sample x<sub>j</sub> = (x<sub>1,j</sub>, ..., x<sub>d,j</sub>)<sup>T</sup> from the model f(· | θ<sub>j</sub>)</td></tr><tr><td>5</td><td>end</td></tr><tr><td>6</td><td>Summary statistics:</td></tr><tr><td>7</td><td>for j ← 1 to N<sub>train</sub> do</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Compute the vector of summary statistics η(x<sub>j</sub>)</td></tr><tr><td>9</td><td>end</td></tr><tr><td>10</td><td>Train a random forest: inputs of the random forest are the {η(x<sub>j</sub>), j = 1, ..., N<sub>train</sub>}, and the outputs are the {T(θ<sub>j</sub>), j = 1, ..., N<sub>train</sub>).</td></tr><tr><td>11</td><td>Approximate ψ<sub>η</sub>(x): give η(x) as input of the trained random forest: the output of the RF allows to obtain easily approximations of ψ<sub>η</sub>(x) = E<sub>π</sub>[T(θ) | η(x)], V<sub>π</sub>[T(θ) | η(x)] and of posterior quantiles (hence credible intervals).</td></tr></tbody></table>

This algorithm automatises the inclusion of summary statistics and do not need the definition of a distance or the choice of a tolerance threshold. It appears mostly insensitive to the presence of non relevant summary statistics. It can then deal with a large number of summary statistics, by-passing any form of pre-selection. However, it is mandatory that relevant ones are present in the initial pool of summary statistics to be considered. To obtain posterior quantiles and then credible intervals, this ABC-RF method uses quantile regression forests to approximate the posterior cumulative distribution function of T(θ) given η(x), F(T(θ) | η(x)). The asymptotic consistency of this approach has been established by Meinshausen [2006], under some conditions. In particular, F(T(θ) | η(x)) should be Lipschitz continuous and strictly monotonously increasing in T(θ). Such conditions are practically unverifiable in an ABC framework where the likelihood is unknown or intractable.
For illustrative purposes, Figures B.1 and B.2 in Appendix B show estimated $\theta$ associated with 3D-confidence ellipsoids obtained using standard ABC, ABC-CNN and ABCD-conformal, and true values. Figure B.1 shows an example where ABCD-Conformal has the smallest confidence ellipsoid with the true $\theta$ inside. Figure B.2 shows an example where ABC-Conformal has the largest confidence region with the true $\theta$ inside, while other methods fail to capture the difficulty of estimating this parameter and provide too small confidence ellipsoids which do not include the true parameter.

<table><caption>Table 4: Lotka-Volterra example, results conditioned by the value of $\theta_3$: for each method and component $\theta_3$, we give the mean and median lengths of the confidence or credible intervals and the coverage. For each method the mean and median volumes of multidimensional confidence or credible ellipses are also given, as well as the associated coverage. For the ABCD-conformal method, the overall uncertainty is used by default. A constant uncertainty has also been used for confidence ellipses. All the indicators are computed on a test set of size $10^3$. The best values for each line are highlighted in bold.</caption><thead><tr><th></th><th>Stand. ABC</th><th>ABC-SMC</th><th>ABC-CNN</th><th>ABCD-Conf. overall</th><th>ABCD-Conf. constant</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">θ<sub>3</sub> &lt; 1, 782 test samples</td></tr><tr><td>mean length CI(θ<sub>3</sub>)</td><td>0.255</td><td>0.164</td><td>0.166</td><td><strong>0.120</strong></td><td></td></tr><tr><td>median length CI(θ<sub>3</sub>)</td><td>0.192</td><td>0.120</td><td>0.148</td><td><strong>0.0967</strong></td><td></td></tr><tr><td>coverage CI(θ<sub>3</sub>)</td><td>98.1%</td><td>94.4%</td><td><strong>97.1%</strong></td><td>94.8%</td><td></td></tr><tr><td colspan="6" style="text-align: center;">1 ≤ θ<sub>3</sub> ≤ 1, 178 test samples</td></tr><tr><td>mean length CI(θ<sub>3</sub>)</td><td>1.922</td><td>0.889</td><td>0.580</td><td><strong>0.560</strong></td><td></td></tr><tr><td>median length CI(θ<sub>3</sub>)</td><td>1.41</td><td>0.745</td><td>0.517</td><td><strong>0.399</strong></td><td></td></tr><tr><td>coverage CI(θ<sub>3</sub>)</td><td>97.8%</td><td>91.0%</td><td>94.4%</td><td>95.5%</td><td></td></tr><tr><td colspan="6" style="text-align: center;">3 &lt; θ<sub>3</sub>, 40 samples</td></tr><tr><td>mean length CI(θ<sub>3</sub>)</td><td>3.89</td><td>1.97</td><td><strong>1.13</strong></td><td>3.11</td><td></td></tr><tr><td>median length CI(θ<sub>3</sub>)</td><td>4.19</td><td>1.98</td><td><strong>1.13</strong></td><td>2.41</td><td></td></tr><tr><td>coverage CI(θ<sub>3</sub>)</td><td>77.5%</td><td>60.0%</td><td>75.0%</td><td><strong>87.5%</strong></td><td></td></tr><tr><td colspan="6" style="text-align: center;">θ parameter (3 dimensional), for θ<sub>3</sub> &lt; 1</td></tr><tr><td>mean volume CE(θ)</td><td>0.00124</td><td>0.00027</td><td><strong>0.00026</strong></td><td>0.00074</td><td>0.0230</td></tr><tr><td>median volume CE(θ)</td><td>0.00013</td><td><strong>0.00004</strong></td><td>0.00008</td><td>0.00007</td><td>0.0230</td></tr><tr><td>coverage CE(θ)</td><td>97.8%</td><td>93.6%</td><td>91.8%</td><td>95.4%</td><td>99.5%</td></tr><tr><td colspan="6" style="text-align: center;">θ parameter (3 dimensional), for 1 ≤ θ<sub>3</sub> ≤ 1</td></tr><tr><td>mean volume CE(θ)</td><td>0.0632</td><td>0.0077</td><td><strong>0.0049</strong></td><td>0.0227</td><td>0.0230</td></tr><tr><td>median volume CE(θ)</td><td>0.0341</td><td>0.0020</td><td><strong>0.0019</strong></td><td>0.0024</td><td>0.0230</td></tr><tr><td>coverage CE(θ)</td><td>94.4%</td><td>89.3%</td><td>80.3%</td><td><strong>94.9%</strong></td><td>94.4%</td></tr><tr><td colspan="6" style="text-align: center;">θ parameter (3 dimensional), for 3 &lt; θ<sub>3</sub></td></tr><tr><td>mean volume CE(θ)</td><td>0.188</td><td>0.028</td><td><strong>0.025</strong></td><td>0.724</td><td>0.0230</td></tr><tr><td>median volume CE(θ)</td><td>0.150</td><td>0.015</td><td><strong>0.019</strong></td><td>0.166</td><td>0.0230</td></tr><tr><td>coverage CE(θ)</td><td>65.0%</td><td>60.0%</td><td>67.5%</td><td>77.5%</td><td>52.5%</td></tr></tbody></table>
In practice, Raynal et al. [2018] showed on several examples that the approximations of posterior expectations obtained by ABC-RF were quite accurate, while posterior variances were slightly overestimated, and confidence sets slightly conservative (larger than the exact ones).

## B Detailed comparison of the four ABC methods used in the examples

Below the four methods used (standard ABC, ABC-CNN, ANC-RF and ABCD-conformal) are compared in details on different points, with a practical point of view. Table 4 summarizes these reflexions.

### Output of the algorithm

First, one big difference between these methods is what is obtained from the algorithm: for standard ABC and ABC-CNN an approximation of the whole posterior distribution of the parameter of interest is obtained. Concerning ABC-RF and ABCD-conformal, these methods focus on the approximation of transforms of interest of the posterior, like posterior mean, posterior variance or posterior quantiles for instance. It is rather common that only transforms of the posterior are of interest for practitioners, and in this case all four methods can be used. But if the whole posterior would be studied, ABC-RF and ABCD-conformal will not be adapted. Note that to obtain a good approximation of the whole posterior, in general more simulations are required than for simply estimating transforms of the posterior.

### Need of relevant summary statistic

Standard ABC need summary statistics to be able to compare simulations and observed data. Sometimes obvious relevant statistics exist (like in the MA2 example for instance), and the best is to have exhaustive statistics. But most of the time, it is not easy to find such relevant statistics and it is a serious bottleneck when performing inference on complex and high-dimensional problems. As explained in section 2.1, it is preferable to have a small number of summary statistics to avoid the burden of multidimensionality, and these statistics should be sufficiently informative [Fearnhead and Prangle, 2012]. Concerning ABC-RF, it enables to automatise the inclusion of summary statistics in an ABC algorithm. But it is necessary to have a set of statistics from which to select relevant ones: relevant statistics should be in the set, and this is a limitation. ABC-CNN and ANCD-conformal do not need any summary statistic, which is an advantage as the choice of these summaries have been proven to be crucial to obtaining good results.

### Need for a distance

A distance to compare simulations and observed data is needed for standard ABC. This distance depends on the summary statistics used and on the problem at hand. Specific distances can be used for genetics problem for instance. In general, a good knowledge of the problem in question and discussions with experts are necessary to define a relevant distance. In case of low-dimensional statistics, relevant distances are generally faster to compute. But in case of high-dimensional statistics, they can be quite long to compute. This is the case of the Lokta-Volterra example: we did not really used summary statistics as datasets and simulations were directly compared using an euclidian distance. But as the dataset was of dimension fourty (two series of twenty times), it was computationally demanding: the standard ABC was the method which took more computing time. Concerning the ABC-CNN algorithm, a summary statistic is obtained through the network, and the distance to be used is in general quite simple and fast to be computed. Of course computing time increases with the dimension of the parameters of interest, but in our examples with a maximum parameter dimension of three, it was very fast. ABC-RF and ABCD-conformal do not need any distance, which can avoid a tricky choice and can save computational time.

### Need for a tolerance threshold

For standard ABC and ABC-CNN, a tolerance threshold is also needed and its choice affects the degree of approximation obtained for the whole distribution. Usually, for practical reasons, quantile-based acceptance thresholds are used (quantiles of the distances between the observed data set and simulations in the training set). ABC-RF and ABCD-conformal do not need a tolerance threshold, as they do not use a distance.

### Dealing with multidimensional parameters

Concerning the parameter of interest, the ABC-RF method can estimate only uni-dimensional transforms
## 4.3 A model for lake ecosystem dynamics

### 4.3.1 The model

In this section, we address the problem of estimating nine parameters in the system of two differential equations introduced in Piccioni et al. [2022]. Although still a simplified model, it is more complex than those presented in the other sections. It involves a much larger number of parameters whose influence is highly non-linear. In addition, it is a non-dominated model; in particular, this means that it does not admit any likelihood with respect to a reference measure. Consequently, estimating the posterior law for such a model using MCMC methods is completely impossible and a likelihood-free approach, such as ABC method or its variants, is required (see for instance Macci and Spizzichino [2015]). It is worth noting that this difficulty is often overcome by adding random noise to the dynamics, whether or not this is justified for the problem under consideration.

The system under consideration represents a lake ecosystem in which the time evolution of phytoplankton and inorganic phosphorus concentrations is described by the following equations:

$$ \frac{dC_P}{dt} = \mu(C_{IP}, T)C_P - d_P C_P \quad (14) $$

$$ \frac{dC_{IP}}{dt} = V_{\min}(C_{OP}, C_{IP}, T) - a_{pc}\mu(C_{IP}, T)C_P + (1 - f_{OP})a_{pc}d_P C_P. \quad (15) $$

For the sake of presentation, we do not provide details here about the meaning of parameters; full information is given in the Appendix C. However, it should be noted that $C_{OP}$ is a linear combination of $C_P$ and $C_{IP}$, $T$ is a time-periodic function, $\mu$ and $V_{\min}$ are non-linear maps. In this model, there are 15 parameters and 2 initial conditions. We fixed these two initial conditions as well as 6 parameters that are not included in the estimation process. Our objective is to estimate the 9 remaining parameters corresponding to the vector $\theta$ defined as follows: $\theta = (\mu_{\max}, K_{IP}, d_P, f_{OP}, V_{\max}, K_{OP}, K_1, \theta_{\min}, \theta_{gr})$. A dataset corresponds to observations of state $(C_P, C_{IP})$ at times 0, 0.5, 1, ..., 364.5 days (i.e. on a one year time-scale, every 12 hours). The prior distributions are independent uniforms.

### 4.3.2 Algorithm parametrization

The Numbers of samples used are given in Table 1. In this model, the different parameters have different scales and units. Therefore, to train the Neural Networks, the training and validation sets were scaled. Concerning standard ABC, there is no obvious choice of summary statistic. Moreover, the trajectories of $C_P(t)$ and $C_{IP}(t)$ have very different scales. As a consequence, for the distance function between a test sample and a training sample, we chose to use the NMAE to get rid of the scale. For ABC-SMC, the same distance is used. The algorithm consists of 6 iterations, in each iteration $5 \times 10^3$ simulations are performed, and the best $10^2$ are kept to be the population of this iteration. Hence a total of $3 \times 10^4$ simulations is performed for each test sample, in order to ensure comparability with standard ABC and ABC-CNN. Concerning ABC-CNN, the distance function between a test sample
of interest, usually a projection on a given coordinate of the parameter. In the discussion of their article, Raynal et al. [2018] explained that their attempts to deal with multidimensional parameters were so far unfruitful. The standard ABC can deal with multidimensional parameter, but it suffers from the burden of multidimensionality, as summary statistics and distances are more difficult to define, and take more time to be computed in high dimension. ABC-CNN, despite having a standard ABC step in the algorithm, suffers less from this burden, as the summary statistic obtained through the CNN is of the same dimension than the parameter of interest (no more than dimension three in our exmaples). Concerning ABCD-conformal, multidimensionality is not a problem, as the computing time for the CNN with droupout and the conformal procedure does not really increases with the dimension (at least on the examples studied, with a maximum dimension of three).

## Choose network or random forest architecture

ABC-CNN and ABCD-conformal are based on neural networks. Hence the architecture of the networks to be used should be chosen: mainly the number of layers, the type of layer, the number of neurons in a layers, the size of the kernel for convolutional layers, and the activation functions. For the Dropout rate, we do not need to choose it when using the concrete Dropout approach of Gal et al. [2017]. In practice, for the examples studied, we did not spend too much time to choose the architecture, as it is quite fast to train different networks and to compare them on the validation set. The most influential parameters for us in our examples were the type of layers (convolutional or dense), and the activation functions. Similarly, to use ABC-RF some parameters should be chosen for the random forest, mainly the number of trees, the minimum node size and the proportion of summary statistics sampled at each split. For this method also, we did not encounter difficulties to set these parameters, standard values making the job on our example.

## Needed datasets

Standard ABC and ABC-RF need a training set of simulations, to be compared with the observed data. ABC-CNN also needs a validation set, which is used to choose the network architecture. Concerning the ABCD-conformal, it also needs a calibration set for the conformal procedure. This calibration set is usually of size 500 or 1000 depending on the accuracy we want for the marginal coverage [Angelopoulos and Bates, 2023]. Therefore, more simulations are needed for ABC-CNN than for standard ABC and ABC-RF, and even more for ABCD-conformal. However, the sizes of the validation and of the calibration sets are generally quite smaller than the size of the training set, and the increase in simulation time is small compared to what is required for the training set. This time can often be recovered in other steps of the algorithm.

## Theoretical justification of the methods

In this article, we are particularly interested by the estimation of parameters of interest, with associated confidence sets. The question of the frequentist coverage of the obtained intervals matters, as the intervals represent our uncertainty of our estimation. We then compared the frequentist coverage of the obtained intervals for the different methods, on three examples, and we studied the bibliography for the theoretical justification of all these methods.

As detailed in section 2.1, Frazier et al. [2018] studied asymptotic properties of standard ABC and gave conditions under which Bayesian confidence sets have valid frequentist coverage levels. Hence, in theory and asymptotically, under some regularity conditions, standard ABC and ABC-CNN methods enable posterior concentration and the *probability matching criterion*. But these conditions can not be verified in general, and this restricts the use and validity of these methods.

Concerning the ABC-RF, as said in section 2.2.2, it is based on quantile regression forests whose asymptotic consistency has been proved by Meinshausen [2006] under some conditions. But these conditions are practically unverifiable in an ABC framework.

On the opposite, the ABCD-conformal method guarantees a non-asymptotic marginal frequentist coverage, without any distributional or model assumption. The test points should solely come from the same distribution as the calibration points. Only the efficiency of the conformal procedure is in question, as the lengths of the confidence sets will depend on the quality of the uncertainty measure used in the procedure. This guaranteed coverage is a great advantage compared to the other methods. However, this coverage is guaranteed marginally, and not conditionally. Nevertheless, some metrics can be computed to check how close we are to conditional marginality, see Angelopoulos and Bates [2023].
and a training sample is defined as the quadratic distance between the parameters predicted by the CNN for the test sample, and the true parameters used to generate the training sample.

Again, the architectures used for the neural networks are similar to the previous applications: 3 convolutional 1D layers with 128 neurons and a kernel size of 2, followed by max-pooling for the first two layers, and by a flatten layer for the third one. Then, we have 3 dense layers of 128 neurons. Each of the previous layers uses the tanh activation function. The last layer which leads to the output of dimension 9 is linear. The raw samples consisting of two time series of lengths 730 (for $C_P$ and $C_{IP}$) are the inputs of the neural network, and the output is the multidimensional associated vector of 9 parameters.

### 4.3.3 Results

All indicators detailed in the beginning of Section 4 are presented in Table 5. In this Table 5, we see that ABCD-Conformal always has a good coverage, while this coverage is sometimes too small for the three other methods.

Looking marginally at each component of $\theta$, when considering NMAE and standard deviations of absolute errors, ABCD-Conformal often outperforms standard ABC, ABC-SMC and ABC-CNN, see for instance $\theta_1$ or $\theta_9$. Concerning the mean and median lengths of the confidence intervals, they are smaller for ABCD-Conformal compared to the other methods, except for $\theta_5$, $\theta_6$ and $\theta_8$. But for these parameters the other methods have too small coverages (except for $\theta_8$ for which ABC-SMC performs particularly well), while the coverage of ABCD-conformal remains good.

Looking at the performances for the 9-dimensional parameter $\theta$, ABCD-conformal gives confidence sets with a quite smaller volume than ABC-standard and ABC-CNN, and with a larger volume than ABC-SMC. However, it is the only method with a satisfactory coverage. Figures C.1, C.2, C.3 and C.4 in appendix C show scatterplots of estimated parameters against the true values for the 1000 test samples and associated confidence sets, for the four methods. This higher-dimensional example clearly illustrates the pitfalls of ABC methods based on summary statistics and threshold selection. Indeed, for standard ABC and ABC-SMC, comparing two large multidimensional time series, even with the informed choice of using the NMAE instead of the naive $L_2$ norm, is not sufficient to capture the effect of each parameter. Regarding ABC-CNN, which circumvents this issue by using a neural network as summary statistics, it still needs to select neighbors to obtain the posterior distribution. As a result, it either overestimates the size of its confidence intervals compared to ABCD-Conformal or fails to ensure coverage, as seen for $\theta_5$, $\theta_6$ and $\theta_7$. Both phenomena hold when considering the multidimensional confidence regions.

In such models, some parameters are often not identifiable. In this particular example, the parameters $V_{max}$, $K_{OP}$ and $K_I$ (i.e. $\theta_5$, $\theta_6$, $\theta_7$) appear in the same function $V_{min}$. Their respective impacts on the dynamics are therefore difficult to capture which explain why they are not well identified.

For models of this kind, it is more common to address the problem of calibration rather than parameter estimation. Calibration consists in finding a set of parameters for which
## Benefits of using concrete Dropout

Apart from the fact that epistemic and aleatoric uncertainties are used to obtain credibility intervals during the conformal procedure, their values themselves and the distinction between these uncertainties are interesting for effective exploration of the uncertainty: which part can be attributed to the model or to the noise in the data, in which circumstances?

Est ce que dans le tableau je rajoute un point "Simplicity of method coding"? Ou bien j'en parle juste dans le texte? Car la simplicité d'utilisation dépend de la complexité de la méthode en elle-même, mais aussi du fait de devoir trouver/calculer des summary stats pour les méthodes ABC standard et ABC-CNN. En gros ABC-CNN et ABCD-conformal peuvent paraitre plus complexes, mais en même temps elles ne nécessitent aucun traitement préalable des données. Donc dur de mettre un smiley généralisant dans le tableau.
the simulated outputs are close to the observed data. It would require a more detailed
analysis, such as that performed in Piccioni et al. [2022]). For illustration, Figure 2 shows
one example of the reconstructed time series of phytoplankton and inorganic phosphorus
concentrations, obtained using the estimated vector of parameters θ (estimations obtained
with ABC-CNN and ABCD-conformal). Note that accurately reconstructing the dynamics
requires simultaneously estimating all parameters correctly. Therefore, the multidimensional
coverage and the size of the confidence region are more meaningful than their marginal
counterparts.

<figure><img src="image_2.png" /><figcaption>Figure 2: Lake ecosystem dynamics example: reconstructions of time series of phytoplankton and inorganic phosphorus, using the estimated θ parameters (estimations obtained with ABC-CNN and ABCD-conformal), for the 10th test sample.</figcaption></figure>

# 5 Discussion and conclusion

In this article we propose a new ABC method that combines several approaches: the ABC framework, Neural Networks with Monte Carlo Dropout and a conformal procedure. This method is free of any summary statistic, distance, or tolerance threshold. Indeed, it does not require samples to be compared. It behaves well in multidimensional settings, and provides exact non-asymptotic confidences sets.

In practice, this method is computationally efficient, and obtains good results. We test
the method on four examples and compare its performances with other approaches: standard
ABC, ABC-SMC, ABC-RF and ABC-CNN. ABCD-conformal often outperforms the other
methods. Even though ABC-SMC, ABC-CNN or ABC-RF may sometimes achieve better
<table><caption>Table 6: Summary table comparing the four methods: standard ABC, ABC-CNN, ABC-RF and ABCD-conformal.</caption><thead><tr><th>Point of comparison</th><th>Standard ABC</th><th>ABC-CNN</th><th>ABC-RF</th><th>ABCD-conformal</th></tr></thead><tbody><tr><td>Approx. of the whole posterior or of transforms of interest</td><td>whole posterior</td><td>whole posterior</td><td>transformed interest</td><td>transformed interest</td></tr><tr><td>No need of relevant summary stat</td><td>✗</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>No need of a distance</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>No need of a tolerance threshold</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td></tr><tr><td>Deal with multidimensional parameters</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>No need to choose network or random forest architecture</td><td>✓</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>Datasets needed</td><td>1: training</td><td>2: training; validation</td><td>1: training</td><td>3: training; validation; calibration</td></tr><tr><td>Justification of the method, guarantees</td><td>asymptotic under conditions</td><td>asymptotic conditions too difficult to check</td><td>asymptotic conditions too difficult to check</td><td>non asymptotic no condition</td></tr><tr><td>Adapted for different types of data and high-dim data</td><td>difficult</td><td>✓</td><td>difficult</td><td>✓</td></tr><tr><td>Computing time</td><td>✗</td><td>difficult to compare, it depends on examples.</td><td>✗</td><td></td></tr></tbody></table>
NMAE results in certain dimensions, we find that the performance of ABCD-conformal is consistently good, among the best if not the best across all parameter dimensions. Regarding ABC-RF, it produced very good results in the unidimensional example, better than ABCD-conformal. However, this algorithm is not adapted for multidimensional parameter estimation. As for ABC-CNN, in low dimensions this method provides good NMAE and is relatively easier to train than ABCD-Conformal. However, it still requires comparing samples using a threshold, which is difficult to calibrate. Moreover, in our experiments, it most often failed to correctly capture the posterior distribution, i.e. its confidence intervals were either too narrow (poor coverage) or, in some cases, too wide, as observed in the lake example (section 4.3). In contrast with the other methods, ABCD-Conformal is the only one with guaranteed coverages for confidence sets. A big advantage of ABCD-Conformal in our practice, is that it always gives both a good estimation accuracy and a good marginal frequentist coverage, which is not always the case for the other methods. In particular in large dimensions, as observed in the lake example (section 4.3), it does not suffer as much as other methods of the curse of dimensionality, providing the best NMAE (within a few percent) and the smallest confidence regions with good coverage. It is an alternative to other methods when there is no obvious summary statistic. The choice of the summary statistics is replaced by the choice of a network architecture. This choice can be guided by common Deep Learning architectures (imaging, time series, ...) and the massive associated literature. We used an almost identical architectures in the four tested examples, i.e. some CNNs with MaxPool to analyse the input features and a few Dense layers to transform back to the parameter space. This is a very different strategy than methods with summary statistics that cannot be applied to different examples. Last but not least, the ABCD-conformal is an amortized procedure, in the sense that once the network has been trained on the reference table and the conformal quantile has been obtained from the calibration set, both this network and this quantile can be reused for any new data sample of interest. There is no need to retrain a network, to recalculate a quantile, or run additional simulations. This nice feature holds for any method that relies on a reference table (such as standard ABC, ABC-CNN or ABC-RF), but not for sequential methods like ABC-SMC.

Nevertheless, the ABCD-conformal approach has certain limitations. One drawback concerns multimodal posteriors. Indeed, the goal of the ABCD-conformal method is to provide point estimates of parameters of interest, along with associated confidence intervals. Unlike standard ABC, ABC-SMC or ABC-CNN, it does not produce an approximation of the full posterior distribution. As a result, in the case of a multimodal posterior, the method is unable to detect such multimodality. Another weakness is that the coverage of the confidence sets is valid marginally, and not conditionally to some values of the parameters of interest. Indeed, on the Lotka-Volterra example (Section 4.2.3) we have seen that the performances of the methods can vary a lot depending on the region of the parameter. However, it should be noted that this limitation applies to all the methods studied. This issue is nonetheless mitigated in the case of ABCD-conformal, thanks to the heuristic uncertainty, which allows the method to adapt its confidence sets according to whether the parameter values are easy or difficult to estimate. In practice, we can note that training neural networks to output both a
parameter value and an associated heuristic uncertainty is more challenging than predicting the parameter value alone. Alternative uncertainty quantification methods beyond using the heteroscedastic loss with MC Dropout could be considered (see below). For a more detailed comparison between the five ABC methods used in our examples, the reader can refer to Appendix F. Table 6 gives a summary of this comparison.

The ABCD-Conformal algorithm proposed is promising. However, several improvements and modifications could be considered mostly at the levels of the uncertainty proxy outputted by the neural network and of the conformal procedure. Concerning the neural network, we can imagine using the Dropconnect technique instead of Dropout. These two techniques prevent “co-adaptation” of units in a neural network. Dropout randomly drops hidden nodes, and Dropconnect drops connections (but all nodes can remain partially active). Dropconnect is a generalization of Dropout since there are more possible connections than nodes in a neural network. Another possibility to be explored, could be to use ensembles of neural networks (or deep ensembles) to obtain a random estimation of the parameter of interest and associated uncertainties, see Lakshminarayanan et al. [2017]. As explained by Srivas-tava et al. [2014c], Dropout can even be interpreted as ensemble model combination. More generally, any method for obtaining an uncertainty associated with the output of a neural network could be used. At the current rate of development of such networks, it seems to us that many methods could be used in the near future. Regarding the conformal procedure, in this article we focused on conformalizing a scalar uncertainty estimate, because the parameter of interest θ is a vector of scalars. But we can also use conformalized quantile regression (see Angelopoulos and Bates [2023] for a presentation of this procedure). Finally, in the method proposed, we have considered a *split conformal procedure*. This is computationally attractive, as the model needs to be fitted only one time. But it requires having a calibration set, in addition to the training and validation ones (even if in general this set is quite small compared to the training set). *A full conformal procedure* could avoid these extra simulations, at the cost of many more model fits, see Angelopoulos and Bates [2023]. Hence, choosing between split or full conformal procedure could depend on the problem at hand.

## Supplementary materials

**Appendix** It contains additional figures and technical details for the applications, an additional example, pseudocodes and comments for ABC-SMC, ABC-CNN and ABC-RF, and a detailed comparison of the five methods used in the applications. It can be read below the body of this manuscript.

**GitLab repository** The source codes, some datasets and pretrained neural networks to run the examples in section 4 are available in the GitLab repository https://forge.inrae.fr/mistea/codes.articles/abcdconformal with the associated README.md file. The results are presented as Quarto (or RMardown) notebooks with all instructions to ensure reproducibility.
**Website** Results and html-rendered notebooks can also be directly consulted on the website https://mistea.pages-forge.inrae.fr/codes.articles/abcdconformal/.

**Julia package** The Julia code of all methods used in this article is contained in the Julia package `ABCMethods.jl` accessible in the following Git repository https://github.com/dmetivie/ABCMethods.jl.

Most of the results presented in this article have been obtained using Julia. However, our ABCD-conformal approach was first developed using R. Hence, even if the R results are not presented, R code is provided in the GitLab repository and in the website.

## Acknowledgments

We thank Nicolas Verzelen for introducing us to the conformal approach during a seminar, and for the helpful discussions that followed. We also thank Mathieu Ribatet for the discussions that led to the initiation of this work.

For the applications, we mainly used the Julia programming language [Bezanson et al., 2017] and its scientific computing ecosystem. We acknowledge here the main open source packages used for this work. Deep learning components were developed using `Lux.jl` [Pal, 2025], with automatic differentiation via `Zygote.jl` [Innes, 2018] and GPU support through `CUDA.jl` [Besard et al., 2019]. Differential equations and jump processes were modeled using tools from the SciML ecosystem, including `DifferentialEquations.jl` [Rackauckas and Nie, 2017], `JumpProcesses.jl`, and `Catalysis.jl`. Data handling and visualization were supported by `DataFrames.jl` [Bouchet-Valat and Kamiński, 2023], `Distributions.jl` [Besançon et al., 2021], and `Plots.jl` [Christ et al., 2023].

Finally, we are grateful to the co-editor, the associate editor, and two anonymous reviewers for their insightful and constructive comments for improving the quality of the article.