<table><caption>Table 5: Lake system dynamics 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 volumes of multidimensional confidence or credible ellipses are also given, as well as the associated coverage. 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.</th></tr></thead><tbody><tr><td colspan="5"><strong>θ₁ parameter</strong></td></tr><tr><td>NMAE(θ₁)</td><td>0.6988</td><td>0.6298</td><td>0.3438</td><td><strong>0.2318</strong></td></tr><tr><td>sd(|θ₁ - ̂θ₁|)</td><td>0.4150</td><td>0.4523</td><td>0.2096</td><td><strong>0.1782</strong></td></tr><tr><td>mean length CI(θ₁)</td><td>2.5674</td><td>2.0100</td><td>1.6073</td><td><strong>1.1071</strong></td></tr><tr><td>median length CI(θ₁)</td><td>2.6670</td><td>2.3142</td><td>1.6073</td><td><strong>1.0822</strong></td></tr><tr><td>coverage CI(θ₁)</td><td>93.9%</td><td>91.6%</td><td>95.4%</td><td><strong>96.0%</strong></td></tr><tr><td colspan="5"><strong>θ₂ parameter</strong></td></tr><tr><td>NMAE(θ₂)</td><td>0.6469</td><td>0.6526</td><td>0.3540</td><td><strong>0.3028</strong></td></tr><tr><td>sd(|θ₂ - ̂θ₂|)</td><td>0.3992</td><td>0.3997</td><td><strong>0.2293</strong></td><td>0.2411</td></tr><tr><td>mean length CI(θ₂)</td><td>2.3185</td><td>2.4052</td><td>1.6189</td><td><strong>1.2964</strong></td></tr><tr><td>median length CI(θ₂)</td><td>2.4200</td><td>2.5227</td><td>1.6079</td><td><strong>1.1795</strong></td></tr><tr><td>coverage CI(θ₂)</td><td>92.6%</td><td>93.7%</td><td><strong>95.1%</strong></td><td>94.3%</td></tr><tr><td colspan="5"><strong>θ₃ parameter</strong></td></tr><tr><td>NMAE(θ₃)</td><td>0.5355</td><td>0.5295</td><td><strong>0.2891</strong></td><td>0.2957</td></tr><tr><td>sd(|θ₃ - ̂θ₃|)</td><td>0.3594</td><td>0.4037</td><td><strong>0.1968</strong></td><td>0.2401</td></tr><tr><td>mean length CI(θ₃)</td><td>1.9546</td><td>1.5897</td><td>1.5722</td><td><strong>1.2270</strong></td></tr><tr><td>median length CI(θ₃)</td><td>1.9927</td><td>1.5891</td><td>1.5667</td><td><strong>1.1233</strong></td></tr><tr><td>coverage CI(θ₃)</td><td>90.7%</td><td>86.3%</td><td><strong>96.3%</strong></td><td>94.7%</td></tr><tr><td colspan="5"><strong>θ₄ parameter</strong></td></tr><tr><td>NMAE(θ₄)</td><td>0.5812</td><td>0.5674</td><td>0.2928</td><td><strong>0.2889</strong></td></tr><tr><td>sd(|θ₄ - ̂θ₄|)</td><td>0.3770</td><td>0.3930</td><td><strong>0.2130</strong></td><td>0.2314</td></tr><tr><td>mean length CI(θ₄)</td><td>2.1368</td><td>1.7773</td><td>1.5487</td><td><strong>1.2091</strong></td></tr><tr><td>median length CI(θ₄)</td><td>2.2571</td><td>1.8893</td><td>1.5386</td><td><strong>1.0860</strong></td></tr><tr><td>coverage CI(θ₄)</td><td>92.8%</td><td>86.4%</td><td><strong>94.2%</strong></td><td>93.0%</td></tr><tr><td colspan="5"><strong>θ₅ parameter</strong></td></tr><tr><td>NMAE(θ₅)</td><td>0.8376</td><td>0.7146</td><td>0.6640</td><td><strong>0.5889</strong></td></tr><tr><td>sd(|θ₅ - ̂θ₅|)</td><td>0.4677</td><td>0.4640</td><td>0.4505</td><td><strong>0.3907</strong></td></tr><tr><td>mean length CI(θ₅)</td><td>2.8402</td><td>2.3203</td><td><strong>1.9357</strong></td><td>2.3954</td></tr><tr><td>median length CI(θ₅)</td><td>2.8960</td><td>2.4176</td><td><strong>1.9476</strong></td><td>2.4514</td></tr><tr><td>coverage CI(θ₅)</td><td>90.2%</td><td>88.4%</td><td>83.6%</td><td><strong>95.2%</strong></td></tr><tr><td colspan="5"><strong>θ₆ parameter</strong></td></tr><tr><td>NMAE(θ₆)</td><td>0.7596</td><td>0.6933</td><td>0.7396</td><td><strong>0.6555</strong></td></tr><tr><td>sd(|θ₆ - ̂θ₆|)</td><td>0.4392</td><td>0.4596</td><td>0.4870</td><td><strong>0.4180</strong></td></tr><tr><td>mean length CI(θ₆)</td><td>2.5613</td><td>2.1558</td><td><strong>1.9344</strong></td><td>2.5650</td></tr><tr><td>median length CI(θ₆)</td><td>2.6478</td><td>2.3235</td><td><strong>1.9508</strong></td><td>2.9140</td></tr><tr><td>coverage CI(θ₆)</td><td>90.0%</td><td>86.9%</td><td>78.6%</td><td><strong>95.4%</strong></td></tr><tr><td colspan="5"><strong>θ₇ parameter</strong></td></tr><tr><td>NMAE(θ₇)</td><td>0.8242</td><td>0.6419</td><td>0.4766</td><td><strong>0.4041</strong></td></tr><tr><td>sd(|θ₇ - ̂θ₇|)</td><td>0.4709</td><td>0.3925</td><td>0.3364</td><td><strong>0.3133</strong></td></tr><tr><td>mean length CI(θ₇)</td><td>2.7555</td><td>2.2553</td><td><strong>1.7518</strong></td><td>1.7600</td></tr><tr><td>median length CI(θ₇)</td><td>2.8220</td><td>2.4063</td><td>1.7186</td><td><strong>1.6961</strong></td></tr><tr><td>coverage CI(θ₇)</td><td>90.2%</td><td>88.1%</td><td>90.4%</td><td><strong>94.4%</strong></td></tr><tr><td colspan="5"><strong>θ₈ parameter</strong></td></tr><tr><td>NMAE(θ₈)</td><td>0.4623</td><td><strong>0.1661</strong></td><td>0.2506</td><td>0.1750</td></tr><tr><td>sd(|θ₈ - ̂θ₈|)</td><td>0.3139</td><td>0.1740</td><td>0.1749</td><td><strong>0.1333</strong></td></tr><tr><td>mean length CI(θ₈)</td><td>2.1980</td><td><strong>0.7461</strong></td><td>1.5547</td><td>0.8071</td></tr><tr><td>median length CI(θ₈)</td><td>2.2147</td><td><strong>0.6423</strong></td><td>1.5445</td><td>0.7914</td></tr><tr><td>coverage CI(θ₈)</td><td>94.6%</td><td>95.5%</td><td><strong>96.6%</strong></td><td>95.2%</td></tr><tr><td colspan="5"><strong>θ₉ parameter</strong></td></tr><tr><td>NMAE(θ₉)</td><td>0.3888</td><td>0.5478</td><td>0.2545</td><td><strong>0.1413</strong></td></tr><tr><td>sd(|θ₉ - ̂θ₉|)</td><td>0.2760</td><td>0.3773</td><td>0.1717</td><td><strong>0.1143</strong></td></tr><tr><td>mean length CI(θ₉)</td><td>1.6639</td><td>1.8246</td><td>1.5408</td><td><strong>0.6059</strong></td></tr><tr><td>median length CI(θ₉)</td><td>1.6985</td><td>1.8451</td><td>1.5252</td><td><strong>0.5815</strong></td></tr><tr><td>coverage CI(θ₉)</td><td>94.9%</td><td>89.9%</td><td><strong>96.9%</strong></td><td>95.0%</td></tr><tr><td colspan="5"><strong>θ parameter, 9-dimensional</strong></td></tr><tr><td>mean length CE(θ)</td><td>2235</td><td><strong>57.12</strong></td><td>823.3</td><td>110.5</td></tr><tr><td>median length CE(θ)</td><td>11683</td><td><strong>14.70</strong></td><td>765.0</td><td>102.8</td></tr><tr><td>coverage CE(θ)</td><td>90.3%</td><td>73.4%</td><td>87.9%</td><td><strong>95.4%</strong></td></tr></tbody></table>
<table><caption>Table 6: Summary table comparing the five methods: standard ABC, ABC-SMC, ABC-CNN, ABC-RF and ABCD-Conformal. A detailed comparison of these methods can be found in appendix F.</caption><thead><tr><th>Point of comparison</th><th>Standard ABC</th><th>ABC-SMC</th><th>ABC-CNN</th><th>ABC-RF</th><th>ABCD-Conf</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>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><td>✓</td></tr><tr><td>No need of a distance</td><td>✗</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><td>✓</td></tr><tr><td>Deal with multidimensional parameters</td><td>✓</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><td>✗</td></tr><tr><td>Datasets needed</td><td>1: training</td><td>Simulations made sequentially</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 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>difficult</td><td>✓</td><td>difficult</td><td>✓</td></tr><tr><td>Computing time</td><td></td><td></td><td>difficult to compare, it depends on examples.</td><td></td><td></td></tr></tbody></table>
The authors report there are no competing interests to declare.

## 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.

Beaumont, M., Cornuet, J., Marin, J. M., and Robert, C. (2009). Adaptive approximate Bayesian computation. *Biometrika*, 96(4):983–990.

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.

Besançon, M., Papamarkou, T., Anthoff, D., Arslan, A., Byrne, S., Lin, D., and Pearson, J. (2021). Distributions.jl: Definition and Modeling of Probability Distributions in the JuliaStats Ecosystem. *Journal of Statistical Software*, 98:1–30.

Besard, T., Foket, C., and De Sutter, B. (2019). Effective Extensible Programming: Unleashing Julia on GPUs. *IEEE Transactions on Parallel and Distributed Systems*, 30(4):827–841.

Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B. (2017). Julia: A Fresh Approach to Numerical Computing. *SIAM Review*, 59(1):65–98.

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.

Bouchet-Valat, M. and Kamiński, B. (2023). DataFrames.jl: Flexible and Fast Tabular Data in Julia. *Journal of Statistical Software*, 107:1–32.

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

Christ, S., Schwabeneder, D., Rackauckas, C., Borregaard, M. K., and Breloff, T. (2023). Plots.jl – A User Extendable Plotting API for the Julia Programming Language. *Journal of Open Research Software*, 11(1).
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.

de Haan, L. (1984). A spectral representation for max-stable processes. *The Annals of Probability*, 12(4):1194–1204.

Del Moral, P., Doucet, A., and Jasra, A. (2012). An adaptive sequential monte carlo method for approximate bayesian computation. *Statistics and Computing*, 22:1009–1020.

Dereich, S., Scheutzow, M., and Schottstedt, R. (2013). Constructive quantization: Approximation by empirical measures. 49(4):1183–1203.

Drenkow, N., Sani, N., Shpitser, I., and Unberath, M. (2022). A Systematic Review of Robustness in Deep Learning for Computer Vision: Mind the gap?

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.

Filippi, S., Barnes, C. P., Cornebise, J., and Stumpf, M. P. (2013). On optimality of kernels for approximate bayesian computation using sequential monte carlo. *Statistical applications in genetics and molecular biology*, 12(1):87–107.

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*.

Gloeckler, M., Deistler, M., Weilbach, C., Wood, F., and Macke, J. H. (2024). All-in-one simulation-based inference. *ArXiv*, abs/2404.09636.

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.

Hyndman, R. J. (1996). Computing and graphing highest density regions. *The American Statistician*, 50(2):120–126.

Innes, M. (2018). Don't unroll adjoint: Differentiating SSA-form programs. *CoRR*, abs/1810.07951.

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.
Kendall, A. and Gal, Y. (2017). What uncertainties do we need in bayesian deep learning for computer vision? *Advances in neural information processing systems*, 30.

Kompa, B., Snoek, J., and Beam, A. L. (2021). Empirical frequentist coverage of deep learning uncertainty quantification procedures. *Entropy*, 23(12).

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.

Lindsay, B. (1988). *Composite likelihood methods*. Statistical Inference from Stochastic Processes. American Mathematical Society, Providence.

Macci, C. and Spizzichino, F. (2015). What about the posterior distributions when the model is non-dominated? In Polpo, A., Louzada, F., Rifo, L. L. R., Stern, J. M., and Lauretto, M., editors, *Interdisciplinary Bayesian Statistics*, pages 1–11, Cham. Springer International Publishing.

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.

Métivier, D. (2025). ABCMethods: A Julia package implementing some Bayesian ABC methods.
Neal, R. M. (2012). *Bayesian learning for neural networks*, volume 118. Springer Science & Business Media.

Pal, A. (2025). *Lux: Explicit Parameterization of Deep Neural Networks* in Julia. Zenodo.

Papamakarios, G. and Murray, I. (2016). Fast $\epsilon$-free inference of simulation models with bayesian conditional density estimation.

Papamakarios, G., Sterratt, D., and Murray, I. (2019). Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In Chaudhuri, K. and Sugiyama, M., editors, *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics*, volume 89 of *Proceedings of Machine Learning Research*, pages 837–848. PMLR.

Picchini, U. and Tamborrino, M. (2024). Guided Sequential ABC Schemes for Intractable Bayesian Models. *Bayesian Analysis*, pages 1–32.

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.

Rackauckas, C. and Nie, Q. (2017). *DifferentialEquations.jl* – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. *Journal of Open Research Software*, 5(1):15–15.

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.

Ribatet, M., Cooley, D., and Davison, A. (2012). Bayesian inference from composite likelihoods, with an application to spatial extremes. *Statistica Sinica*, 22:813–845.

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.

Stoehr, J. and Friel, N. (2015). Calibration of conditional composite likelihood for bayesian inference on gibbs random fields. *Journal of Machine Learning Research: Workshop and Conference Proceedings*, 38:921–929.

Van der Vaart, A. (1998). Asymptotic statistics. cambridge series in statistical and probabilistic mathematics.

Varin, C., Reid, N., and Firth, D. (2011). An overview of composite likelihood methods. *Statistica Sinica*, 21(5–42).

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.

Wilkinson, D. (2024). smfsb: Stochastic modelling for systems biology.

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 Application: Moving Average 2

Complementing Section 4.1 of the article, additional figures are given for illustration purposes. Note that all these figures can be found in Notebooks associated with this publication.

* Figure A.1 shows, for a test sample, the true parameter $\theta$, as well as estimated $\hat{\theta}$ using standard ABC, ABC-CNN and ABCD-conformal and associated confidence ellipses.

* Figure A.2 illustrates that the true value of $\theta$ is close to the associated exact posterior mean. This last one is calculated from explicit formulas from Marin and Robert [2007][Chapter 5] and requires the inversion of a large matrix, which is numerically costly and a source of possible additional numerical error.

# B Application: Lotka Volterra

Complementing Section 4.2 of the article, additional figures are given for illustration purposes. Note that all these figures can be found in Notebooks associated with this publication.

* Figure B.1 shows for the test sample 997, estimated parameter $\theta$ with associated 3D confidence ellipsoids obtained using standard ABC, ABC-CNN and ABCD-conformal, and true value of $\theta$.

* Figure B.2 shows for the test sample 667, estimated parameter $\theta$ with associated 3D confidence ellipsoids obtained using standard ABC, ABC-CNN and ABCD-conformal, and true value of $\theta$.

# C Application: Lake Ecosystem dynamics

## C.1 Detailed description of the model

The system we are considering here corresponds to a lake ecosystem for which the time evolution of the phytoplankton and the inorganic phosphorus concentrations is represented through the following equations:<fnref n="1" />

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

$$ \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 (17) $$

where $C_P$ and $C_{IP}$ respectively represent the phytoplankton concentration [mg C.L⁻<fnref n="1" />] and the inorganic phosphorus concentration [mg P.L⁻<fnref n="1" />], $C_{OP}$ is the organic phosphorus concentration

<footnote n="1">¹mg, P, C and L stand for milligram, phosphorus, carbon and liter respectively.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure A.1: Moving Average 2 example: for the 77th test sample, the true value of $\theta$ is represented (up-triangle), as well as the estimated $\hat{\theta}$ using standard ABC, ABC-CNN and ABCD-conformal and associated confidence ellipses (circles and ellipses). The exact posterior mean is also represented (down triangle), for comparison with the true value.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure A.2: Moving Average 2 example: each red (resp. green) dot represents the exact posterior mean $\mathbb{E}(\theta_1|\mathbf{x})$ (resp $\mathbb{E}(\theta_2|\mathbf{x})$) in the x-axis, against the true value of $\theta_1$ (resp. $\theta_2$) in the y-axis, for each of the 1000 test samples.</figcaption></figure>

$[\text{mg P.L}^{-1}]$ that is assumed to evolve slowly compared to the other variables, and that is given, due to the law of conservation of mass, by:

$$C_{\text{OP}} = \varepsilon(K_0 - a_{\text{pc}}C_{\text{P}} - C_{\text{IP}}) \quad (18)$$

where $K_0 = a_{\text{pc}}C_{\text{P},0} + \frac{1}{\varepsilon}C_{\text{OP},0} + C_{\text{IP},0}$ [$\text{mg P.L}^{-1}$] represents the initial total phosphorus concentration in the water (which is conserved) and $C_{\text{P},0}$ [$\text{mg C.L}^{-1}$], $C_{\text{IP},0}$ [$\text{mg P.L}^{-1}$] and $C_{\text{OP},0}$ [$\text{mg P.L}^{-1}$] are the initial concentrations of phytoplankton, inorganic phosphorus and organic phosphorus in the water.

In the model, the parameter $d_{\text{P}}$ represents the death rate of phytoplankton [$\text{day}^{-1}$], $f_{\text{OP}}$ is the fraction of dead phytoplankton going to the organic phosphorus pool [-], $a_{\text{pc}}$ is the phosphorus-to-carbon ratio of the phytoplankton [$\text{mg P.(mg C)}^{-1}$] and $\varepsilon$ [-] is an adimensional parameter reflecting the slow evolution of organic phosphorus compared to the other variables.

The time-dependent temperature of the water $T$ is assumed to be a sinusoidal function
<figure><img src="image_5.png" /><figcaption>Figure B.1: Lotka-Volterra example: for the test sample 997, estimated parameter $\theta$ with associated 3D confidence ellipsoids obtained using standard ABC, ABC-CNN and ABCD-conformal, and true value of $\theta$.</figcaption></figure>

of the form:

$$T(t) = \left( \cos \left( \frac{\pi}{f} t + \pi \right) + 1 \right) \frac{T_{\max} - T_{\min}}{2} + T_{\min} \quad (19)$$

where $T_{\max}$ [° C] and $T_{\min}$ [° C] are the maximal and minimal temperatures reached by the water temperature during the year and $f$ [day] is the period of the sinusoidal function.

The growth rate of the phytoplankton $\mu$ is expressed by:

$$\mu(C_{\text{IP}}, T) = \mu_{\max} \frac{C_{\text{IP}}}{K_{\text{IP}} + C_{\text{IP}}} \theta_{gr}^{T-20} \quad (20)$$

where $\mu_{\max}$ is the maximal growth rate of the inorganic phosphorus uptake [day$^{-1}$], $K_{\text{IP}}$ is the half saturation constant of the phytoplankton growth function [mg P.L$^{-1}$], and $\theta_{gr}$ is the temperature coefficient for phytoplankton growth [-].

Finally the mineralization speed of the organic phosphorus into inorganic phosphorus $V_{\min}$ is given by:

$$V_{\min}(C_{\text{OP}}, C_{\text{IP}}, T) = V_{\max} \frac{C_{\text{OP}}}{K_{\text{OP}} \left(1 + \frac{C_{\text{IP}}}{K_{\text{IP}}}\right) + C_{\text{OP}}} \theta_{\min}^{T-20} \quad (21)$$

with $V_{\max}$ the maximal mineralization rate [mg P.L$^{-1}$.day$^{-1}$], $K_{\text{OP}}$ the half saturation constant of the mineralization function [mg P.L$^{-1}$], $K_{\text{I}}$ the inhibition coefficient of the enzymatic reaction [mg P.L$^{-1}$] and $\theta_{\min}$ the temperature coefficient for mineralization [-].
<figure><img src="image_6.png" /><figcaption>Figure B.2: Lotka-Volterra example: for the test sample 667, estimated parameter $\theta$ with associated 3D confidence ellipsoids obtained using standard ABC, ABC-CNN and ABCD-conformal, and true value of $\theta$.</figcaption></figure>

In the model, there are 15 parameters and 2 initial conditions. Among the 15 parameters, the following 6 parameters have a fixed value and are not included in the estimation process:

$$a_{\text{pc}} = 0.00625 \, [\text{mg P.}(\text{mg C})^{-1}], \quad T_{\text{min}} = 5 \, [^{\circ}\text{C}], \quad T_{\text{max}} = 30 \, [^{\circ}\text{C}], \\ f = \frac{365}{2} \, [\text{day}], \quad \varepsilon = 10^{-6} \, [-], \quad C_{\text{OP},0} = 1 \, [\text{mg P.L}^{-1}].$$

The objective is to estimate the 9 remaining parameters that is the vector $\theta$ given by: $\theta = (\mu_{\text{max}}, K_{\text{IP}}, d_{\text{P}}, f_{\text{OP}}, V_{\text{max}}, K_{\text{OP}}, K_{\text{I}}, \theta_{\text{min}}, \theta_{\text{gr}})$

The initial conditions are fixed to $C_{\text{P},0} = 0.01 \, [\text{mg P.L}^{-1}]$, and $C_{\text{IP},0} = 0.01 \, [\text{mg P.L}^{-1}]$. A dataset corresponds to observations of state ($C_{\text{P}}, C_{\text{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 $U[\theta_i^m, \theta_i^M]$ for each parameters $\theta_i$ of $\theta$, with $\theta^m = (1, 0.0009, 0.03, 0, 0.05, 0.03, 0, 1, 1)$ and $\theta^M = (3, 0.052, 0.7, 1, 0.22, 3, 0.04, 1.1, 1.1)$.

## C.2 Additional figures

Complementing Section 4.3 of the article, additional figures are given for illustration purposes. Note that all these figures can be found in Notebooks associated with this publication.
Figures C.1, C.2, C.3 and C.4 show for each component of $\theta$, the estimated values against the true value for the 1000 test samples, as well as associated confidence sets. The estimations are obtained using standard ABC, ABC-SMC, ABC-CNN and ABCD-conformal respectively.

<figure><img src="image_7.png" /><figcaption>Figure C.1: Lake ecosystem dynamics example: for each component of normalized $\theta$ denoted $\tilde{\theta}$, estimated values against true values, for 1000 test samples, as well as associated confidence intervals (in gray, for the overall uncertainty). The estimated values are obtained using standard ABC.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure C.2: Lake ecosystem dynamics example: for each component of normalized $\theta$ denoted $\tilde{\theta}$, estimated values against true values, for 1000 test samples, as well as associated confidence intervals (in gray, for the overall uncertainty). The estimated values are obtained using ABC-SMC.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure C.3: Lake ecosystem dynamics example: for each component of normalized $\theta$ denoted $\tilde{\theta}$, estimated values against true values, for 1000 test samples, as well as associated confidence intervals (in gray, for the overall uncertainty). The estimated values are obtained using ABC-CNN.</figcaption></figure>

# D Application: Gaussian random field

## D.1 The model

We study stationary isotropic Gaussian random field on the domain $[0, 5] \times [0, 5]$ with a regular grid of size 100×100, with exponential covariance functions [Wood and Chan, 1994].
<figure><img src="image_10.png" /><figcaption>Figure C.4: Lake ecosystem dynamics example: for each component of normalized $\theta$ denoted $\tilde{\theta}$, estimated values against true values, for 1000 test samples, as well as associated confidence intervals (in gray, for the overall uncertainty). The estimated values are obtained using ABCD-conformal.</figcaption></figure>

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.

## D.2 Algorithm parametrization

Number of used samples are given in Table 1 of the article. For 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-standardized, 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 ABCD-conformal, we used 3 convolutional 2D layers having 32, 64 and 64 neurons, kernels of size 3 and the **relu** activation function, followed by max-pooling for the first two layers and by a flatten layer for the third one. Then, 2 dense layers having both 64 neurons and the **relu** activation function. The last layer to obtain the output of dimension 1 is linear. 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$.

Concerning ABCD-conformal, we look at its performances using two different measures of uncertainty: the overall uncertainty and the epistemic uncertainty, both returned by the Dropout procedure.

**Remark 4.** This example is the only one where we use the ABC-RF which can be applied only on uni-dimensional parameters. For the other examples with multidimensional parameters, we only considered methods that can deal with multidimensionality.

**Remark 5.** In this example, we do not show results of ABC-SMC. The ABC-SMC was implemented, but using *R* the algorithm took too long to be run for a sufficient number of iterations and satisfactory population sizes. Indeed, ABC-SMC is the only method that do not rely on a reference table and by consequence requires many simulations for each test sample, see section E.1 for a discussion about that.

## D.3 Results

All indicators detailed in the beginning of Section 4 of the article are presented in Table 7. Figure D.1 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 7, we see that in this example, ABC-RF outperforms the other methods, with the smallest NMAE and $sd(|\theta - \hat{\theta}|)$, while having a very good coverage, 98%, which is better than the expected coverage of 95%. The standard ABC and ABC-CNN do not give
<table><caption>Table 7: Two dimensional Gaussian random field: for each method, we give the NMAE, $sd(|\theta - \hat{\theta}|)$, the mean and median lengths of the confidence or credible intervals and the coverage. These indicators are computed on a test set of size $10^2$. The best values for each line are highlighted in bold.</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><strong>0.0151</strong></td><td>0.0300</td><td>0.0300</td></tr><tr><td>sd(|θ - ̂θ|)</td><td>0.0229</td><td>0.0106</td><td><strong>0.0076</strong></td><td>0.0141</td><td>0.0141</td></tr><tr><td>mean length CI(θ)</td><td>0.1313</td><td><strong>0.0307</strong></td><td>0.0405</td><td>0.0684</td><td>0.0711</td></tr><tr><td>median length CI(θ)</td><td>0.1341</td><td><strong>0.0298</strong></td><td>0.0391</td><td>0.0637</td><td>0.0672</td></tr><tr><td>coverage CI(θ)</td><td><strong>100%</strong></td><td>88%</td><td>98%</td><td>93%</td><td>94%</td></tr></tbody></table>

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 confidence intervals. In the case of the ABC-CNN, like in the previous example, the NMAE, $sd(|\theta - \hat{\theta}|)$ and the lengths of the confidence sets are quite small, but it is counterbalanced by a too small coverage (88%). Finally, the results given by the ABCD-Conformal approach are satisfactory, with coverages close to 95%. 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 coverages. Standard ABC can be discarded compared to other methods in view of its poor results, except for the coverage.

# E Pseudocodes for ABC-SMC, ABC-CNN and ABC-RF and some comments

## E.1 ABC-SMC

Like the standard ABC algorithm, the ABC-SMC requires summary statistics and a distance to compare two data samples. While standard ABC only requires a single threshold, ABC-SMC needs a sequence of decreasing thresholds (or a way to compute it) at each iteration. Other hyperparameters include the size of the sequential populations $N$, the number of iterations $T$, as well as perturbation samplers (transition kernels) $K$ used to move from one population to another. Therefore, compared to standard ABC, ABC-SMC involves many more calibration parameters. Generally speaking, ABC-SMC yields better results than standard ABC, but requires a bigger expertise.

That said, assuming an informed user knows how to set up an ABC-SMC algorithm, one major difference between ABC-SMC and all the other methods used in this paper is that simulations must be carried out sequentially. In contrast, standard ABC, ABC-CNN, ABC-
<figure><img src="image_11.png" /><figcaption>Figure D.1: Two-dimensional Gaussian random field example: 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>
RF, and ABCD-conformal benefit from the use of a reference table, which allows simulations to be generated in advance, potentially through parallelization. Additionally, the same reference table can be reused to estimate parameters for different data samples of interest. This enables these algorithms to be calibrated and tested much more efficiently. For example, for a test set of 1000 samples, it is not necessary to generate a new reference table for each test sample while the ABC-SMC, on the other hand, requires a new set of simulations for each test sample.

The version of the ABC-SMC algorithm used for the numerical simulations is given below in Algorithm 3. Its implementation is very similar to that of the R package smfsb Wilkinson [2024]. There are multiple ways to define a stopping criterion for the algorithm, for instance, one can choose to fix the total number of simulations performed during the process, to prevent excessively long runs (this was our choice). However, in this case, the final tolerance level is not guaranteed. Conversely, the user may decide to fix a final tolerance level that needs to be reached, but in this case, the total number of simulations is not known in advance. This is just one example; many different versions of this algorithm have been proposed since Sisson et al. [2007]; Beaumont et al. [2009] and Del Moral et al. [2012].

To complement the Algorithm 3, here are a few remarks:

* Regarding the transition kernel $\mathcal{K}$ used to perturb a random draw from a previous population (line 10 of Algorithm 3), we chose the classic Gaussian kernel (e.g., Filippi et al. [2013]). At iteration $t$, this kernel is centered on a random draw from the weighted population ($t-1$) and has a covariance matrix equal to the empirical weighted covariance matrix of this population, $\Sigma_{t-1}$. The importance weights are then given by (with the usual abuse of notation): $\tilde{w}_i^{(t)} = \pi(\theta_i^{(t)}) / \sum_{j=1}^N w_j^{(t-1)} \mathcal{N}(\theta_i^{(t)}; \theta_j^{(t-1)}, 2\Sigma_{t-1})$, where $\mathcal{N}(\theta; m, \Sigma)$ denotes the probability density function of a Gaussian distribution with mean $m$ and covariance matrix $\Sigma$ evaluated at $\theta$. Note that we experimented with the popular choice of using a factor of 2 for the covariance matrix [Filippi et al., 2013], with no clear advantage.

* At line 17 of Algorithm 3, the new set $\theta^{(t)}$ is, by construction, of size $N$. In the case of ties for the $N$-th smallest distance, we arbitrarily select from the associated particles $\theta_j$ to ensure the set contains exactly $N$ particles. For each iteration $t$, this way of selecting particles corresponds to a tolerance level $\epsilon_t$ defined by $\Delta_{(N)}$.

* For fair comparisons with other methods, we ensure that the total number of calls to the model $f$ is equal to $N_{train}$, i.e., $N \times \Psi \times T = N_{train}$. Note that this means we must generate $N_{train}$ samples from the model for every test sample, which is significantly more time-consuming than other methods that use a common reference table.

* To guarantee fairness in the Lotka-Volterra example (Section 4.2 of the article), we require the draw at line 12 of Algorithm 3 to yield a valid, non-zero trajectory. If an extinction occurs, the draw is discarded and we return to line 9 of the algorithm.
**Algorithm 3:** Pseudocode for the ABC-SMC sampler [Sisson et al., 2007; Beaumont et al., 2009; Del Moral et al., 2012]

**Input**: A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, the number of samples kept per iteration $N$, the number of steps $T$, summary statistics $\eta: \mathbb{R}^d \rightarrow \mathbb{R}^l$, a (pseudo)-distance $d$ on $\mathbb{R}^l$, a perturbation sampler $\mathcal{K}_t$ and an integer factor $\Psi > 0$.

**Output**: A weighted set $(\theta_1^{(T)}, \tilde{w}_1^{(T)}), \dots, (\theta_N^{(T)}, \tilde{w}_N^{(T)})$ whose distribution is approximately $\pi(\theta | \mathbf{x})$

1 **Initial population**:

2 **for** $i \leftarrow 1$ **to** $N$ **do**

3 Set $\theta_i^{(0)} \sim \pi$;

4 Set $w_i^{(0)} := 1/N$.

5 **end**

6 **Sequential populations**:

7 **for** $t \leftarrow 1$ **to** $T$ **do**

8 **for** $j \leftarrow 1$ **to** $N \times \Psi$ **do**

9 Sample with replacement $\theta_j^*$ from the weighted set $\{\theta_i^{(t-1)}, w_i^{(t-1)}\}_{i=1,...,N}$;
10 Perturb $\theta_j^*$: $\theta_j^{**} \sim \mathcal{K}_{t-1}(\cdot | \theta^*)$;
11 if $\pi(\theta_j^{**}) = 0$ go to line 10, otherwise continue;
12 Draw synthetic sample $\mathbf{x}_j^{**}$ from the model $f(\cdot | \theta_j^{**})$;
13 Compute the summary statistics $\eta(\mathbf{x}_j^{**})$;
14 Compute the distance $\Delta_j = d(\eta(\mathbf{x}_j^{**}), \eta(\mathbf{x}))$;

15 **end**

16 Let $\Delta_{(1)} \le \Delta_{(2)} \le \dots \le \Delta_{(N \times \Psi)}$ be the sorted distances;

17 Set the new population $\{\theta_i^{(t)}\}_{i=1}^N$ to be the particles $\theta_j^{**}$ for which $\Delta_j \le \Delta_{(N)}$;

18 Set $\tilde{w}_i^{(t)} = \pi(\theta_i^{(t)}) / \sum_{j=1}^N w_j^{(t-1)} \mathcal{K}_{t-1}(\theta_i^{(t)} | \theta_j^{(t-1)})$;

19 Normalize the weights $w_i^{(t)} := \tilde{w}_i^{(t)} / \sum_{j=1}^N \tilde{w}_j^{(t)}$;

20 **end**
## E.2 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 in the article) 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 4. The output of the algorithm is similar to the one of the standard ABC: for a given data sample **x** and a tolerance threshold, we hope that the parameters **θ** accepted during the process are approximately distributed from the posterior distribution π(**θ** | **x**). Hence, to have asymptotic consistency and valid frequentist coverages, the same kind of conditions are needed than for standard ABC (Algorithm 1 in the article). 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 α 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 standard ABC (Algorithm 1 in the article). However, this two-stage procedure does reduce the Monte Carlo error inherent in estimating the approximate posterior distribution π<sub>α,d</sub>(.|η(<u>x</u>)) 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.

Another difference with the standard ABC is that a validation dataset should be gen-
erated, 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 sam-
ples in the reference table (lines 9 of Algorithm 1 and 15 of Algorithm 4). This comparison
can be time-consuming, depending on the distance used between summary statistics and the
**Algorithm 4:** Pseudocode 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;&nbsp;Draw $\theta_j \sim \pi$;

4. &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)$;

5. **end**

6. **Generation of a validation dataset:**

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

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

9. &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)$;

10. **end**

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

12. **Summary statistics and distances:**

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

14. &nbsp;&nbsp;&nbsp;&nbsp;The vector $\mathbf{x}_j$ is given as input to the trained CNN to obtain $\eta_{\text{CNN}}(\mathbf{x}_j)$;

15. &nbsp;&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.
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 generalize, and it clearly depends on the problem at hand.

## E.3 ABC-RF

Raynal et al. [2018] 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. Their approach relies on the random forest (RF) methodology of Breiman [2001] applied in a regression setting.

<table><caption>**Algorithm 5:** Pseudocode for an Approximate Bayesian Computation - Random Forest sampler</caption><tr><td>Input</td><td>: A Bayesian parametric model {$f(\cdot | \theta), \pi$}, a data sample $\mathbf{x}$, $N_{\text{train}}$ the size of the training set, summary statistics $\eta: \mathbb{R}^d \to \mathbb{R}^l$.</td></tr><tr><td>Output</td><td>: For $T(\theta)$ a scalar transform of $\theta$, approximations of the posterior expected value $\psi_{\eta}(\mathbf{x}) = E_{\pi}[T(\theta) | \eta(\mathbf{x})]$, the posterior variance $V_{\pi}[T(\theta) | \eta(\mathbf{x})]$, and posterior quantiles.</td></tr><tr><td>1</td><td>Generation of a reference table</td><td>:</td></tr><tr><td>2</td><td>for $j \leftarrow 1$ to $N_{\text{train}}$ do</td><td></td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Draw $\theta_j \sim \pi$;</td><td></td></tr><tr><td>4</td><td>&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)$</td><td></td></tr><tr><td>5</td><td>end</td><td></td></tr><tr><td>6</td><td>Summary statistics</td><td>:</td></tr><tr><td>7</td><td>for $j \leftarrow 1$ to $N_{\text{train}}$ do</td><td></td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Compute the vector of summary statistics $\eta(\mathbf{x}_j)$</td><td></td></tr><tr><td>9</td><td>end</td><td></td></tr><tr><td>10</td><td>Train a random forest</td><td>: the training pairs are the {$(\mathbf{x}_j, \theta_j), j = 1, \dots, N_{\text{train}}$}.</td></tr><tr><td>11</td><td>Approximate $\psi_{\eta}(\mathbf{x})$</td><td>: give $\eta(x)$ as input of the trained random forest: the output of the RF allows to obtain easily approximations of $\psi_{\eta}(\mathbf{x}) = E_{\pi}[T(\theta) | \eta(\mathbf{x})]$, $V_{\pi}[T(\theta) | \eta(\mathbf{x})]$ and of posterior quantiles (hence credible intervals).</td></tr></table>

This algorithm automatizes 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(\theta)$ given $\eta(\mathbf{x})$, $F(T(\theta) | \eta(\mathbf{x}))$. The asymptotic consistency of this approach has been established by
Meinshausen [2006], under some conditions. In particular, $F(T(\theta) | \eta(\mathbf{x}))$ should be Lipschitz continuous and strictly monotonously increasing in $T(\theta)$. Such conditions are practically unverifiable in an ABC framework where the likelihood is unknown or intractable.

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).

# F Detailed comparison of the five ABC methods used in the examples

Below the five methods used (standard ABC, ABC-SMC, ABC-CNN, ANC-RF and ABCD-Conformal) are compared in details on different points, with a practical point of view. Table 6 in the article summarizes these reflections.

## Output of the algorithm

First, one big difference between these methods is what is obtained from the algorithm: for standard ABC, ABC-SMC and ABC-CNN an approximation of the whole posterior distribution of the parameter of interest is obtained. The ABC-RF and ABCD-Conformal 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 statistics

Standard ABC and ABC-SMC 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 of the article, 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]. ABC-RF enables to automatize 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 ABCD-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 and ABC-SMC. 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 Lotka-Volterra example: we did not really use summary statistics, as datasets and simulations were directly compared using an Euclidian distance. But as the dataset was of dimension forty (two series of twenty times), it was computationally demanding: the standard ABC was the method which took more computing time. Regarding 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. 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 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-SMC requires a decreasing sequence of tolerance levels, or at least a first threshold and a way to decrease it between each population (usually using percentiles). Like for standard ABC and ABC-CNN, the choice of this sequence can affect a lot the degree of approximation obtained. On the opposite, 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 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 and the ABC-SMC 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 as the parameter of interest. For ABCD-Conformal, multidimensionality is not a problem, as the computing time for the CNN with Dropout and the conformal procedure does not really increase with the dimension (at least on the examples studied).

## 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 layer, 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. The ABCD-Conformal 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. Concerning the ABC-SMC, it cannot rely on datasets generated in advance.
Simulations must be performed sequentially, and are specific to each data sample of interest.
In practice, ABC-SMC therefore requires significantly more time for calibration and testing
compared to the other methods.

## Theoretical justification of the methods

In this article, we are particularly interested in 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 four examples,
and we studied the bibliography for the theoretical justification of all these methods.

As detailed in section 2.1 of the article, Frazier et al. [2018] studied asymptotic proper-
ties 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 ABC-SMC, the convergence of importance sampling-based algorithms has
been extensively studied, as many different algorithms have been proposed, see for instance
Beaumont et al. [2009]; Del Moral et al. [2012]. Convergence requires specific conditions,
in particular regarding the sequence of tolerance thresholds, and the choice of perturbation
samplers (transition kernels).

As said in section 2.3.2 in the article, ABC-RF 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 samples
should solely come from the same distribution as the calibration samples. 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].

Benefits of using concrete Dropout

Apart from the fact that epistemic and aleatoric uncertainties are used to obtain confidence
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?