Multi-objective Bayesian optimization for Likelihood-Free inference in sequential sampling models of decision making
Authors: David Chen
Xinwei Li
Eui-Jin Kim
Prateek Bansal
David J Nott
Venue: Accepted by TMLR ()

Abstract: Statistical models are often defined by a generative process for simulating synthetic data, but this can lead to intractable likelihoods. Likelihood free inference (LFI) methods enable Bayesian inference to be performed in this case. Extending a popular approach to simulation-efficient LFI for single-source data, we propose Multi-objective Bayesian Optimization for Likelihood Free Inference (MOBOLFI) to perform LFI using multi-source data.  MOBOLFI models a multi-dimensional discrepancy between observed and simulated data, using a separate discrepancy for each data source. The use of a multivariate discrepancy allows for approximations to individual data source likelihoods in addition to the joint likelihood, enabling detection of conflicting information and deeper understanding of the importance of different data sources in estimating individual parameters. The adaptive choice of simulation parameters using multi-objective Bayesian optimization ensures simulation efficient approximation of likelihood components for all data sources.  We illustrate our approach in sequential sampling models (SSMs), which are widely used in psychology and consumer-behavior modeling.  SSMs are often fitted using multi-source data, such as choice and response time.  The advantages of our approach are illustrated in comparison with a single discrepancy for an SSM fitted to data assessing preferences of ride-hailing drivers in Singapore to rent electric vehicles.

Published in Transactions on Machine Learning Research (06/2025)
Multi-objective Bayesian optimization for Likelihood-Free in- ference in sequential sampling models of decision making
David Chen
e e10396889a. ndu
Department of Statistics and Data Science National University of Singapore 117546, Singapore
Xinwei Li
Sinwei.Wiaa edu
Department of Civil and Environmental Engineering National University of Singapore 117576, Singapore
Eui-Jin Kim
euijin@ajov.ac.. kr.
Department of Transportation Systems Engineering
Ajou University
16499, Korea
Prateek Bansal
prateckb @rous.edee edu.sg sg
Department of Civil and Ennirommental Engineering National University of Singapore 117576, Singapore
David Nott
Blandj@ruus.. www.errwuuuu.....
Department of Statistics and Dala Science National Universily of Singapore 117546, Singapore
Reviewed on OpenReview: https: // openreview net/ forum? id- hQjwDq (Saj
Abstract
Statistical models are often defined by a generative process for simulating synthetic data, but this can lead to intractable likelihoods. Likelihood free inference (LFI) meth- ods enable Bayesian inference to be performed in this Extending popular approach to simmlation-ccficient LFI for single-source data, propose Multi-objective Bayesian Op- linization for Likelihood Free Inference (MOBOLFI) to perform LFI using multi-source data. MOBOLFI models lalti-dimennsonal discrepancy between observed and simulated data, using a separate discrepancy for each data source. The use of a multivariate dis- crepancy allows for approximations to individual data source likelihoods in addition to the joint likelihood, enabling detection of conflicting information and deeper understanding of the importance of different data sources in estimating individual parameters. The adaptive choice of simulation parameters using multi-objective Bayesian optimization ensures simula- lion elficient approximation of likelihood components for all data sources. We illustrate our approach in sequential sampling models (SSMs), which are widely used in psychology and consumer-belavior modeling, SSMs are often fitted using multi-source data, such as choice and response time. The advantages of our approach are illustrated in comparison with a single discrepancy for an SSM fitted to data assessing preferences of ride-hailing drivers in Singapore to rent electric vehicles.
*David Chen and Xinwei contributed equally this work.

---

Published in Transactions on Machine Learning Research (06/2025)
1 Introduction
In many scientific applications, the most natural way to incorporate subject matter knowledge into statistical model is by specifying generative algorithm (or simulator) for synthetic data. In this case, the likelihood function can sometimes be difficult to compule, and this makes (raditional statistical inference methods hard to apply. In response to this challenge, researchers have developed likelihood-free inference (LFT) methods for conducting Bayesian inference using only simulation from the model, without any likeli- hood evaluations. An alternative name for LFI is simulation-basse inference (SBI), which stresses the role of simulating synthetic data. The phrase "likelihood-free should not be taken to mean that there is no likelihood. Rather, it means that approximate likelihoods or approximate posterior sampling methods can be constructed from model simulations, without evaluating the likelihood.
Our work considers settings where it is important to obtain accurate approximations to the posterior density using as few simulations as possible, and where the computationally intractable likelihood function factorizes into terms for multiple data sources. We extend a widely-used simulation efficient method for single source data, Bayesian optimization for likelihood-free inference (BOLFT) (Gutmann & Corander 2016). BOLFI uses Bayesian optimization to sequentially choose the simulation paramcters to make best use of limited computational budget. In the case of multi-source data, we consider a multi-objectivve extension which call MOBOLFI, which measures similarity between synthetic and observed data by a multi-dimensional discrepancy with a component for each data source. The multi-dimensional discrepancy is modelled, and this allows approximation of likelihood terms for different data, which is important for understanding the importance of different data sources for inference about individual parameters and for model checking. While BOLFI provides a closed form approximate likelihood for inference tasks, it does not provide approximations of likelihood terms for individual data sourees for multi-source data. In the multi- source setting, it is often natural to choose summary statistics and discrepancy metrics separately for each data source, and with the classic BOLFI approach we then have to combine all discrepancies into single one. If this is not done carefully, the information in the combined discrepancy may not reflect correctly the relative importance of different data sources, resulting in information loss and highly conservative uncertainty quantification.
There are many existing methods for LFI. One of the most well-established methods is approximate Bayesian computation (ABC, Sisson et al., 2018), which approximates the posterior by parameter samples for which a simulated dataset has a discrepancy from the observed dataset below a certain tolerance. See Scction 2.1 for more detailed discussion. ABC implicitly uses kernel density approximation to the likelihood for data summaries. A competitive approach, probability density approximation (PDA, Turner & Sederberg 2014), has also been widely used for over a decade. However, with high-dimensional summary statistics, both ABC and PDA can perform poorly, making then unsuitable for many highly parameterized models.
Structured density estimation approaches which deal better with high-dimensional summary statistics include synthetic likelihood, which is based working normal model for summary statistics (Wood] 2010:
Frazier et al. 2022). flexible extensions of synthetic likelihood based on copulas empirical saddle point
approximation Fasiolo et al. 2018 An et al.) 2020) and neural density estimation approaches (Papamakarios
& MMrray 2016 Papamakarios et al., 2019 Greenberg et al. 2019). However, these structured density estimation methods require a large number of model simulations to obtain adequate approximations, which may be infeasible if simulation from the model is computationally expensive. A preferred approach in such settings is to use active lcarning approaches such as Bayesian optimization for Likelihood inference (BOLFI, Gutmaan & Corander 2016), and extending these methods is the focus of our work. We review the BOLFI approach in Section 2.3 and related existing literature in Section A.2 of the Appendix.
Another challenge for LFI methods is model misspecification (Wilkinson 2013 Frazier et al., 2020 Frazier & Drovandi 2021 Ward et. al. 2022). There is an extensive literat ure on this issue, but here we focus on model checking with multi-source data, where the different data sources may bring conflicting information about the model parameter (Presanis et al. (203)). Modern LFI mothods have been used in the context of Sequential Sampling Models (Radev et al. [2023) and for multi-source data (Schmitt et al. 20233] Bahg et al., 2020), but questions of model adequacy have been previously addressed in global way. The MOBOLFT
2

---

Published in Transactions on Machine Learning Research (06/2025)
method proposed here has the capability to approximate both the joint likelihood as well as likelihoods for individual data sources, enabling checks of the consistency of information supplied by different parts of the data. Examining this consistency is essential for trustworthy analysis of complex multi-source data such as those arising for Sequential Sampling Models (SSMs). It is also valuable for understanding which data sources are most informative for the estimation of individual SSM parameters.
This work makes three main contributions. First, the multi-objective aspect of our MOBOLFI method ensures efficient exploration of the high-likelihood region for the separate likelihoods of all data sources. Importantly, in the case of computationally expensive simulator, the MOBOLFI approach can explore high-likelihood regions for all data source likelihoods in single run. This is particularly difficult to do with other approaches when the high-likelihood regions for individual data source likelihoods are disjoint. Our second main contribution relates to the difficulty of combining multiple discrepancies from individual data sources into a single discrepancy, with consequent information loss and highly conservative uncertainty quantification poor choices are made. MOBOLFI avoids the need to do this, making it casier to use in practice. Finally, MOBOLFI is able to leverage likelihood approximations for individual data sources to check their consistency and understand their importance for estimating individual SSM parameters.
The paper is organized as follows. Section 2 gives background on Bayesian optimization and the BOLFI method of Gutmann & Corander (2016). Section then discusses multi-objective optimization, the motiva- tion for the MOBOLFI method, and its implementation. Section 4 presents two instructive examples, one involving a Brownian motion and the second involving a SSM. These examples illustrate the advantages of MOBOLFI compared to the BOLFI method with single discrepancy in settings involving multi-source data. In the second example, use an empirical choice-response time (choice-RT) dataset consumer preferences for electric vehicles, and demonstrate MOBOLFT's application in establishing the importance of different data sources for estimation of parameters of an SSM. Key takeaways and avenues for future research are discussed in section 5
2 Method
Before explaining the MOBQLFI method, we give some necessary background on Bayesian optimization, ABC and the traditional BOLFI approach. Describing our method requires introducing some complex notation in Sections 2 and 3. To help the reader, a glossary of notation given in Appendix A.
2.1 Approximate Bayesian computation
ABC methods approximate the likclihood by the probability of a synthetic datasct being close to the observed data in terms of some discrepancy. The discrepancy is often defined in terms of data summary statistics. It is important to explain the ABC likelihood approximation, since this motivates the BOLFI approach discussed in subscction 2.3 For more extensive discussion of ABC than we give here see Sisson et al. 2018).
Let be parameters in a statistical model for data y E y to be observed. Write yobs for the observed value and p(yor=l6) for the likelihood, where p(y]@) the density of y/O. In ABC, inference is usually basod on a diwer-dieeniional summary of the original data of dimension d say, defined by a mapping S y Rd. The observed summary statistic value S(yobs) will be denoted by Sobs-
In ABC we first approximate the posterior density p(0)Mobs)) & p(0)p(uora(0)) by the partial posterior density which conditions on Sobs rather than Yobs,
S(0)S x ((0)p(SStmml))
If the summary statistics S are sufficient, then there is no loss in replacing Bobs with Sabes but non-trivial sufficient statistics will not exist in most complex models of interest, and we only require summary statistics to be informative about 0. Next, since pl p(S (0) is infeasible to compute if p(yota is, we replace p(Sotaa
3

---

Published in Transactions on Machine Learning Research (06/2025)
with the "ABC likelihoodd denoted (r (0) for some tolerance t 0:
(s ((S)))((A6(SS55 ds,
(1)
where Sots) 1 ||S Sobell and ||-II is some distance measure. Ao Sobs) is a mcasure of the discrepancy between the simulated and observed summary stalistics. As discussed further below, the ABC likelihood at can be thought of as proportional to the probability that a synthetic dataset p(5(0) within of the observed data in terms of the discrepancy. The ABC posterior is
Pe (6 x p(0) (Sobs 10),
(2)
and there are a variety of methods for sampling from equation 2 which are based on the fact that the right- hand side of equation 1 can be estimated unbiasedly by I(4o Sobe) < t) for S < p(S)0). The discussion above can be generalized by replacing the indicator function with a more general kernel.
2.2 Bayesian optimization
Next we briefly describe Bayesian optimization (BO), which is used in the BOLFI method discussed in subsection 2.3 A more detailed introduction can be found in Garnett (2023). Bayesian optimization is used for finding global optimum of function f()., 0 θ S R", where derivatives of f()) are not available and evaluations of f()) may be corrupted by noise. Here we consider minimization problems, but changing sign of the objective function turns minimization into maximization. BO models noisy evaluations of f(0) with a surrogate model deseribing our uncertainty about f(-) given the function evaluations made so far. This surrogate model usually chosen to be a Gaussian process (Rasmussen) 2003). and it guides the decision of where further function observations should be made.
Since the surrogate model we use for BO is Gaussian process (GP), we need some background about GPs. A random function f(-) defined on θ is GP with mean function μ ② + R and positive definite covariance function C θ x θ - I if, for any n, and any 0,....,,, Gn E O, the random vector f(11..)= = ((O)),...,,((O))) is multivariate normally distributed, with mean vector (u(0)),...,,u(0) u((n)), and covariance matrix C(f fmm) = [C(0,, 0 Suppose we observe the Gaussian process f(-) with noise at points O On E O. The noisy observations are
2 f(0,) + Ci, i= i=.,...,
where e, did (0, o²), for some variance 0² O. Write agn (2I,..., We are interested in describing uncertainty about f(0*) for some 0* O, given the noisy observations Zgn- The distribution of f(0*) [=<n is Gaussian, ((p, (0")))) where the form of pm (8*) and of (0*) are given in Section A.1 1 in the Appendix.
The uncertainty quantification provided by the Gaussian process surrogate can be used to decide which 0* should be used to obtain a further noisy observation
z" ==0 +C, C~N(000").
in our search for the minimizer of f(-). 0" is usually chosen to minimize so-called acquisition function. As an example, in the BOLFI method described next, Gutmann & Corander (0016) suggested using the lower
confidence bound acquisition function (Cox & John 1997 Srinivas et al. 2012
A,) (0) V (0), where 1/ = 2 log n F+2 m
3
with the default value G, = 0.1. Intuitively, (0) is an estimate of f(0) from the noisy observations so far, and choosing 0* to mininize (o) encourages choosing 0 where this estimate is small, or where (0) is large and we are highly uncertain about f(0). Hence the Gaussian process model can help to manage an "oxploation-cxxlootatio trade-off" in searching for the minimum. Other acqusstiin functions can also be used (Järvenpää al., 2019). Forms for the mean and covariance function need to be specified, and any parameters, including the noise o², estimated. This is usually done using marginal maximum likelihood.
4

---

Published in Transactions on Machine Learning Research (06/2025)
2.3 BOLFI
We now discuss the use of Bayesian optimization in the BOLFI method of Gulmann & Corander (0166) A discussion of more recent work improving the BOLFI method is given in Section A.2 of the Appendix. We can rewrite the ABC likelihood equation
as
(r (0) x (r(Ao(S, Sobs) t),
for S ~ p(510). Hence the ABC likelihood can be approximated if we know the distribution of (S, Sobs). ) S ~ 0) This suggests we may be able to approximate the distribution of (s Sobs) as a function of 0 using regression, in order to obtain an approximation of the ABC likelihood. This is what the BOLFI method does, while using a sequential design approach based on BO to choose which to simulate from next for maximum benefit. The sequential design aspect allows simulation efficient exploration of the high likelihood region, making BOLFI highly suited for the case of computationally demanding simulation models.
Gutmann & Corander (2016) propose selecting parameter values for simulation using a BO algorithm to minimize the expected discrepancy function:
D(0) - ((Ao(S,,So.))
(6)
where the expectation is taken with respect to S p(S(0). In their approach, first some initial set of locations 0g θ i I,..., Po, are chosen according to some space-filling design. Synthetic data are then simulated at these locations to obtain discrepancy values A; = Ao, (S, Sohn), S p(S)0.), = I,..., no- We can write
A = D(0.) + + i= 1, ,,...,,.0 n0.
(7)
where the 6 are mean error terms. In BOLFT it is usually assumed, perhaps after transformation of the A;, that E,
o²) for some variance parameter o².
Next, model D(0) as Gaussian process, and denote the training data used to fit the Gaussian process model required for BO as Tr, = ((,,, = 1....,,,,,pp nu]. A new location 0n is then chosen to simulate the next sunmary statistic value and obtain the discrepancy Ano+1- This is done by optimization the BO acquisition function, which uses the uncertainties of the Gaussian process model to define the benefit of simulating new discrepancy at any 0. The Gaussian process is then refitted with training data TRO+1 The (0 The process of optimization of the acquisition function, simulation and retraining is repeated until some computational budget n, >> no of simulations has been exhausted. The final Gaussian process model is then fitted to the training data Tzy. Using the Gaussian process model and the noise assumption for equation 7 we can approximate the distribution of a discrepancy value observed at any 0. we assume that A = D(A) + e, and write (0), o² (#) for the mean and variance of the predictive distribution of D(0) given Tny, then D(0) N (un, (0), o (0)) and
N0, independently. We can approximate the distribution of A as A ~ NHn, (0), o, (0) + 0²). Further details about the derivation and form of pny(0) (0) and o² 0) are given in Section A.1 in the Appendix. Employing this Gaussian approximation we obtain
to calculate the tail probability on the right-hand side of equation
1- 1-Mu, (0) 1022, (0) a2
Pr (Sobb (0) X +b
where & denotes "approximately proportional to". This likelihood approximation can be used with MCMC sampling to draw approximate posterior samples from (o)s)))). The choice of tolerance I in our examples is discussed later.
3 The MOBOLFI method
Next we describe LFI with multi-source data with a discrepancy for cach data source and the method of approximating the likelihood we consider. This is followed by background on multi-objective optimization, and finally the description of the new MOBOLFF method.
5

---

Published in Transactions on Machine Learning Research (06/2025)
3.1 Likelihood approximation with multiple discrepancies
To ease notation, consider the case of two data sources but the extension to three or more sources is immediate. Suppose the data comprises y = (x = E y = X' x W and decompose the joint donsity for y/ f as p(y(0) = p(z(0)p(wlzz 0). The observed data is Yobs = (T w p(yobal0) = p( (0)p
wobe) the likelihood is 8), and now there are summary statistics S = (T) .U )= : + R. where S concatenates les summary statistic mappings T and U for the data sources x and w respectively, T X R U R d b+c Write Sobs (yobs), Tobs (xobs) and Uobs (wobb)).
Similar to previous discussion of ABC, replace the likelihood p(yobs 10) with the summary statistic likelihood
p(Sohx = p(Tober Uabu = (TTotl0))p(L Tobs 0).
(8)
We then consider two discrepancies, (T, Toba) and Mo(U, Uobs), for simulated summary statistic values S = (TT, U and formulate an ABC likelihood approximating equation
Pe(Sobal0) x ( (A < (V), (Ubs) < dS,
(9)
where t = (21, is discrepancy vector with t, > 0 Another way of writing equation
(Sobbal)) x P(Aa(T, Toba) < t1, (U, Uoba) <t2)
= P(A <t) x PUVo (abu)) < (o, (T, Tobe) < tt),
(10)
for S = 3S=TTT,UT) p(S10).
10 approximates p( while the second approximates We could also approximate P(Uob=l0) and p(TobalUube 0) by switching the two discrepancies in equation 10 Our extension of BOLFI will model (s Sohs) for
The first term on the right- right-hand side of equation
p(Uobe|Tote: 0) (up to constants of proportionality)
p(S]0) as bivariate Gaussian process. We will use this bivariate process for sequential design using multi- objective Bayesian optimization, and also for approximation of data-source specific likelihood contributions such as those shown in equation I0
3.2 Multi-objective optimization
Our MOBOLFI exlension of BOLFI uses multi-objective optimization, and we describe this now. Let f(0) = (f. (()) be multivariate function, and suppose we wish to minimize the components of f()). There need not be any common 8" E @ for all components where a minimum is achioved, and multi- objective optimization methods approximate the set of "nondominated" solutions which are not obviously inferior any other solution. We say that a value with the inequality being strict for least one j. The dominated solution is inferior for minimizing f() along some dimensions and no better for other dimensions. Multi-objective optimization algorithms try to find the set of nondominated points in θ, the "Pareto optimal set.".
dominates 0' e if f,(0) S
K,
The Pareto optimal set is infinite in general, and numerical multi-objective optimization methods obtain finite approximations to it. The Pareto set is mapped by f(-) onto the Pareto frontier, the set of optimal function values obtained by the points in the Pareto set. Multi-objective Bayesian optimization (Garnett, 2023, Section 11.7) uses surrogate models to implement multi-objective optimization for expensive to eval- uate functions, possibly observed with noise. Similar Bayesian optimization with scalar objective, the representation of uncertainty given by the surrogate is used to efficiently decide where perform the next function evaluation.
The surrogate model in our work is a multivariate Gaussian process, and we need to explain what this means. Let f(0) = (x)
e O, be a multivariate random function. It a multivariate Gaussian process with mcan function / : - RRKK 100) = (0)....... (O)),,,, and positive definite covariance function C e x θ R R x if for any n, and 01.....0, An E θ, f(1..) (f(01) (.ee) f(.), is multivariate Gaussian with mean ((01..)) (u(O1)T ,-... ,,,,)) and covariance matrix C(01.n; 01.))
6

---

Published in Transactions on Machine Learning Research (06/2025)
())),,, with K x K block elements C(O,, 0)). Once again extending the discussion of subsection 2.2 suppose we observe values of [()) with noise at 0l 01,..., On O to obtain
2, = [(0.) +
(11)
where now 2; E RK and ej äid N(O, Σ) where E E RRKK is some positive definite covariance matrix. For any Hr E θ and writing 2gn (...,,121 2,TTT T. The distribution of f(g)) is multivariate Gaussian, (0)) (f")), where the form of (0") and (0") are given in Section A.3 in the Appendix. Given the uncertainty quantification provided by the multivariate Gaussian surrogate, an acquisition function can be defined. If there is a finite set of points, say . 8, approximating the Pareto set, with corresponding approximation fn of the Pareto frontier, one measure of performance that has been used is the volume of the space dominated by the current approximation of the Parcto fronticr and bounded below by a reference point, the so-called Pareto hypervolume. Expected hypervolume improvement (EHVT) was first used as an acquisition function in multi-objective Bayesian optimization by Emmerich (2005). For implementing the MOBÖLFI method of the next section, we use the noisy expected hypervolume improvement (NEHVI) method of Daulton et al. (2021) which copes well with noisy function evaluations and the trialing of batches of solutions in parallel with reduced computational demands. The method is implemented in the open source python package BoTorch (Balandat et al.] 2020). To the best of our knowledge, no asymptotic guarantees of recovery of the Parcto front have been proven for the NEHVI approach. [Daulton et al. (2021) derive a regret bound when trialing batches of samples chosen in a greedy fashion, and they also study theoretically the effects of approximating the acquisition function with sample average function approximation to it. Some random scalarization algorithms for multi-objective Bayesian optimization do come with theoretical guarantees (e.g. Zhang & Golovin 2020) and exploring such an approach in practice for LFI approximations with multi-source data is an interesting direction for future work. We will not discuss further the extensive literature on multi-objective Baycsian optimization, but refer the reader to Garnett 2023 Section 11.7) for an accessible introduction.
3.3 MOBOLFI
While BOLFI provides a closed form approximate likelihood for inference tasks, cannot approximate likelihoods for individual data sources for multi-source data. Multi-source data has become increasingly common in SSM design in recent years. Often it is natural to choose summary statistics and discrepancy metrics separately for each data source, and with the classic BOLFI approach we then need to combine all discrepancies into a single one. If this is not done carefully, the information in the combined discrepancy may not reflect correctly the relative importance of different data sources, resulting in information loss and highly conservative uncertainly quantification.
Motivated by these issues for multi-source data, we develop our MOBOLFI extension of the original BOLFI method. It achieves simulation efficient likelihood approximations LFI for multi-source data by applying multi-objective BO methods a vector of dala-source specific discrepancy functions. Consider again the bivariate setting and notation of subsection 3.1 for simplicity, and define the vector of expected discrepancies
D(0) E(Ao Soba)) = (D), (2)) T2)
where (0) (T, D2(0)= E( o (U, Uobs)), S (T,U) ~ p(S)0). A multi-objective Bayesian optimization algorithm applied to ((0) efficiently explores the set of 0 where both of the data source dis- crepancy components is likely to be small, leading to efficient approximations to the likelihood contributions from multiple data sources.
The optimization algorithm proceeds similarly to the case of a univariate objective. Firstly, we choose some initial set of points O,,..., On E θ according to some space filling design. We then simulate discrepancy values A (St, Sobs), St p(S)0,), i 1...., no, and we can write
A = D(,))) +
(12)
where the 6 are zero mean independent errors. It will be assumed that &, ~ N(O, E), for some covariance matrix E. Assuming a Gaussian process model for D(:), we fit Gaussian process surrogate model to the
7

---

Published in Transactions on Machine Learning Research (06/2025)
training data Tno {(0,, L,...,, nu). learning all hyperparameters including Σ from the data. We can then choose the next observation point 2 to minimize the NEHVI acquisition function (or perhaps choose a batch of points), retrain the GP and acquire new points, continuing until we have ny training points for the final fitted GP. In GP fitting, we use a constant mcan function and Matérn covariance kernel.
From the model equation 12 and the assumed (O,, Σ) distribution of the errors, we can approximate the ABC likelihood equation 10 up to proportionality constant by a bivariate Gaussian probability given well-chosen vector-valued tolerance = t2):
(S 10) = = ((1,, (2) Pay (0), (0) + E).
(13)
where ((: (0), En, (0)) denotes the edf of a normal NHH, (0), En, (0)+ distribution. Henceforth we will omit the "up to a proportionality constant" qualification when we talk about likelihood approximations. It is also possible to decompose the bivariate probability equation 13 into marginal and conditional components, approximating the terms (T, Tobs) < (1) and P(Vo(U, Uobs) < (2]
< 41) which in turn 8). The discrepancies can also be swapped in the above expressions.
approximate p(Tobal0 (0) and p(uobel Tobs.
Write Nhg (0) = (#m1 (0), (n)( ())) and write (0) and Egg for the j)th entries of ,, (0) and E respectively, j 1,2. We approximate pTooba by
41 Mmn (0) EEni (0) + E11
Pe (Tobbbl0)) = +b
(14)
and p(Unta Tobs, 8) by
Pe (Uobal Tobs: 0) = ( f/ - Mn: (0)
(15)
En221 (0)
where
En)2(0) (0) E12
(0) = (n2(()) + En ( + Ez2
(Tobs - Mun Un)(
and
10) 10) + (2n12(0) + E12)
Enzz(1) + Ez2
The likelihood for p(Untal0) can be approximated by
f2 pn (0)
Pe (Usbee =
(16)
(n22 (0) + 222
Gutmann & Corander (0166) chose the tolerance in the univariate BOLFI method the q-quantile of A1,..., Anys where y E 1). In the bivariate MOBOLFI method, we extend the choice of tolerance / = (f1,t2) to the 2-dimensional vector q-quantile of A,,.., Any, where E (0, 1)2. The approximate likelihood in oquation 13 is sensitive to the value t. For results in this paper, we set = 0.05. A comparison of different g-quantile tolerances is given for 3 different examples in the Appendix. The dependent noise covariance matrix Σ estimated by the covariance of simulated sample {Aig,))==, where Aiy.J = D(0rg) + 6ig are simulated noisy discrepancies for some 0gg. For results in this paper, we set ny 100 and 0:g arg arg min (A,- (0,) En, (6,)) '(A, (n,)))) For better performance and numerical stability, apply (9..A.)ET..
scaling to control the magnitude of A. Details of scaling provided in Section A.3 of the Appendix. We have described the MOBOLFI method in detail for the case of two data sources. While the extension to more than two data sources is conceptually simple, the multivariate Gaussian process surrogate caleulations in multi-objective BO can become more expensive when there are many discrepancies which need to be modelled jointly. With many data sources with corresponding discrepancies, and if model simulation is no longer the dominant computational expense for MOBOLFI, then BOLFI could be preferred if good way of combining the discrepancies into a single one can be found, and if we are confident there is no conflict between different data sources.
8

---

Published in Transactions on Machine Learning Research (06/2025)
Estimating covariance matrix between discrepancies in MOBOLFI in the noise model for the multivariate GP surrogate allows to capture dependence between discrepancies in flexible way. This could be especially important in situations where the data sources are dependent, as in our motivating SSM application with choice and response data. It is difficult to give a thcoretical analysis of the choice of distance in the BOLFI and MOBOLFI methods. However, BOLFI attempts to approximate an ABC method, and the case of ABC there is some relevant theory demonstrating the importance of scaling the ABC distance using the summary statistic covariance matrix. See Li & Fearnhead (2018) Section 3, Theorem 1) where it is demonstrated that without appropriate scaling, the limiting distribution of the ABC posterior mean may not coincide with the limiting distribution of the true posterior mean given the summary statistics, in cases where the number of summary statistics is larger than the number of parameters.
Estimating summary statistic covariance matrices can be difficult to do reliably when the number of summary statistics is large, since we need to simulate a large number of summarics at a point estimate of the parameter, which is undesirable in the setting of our work with a computationally expensive simulator. Our approach is to scale the summary statistics that they have roughly equal variance, and to caplure the dependence between the summary statistics through the covariance matrix of the vector of discrepancies, rather than the vector of summary statistics. The covariance matrix of discrepancies is low-dimensional, making estimation easier for a computationally expensive simulator. If we combined information from vectors of summary statistics by a linear weighting of their corresponding discrepancies to obtain a scalar joint discrepancy, this would not be a flexible approach to allowing for the dependence between discrepancies. For example, with k discrepancies we have & linear weights, but the covariance matrix of the discrepancies used in the MOBOLFI surrogate has k(k 1)/2 distinct parameters and provides greater flexibility and many transformations can be used in the noise model. Furthermore, it is not obvious how to estimate good lincar weights for combining discrepancies, whereas estimating the covariance matrix of discrepancies in MOBOLFI by simulation al a point estimate is relatively simple.
The above discussion is for the case where the model is correctly specified. In the case where the model is misspecified and indivdual data source likelihoods are peaked in different parts the parameter space, then BOLFI with a single discrepancy formed as a weighted sum of data source discrepancies will not explore the regions of high likelihood for all the data source likelihoods in single run. Hence single joint discrepancy is undesirable in that setting, matter how it might be chosen, which is strong motivation for the multi-objective approach. This is demonstrated in an example in Section 4.
3.4 Checking consistency of different data sources
There are number of advantages in separately approximating likelihood contributions in problem with multi-source data. One of the most important is the ability to detect conflicting information about the parameter from different parts of the data. Given a prior p(0) we can use the likelihood approximation equation 13 to sample from the approximate posterior
(o) & (So)):(S.tn.0))
(17)
Using MCMC to sample from this does not involve any further (computationally expensive) simulation from the model. This is true for the other posterior approximations discussed below also. If the two components of the likelihood are in conflict and induce modes in different regions of the parameter space, our method of multivariate approximation can capture some of that complexity by separately considering the likelihood contributions for the different data sources.
we want to know what information is contained in the first data source only, we can compute the approximate posterior
(B x (Tosp)(Taab))))
(18)
Comparing the posterior densities equation 17 and equation I8 tells us how the second data source changes the inference. We could also consider a weakly informative prior pw (0) and consider the information about the parameter contained in each data source through the posterior approximations
Pe (0 x (O))oe
(19)
9

---

Published in Transactions on Machine Learning Research (06/2025)
and
(ob(b)))) x pw (0)pe(Uobsl)))
(20)
The purpose of using a weakly informative prior here is to separate the information in the data from that contained in the prior, in so lar as that is possible.
In this work we consider only informal comparisons of posteriors based on different data and prior choices, but there are more formal methods check for conflict between priors and data, more generally between
information in different parts of hierarchical model. We refer the reader to Evans & Moshonov (2006)
Marshall & Spiegelhalter (0077) and Presanis et al. (2013) for further discussion of these.
4 Examples
In this section, we implement MOBOLFI in several examples. First, we apply MOBOLFI a simple synthetic data example where can illustrate the advantages of MOBOLFI compared to BOLFI for multi- source data, including the situation where misspecification is present. Following this, use MOBOLFI to estimate soquential sampling model (SSM), the Multi-attribute Linear Ballistic Accumulator (MLBA), to highlight MOBOLFT's superior design for choice-RII joint data sources compared to BOLFI. MLBA is chosen from other SSMs because of its generality for diverse choice situations and its closed-form likelihood, which enables gold-standard comparisons to be made for the likelihood approximations used in MOBOLFI. MOBOLFI requires training a bivariate GP, which involves higher computational cost than the BOLFI training of univariate GP. However, if model simulation is computationally expensive, this will dominate the computation time for both methods. Although the likelihood is tractable for the MLBA model, summary statistic based Likelihood- Free inference can be of interest when the assumed model is misspecified. In this case, conditioning on insufficient statistics that discard information which cannot be matched by the assumed model while matching important features can be used to develop models that are "fit for purpose" see for example Lewis et al. (2021) for further discussion.
The detailed setup of our experiments and further analyses are described in Sections Section D in the Appendix discusses an example on bac-
B. and C of the Appendix. terial transmission in day care centres (Gutmann & Corander 2016 Numminen et al. 2013
Multi Code and detailed results are submitted to github: https: /ttps:///tttttt..t htt/cibbbcccm//0ss//
Likelihood Free-
4.1
Toy example
Our first synthetic example uses the MOBOLFI method to infer a parameter 0 = (01....,010) 010) E R for model with two data sources.
The example is modified from Schmitt et al. (2023 The first data source consists of N 20 independent The second data source .....,50. 50. The Wm are obtained by observing
10-dimensional normal observations with mean 8, X,, ~ N(0, I), consists of M = 50 10-dimensional observations Wm, m = 1,..., discretely the following ten-dimensional Brownian motion (BM) with drift:
N. 1....,N..
dw(t) = Adt + ot)( E { [0,3].
where wo(0) is a ten-dimensional vector of zeros, W(!) is a standard ten-dimensional BM, and 0.5, leading to Wm w(mm 1)0), ó 3/(M - 1), M..
For the dataset X = (X,, N-1 Σ., Xn, X = N. X and the discrepancy for data source X is (X, X²) = I X'1, where || denotes the Euclidean distance. For data W = {Wm) IM and observed data W = {Wo JM
and corresponding observed data denoted {X #n we write X =
m==
we write AW,,. Wm+1 - m, m I,...,AA . AWo m W W m 1,.... M 1, AW
(M 10--1)7 11)- M AW and AW = (M-1)-1 > m= AW Then the discrepancy used for data source Y is 42 W AW A.
10

---

Published in Transactions on Machine Learning Research (06/2025)
1.2
MOROLFI BOLF prior
1.0
XOBOLFIX DOLFI
1.0
NOBOLFI W
BOLFI
peior
prior
0.8
0.8
0.8
Brue
Piue
Dmnyy
0.6
Demty
0.6
0.6
0.2
0.2
0.2
0.0
0.0
0.0
0
2
0
2
0,
0,
(a) Joint Likelihood
(b) Information in X
(c) Information in W
Figure 1: Approximate posteriors for the toy example. The left column shows approximate postcriors for the joint likelihood. The middle/right columns show approximate posteriors conditioning only on X/W. The blue and orange curves are kernel density estimates obtained from HMC samples for MOBOLFI and BOLFI posteriors respectively. The green curves shows the true posteriors (joint/X/W X/W on left/ heft/middlo//ighhee column respectively), and the red dash line is 0true -077.
We generate the observed data (X", Wa) for the analysis by simulation with Øtre (-0.7, 0.7, -0.7, For implementing BOLFI and MOBOLFI, 100 initial prior samples (00), (x((, ((),)))))0 are drawn, where x = (x!1) In=1: W( = (W Im=l: 0(0) ~ (0, I), and X( and W(O are the ith simulations for and W given O6). The training data is (A)(())))))) X²), A2((0,, Wo)) For implementing the BOLFI approach, we get a single discrep- ancy by weighted sum of the two data. source specific discrepancies, (X, Xº) + A2WW, Wo). The weight of 0.4 on the first discrepancy is chosen to put the two data source specific discrepancies on similar scale. For implementing MOBOLFI, we use the vector discrepancy (A1, A2). A2). For both approaches, 150 acquisitions are made in the BO algorithm, and 250 model simulations are needed in total.
Given the symmetry of the model in the way that the components of 8, X and W are generated, without loss of generality we present only results for inference about 0₁. Figure [ compares MOBOLFT and BOLFT approximate posteriors of 0j. We make three observations. First, the true posterior has smaller variance than the approximate posteriors. This is expected and mostly reflects the finite tolerance used and uncortainty in the likelihood approximation from the surrogate model. Secibd, the MOBOLFI approximations to the posterior conditional on X or W only show similar posterior location to the corresponding true posteriors, but with somewhat larger variance. The MOBOLFI posteriors conditional on individual data sources are obtained without significant additional computation after the posterior approximation for the joint posterior has been obtained. Third, in the left column of the figure, Hamiltonian Monte Carlo (HMC) samples from the MOBOLFI approximate posterior exhibit lower variance compared to those from the BOLFI postcrior, and the MOBOLFI posterior is closer to the true posterior. This suggests that there is information loss from combining the discrepancies from different data sources compared to MOBOLFI with multiple discrepancies unless the combined discrepancy is constructed very carefully. Section A.3 discusses the method used to scale discrepancies in the BOLFI and MOBOLFI algorithms. Figure 2 shows the effects of changing the scaling and for BOLFI (left) there is greater sensitivity than for MOBOLFT (right). In combining discrepancies for BOLFI with multi-source data, a linear combination of the component discrepancies may be sub-optimal, and obtaining a good combinination may not be as simple as choosing a lincar weight.
Additional findings for this example are described in Section B.2 of the Appendix, where we explore the effect of tolerance and the number of BO acquisitions on the inference. Another experiment described there investigates the performance of MOBOLFI when parameters are present in only subset of the data source-specific models.
4.1.1 Misspecified simulator
We continue the above example by modifying to introduce misspecification. Our purpose demonstrate that MOBOLFI is useful for understanding misspecification when different data sources pro-
11

---

Published in Transactions on Machine Learning Research (06/2025)
1.2
MOBOLFI scale 0.1
0.6
BOLFI true
1.0
MOBOLFI scale 0.7
MOBOLFI scale 1.0
MOBOLFI scale 0.4
0.8
0.4
Denstt
0.6
0.4
0.2
0.2
0.0
0.0
-2.0
-1.5
-1.0
0.5
0
0
(a) BOLFI with different scaling
(b) MOBOLFI with different scaling
Figure 2: Approximate posterior in the toy example given different scaling weights to A1. The left column shows the BOLFI approxinate posteriors given weight = 0.1 0.1/0.4/0.77 /0.4/007/1.. scaling, while the training objec- tives for BO is defined to be weight * Å + A2. The right column presents the MOBOLFI approximate posteriors given the same scaling.
vide conflicting information. In the modified example the observations are one-dimensional, Xn ~ No0 1), n 1, .N, and Wm, m I,...,M,, M, with Wm w((m 1)6), 8 3/(M 1) with w(t), E [0,3] a univariate process defined by (w( = Odt + (t), w(0) = 0,
with (() a univariate standard BM, and N, M and as before.
The true data gencrating process (DGP) is not the above assumed model. Instead, in the true DGP there are different values of in the models for X and W, so that X, ~ (0x,, 1) and dw<(t) = dt (t), where 0x # lw. We simulate the observed data for the analysis using 0x = 0.3, Ow = -0.7. In the misspecified assumed model the likelihood contributions for different data sources produce conflicting information about 0.
Approximate posteriors for are presented in Figure 3 The MOBOLFI and BOLFT approximate posteriors are similar when conditioning on both X and W, with the posterior mode lying between 0x and Ow- The middle/right columns shows the MOBOLFI approximate posterior conditioning on only X, or only W, compared to the corresponding true posterior and the BOLFI posterior conditioning on both and W. We make two observations. First, the MOBOFLI approximate postcriors conditioning on X only and on W only produce good estimates of 0x and Ow, and the conflicting information in the two data sources is evident. Second, the approximations of the individual data source posteriors in MOBOLFI are obtained without substantial additional computation, since the joint likelihood and likelihoods for X and Y are obtained simultancously in the MOBOLFI approach.
4.2 Multi-attribute Linear Ballistic Accumulator
MLBA is a state-of-the-art SSM for understanding the human decision decision-making process. We will apply MOBOLFI to infer parameters of one variant of MLBA (Hancock et al. (021b) using a simulated datasct where 0 is known, and for real-world dataset on preferences of ridehailing drivers to rent electric vehicles in Singapore.
4.2.1 Simulator introduction
MLBA assumes that decision-making can be seen as an evidence accumulation process, illustrated in For a prcision-maker facing M alternatives, MLBA posits there exist M independent evidence accumulators starting to collect supporting evidence for each corresponding alternative a from a random indep. starting point Ga i.id. U(00 A). Each alternativo-ppccifie accumulator evolves linearly drift rate v N(d,, s²), where 1, ... d, depends on the pairwise attribute comparisons for all K attributes of alternative with the remaining alternatives for each observation. A and x² are constant over observations to
Figure 4
12

---

Published in Transactions on Machine Learning Research (06/2025)
1.0
NOBOUPI BOLFI
1.0
MOBOLFI
1.0
MOBOLF W
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.2
0.2
0.0
0.0
0.0
(a) Joint Likelihood (mis)
(b) Information X (mis)
(c) Information in W (mis)
Figure 3: Approximate posteriors for the toy example under misspecification. The left column shows MOBOLFI and BOLFI approximate posterior densities conditional on both X and W. The middle/rige presents approxinate posterior densities conditioning only on X/W. The blue and orange curves are kernel estimates of posterior densities from HMC samples for MOBOLFI and BOLFI respectively. On the left column, the green curve shows the prior while the red curve shows the true posterior; on the middle/right columns, the green curve shows the posterior conditional on X/W while the red curve shows the prior. The two dash lines are values of parameter produced by data sources with conflict information, labeled as true and mis(specifiedl) 0x 0.3 and Ow -077.
mcasure the scale of initial preference and unobscrved accumulation process noise, respectively. Finally, the alternative whose evidence accumulator reaches the common threshold x first is considered the final choice, and the deliberation time is the same this alternative's evidence accumulation time x=3, Strictly, the response time (RT) is defined as the sum of the deliberation time and the non-decision time, which stands for information encoding and decision execution time. However, for complex choice tasks with multi-attribute information, non-decision time can be omitted if it does not contribute significantly to RT.
Fridence conmulation Amount
choscn. Fme
Decisian Thrmdoll, z
Relouse (8L)
(imt sufficent time P(ATIRC -
#
Starting point upper bound,
Deift rate uy. Níd,, s) independently.
Starting peints, A) L.d.
Non Time (c)) Desiion Time (DT)
The simulation process of MLBA for three alternatives (M = 3). In this case, the final choice outcome is Alternative (in blue and accumulator is bold). The RT is the time for Alternative 1 reach the threshold with non-decision time To. The directed solid lines demonstrate the realized evidence accumulation trajectory for each alternative. The ML.BA simulator contains three steps annotated above. After iterating sufficient times, the choice proportion and its conditional RT distribution can be acquired.
Figure 4:
13

---

Published in Transactions on Machine Learning Research (06/2025)
The MLBA model can be described in two parts. The back-end part, Linear Ballistic Accumulator (LBA), was originally proposed by Brown & Heathoote (2008) for implicit choice situations (i.e., choice tasks without listed attribute values). To extend LBA to MLBA, Trueblood et al. (2014) added the front-end part to LBA by specifying alternative-specifee drift rates. Several variants of the frout-end part have been investigated to define the drift rate mean of MLBA. We choose HHaccock et al. (2021b)'s 's specification because it can handle more than two attributes. For a choice set C, the drift rate mean da for the alternative a is:
de max - 8, +
Σ
Wask BabkBB(AK (R -6b),
(21)
bEC,byaak 1=1
where Xak is the value of altribute k of alternative a and
-1 (Xas - -xxx]3 (x Xok) 2 0
Wabk
(22)
exp(->2 3a Xbx)1} Bx Xna) < O.
The final simulation output is a concatenation of RT and choice outcome:
(RT, CH) (((x (cn)/eem) + To. argmin(( - gn)/
(23)
aEC
: (A1,A2) is interpreted a vector of sensitivity parameters of positive/ aositive/negative negative difference of attribute pairwise-comparison respectively; 1o is a common drift rate mean constant that prevents negative drift rate; g = ($1....9k] are scaling parameters for the attributes; and 5 (f..... 5M) are alternative-speeific constants for the drift rate mean. Thus, the estimable parameter is A (A. To, 3, 1-11 x), while (A, To,61) 61) are fixed for parameter identification and 5_1 = (82, ..., ".. Following Brown Hcathcotc (2008) and Terry et al. (2015), the joint likelihood of observing the response time RT and the choice made CHH of multiple candidates has closed form expression. It is presented in Section C.1 of the Appendix. We use the closed-form likelihood as gold-standard to benchmark approximate posteriors obtained by LFI in subsequent experiments.
4.2.2 Synthetic data
We use the survey question design and alternative-specifi attribute matrix for the real example of soction 4.2.4 A choice question with M = 3 alternatives is provided to candidates. Each alternative offers values of K = 3 attributes, and a respondent will make decisions by evaluating the alternative and its attributes. For data generation, the synthetic output data is simulated using 0 = (A, To, 6, -5, -6), (0, 1.5),, resulting in 1280 = 320- (3+1) 1) observations, after concatenating CH (in ono-hot encoding) and RT for 320 candidates. The size of the alternative-speeific attribute matrix ftuue
(.,,,,
X is 320 9. For parameter estimation, the true values of parameters we larget. in inference (A1, B2, 82, dz, log(x A))true (0.1, 3,,,,, 1.5, 19g(9))). The parameters A, To, dr are fixed 1, 1, O, and 0 for identifiability, and A2 = 0.8, lo 2, B3 -6 are assumed to be known for reducing the computational cost in further training such as facilitating scaling of the estimable parameters. This implementation is common in empirical studies since MLBA is designed bottom-up so that some parameters are highly correlated with each other. Lastly, we have log-transformed x A in order to increase sampling efficiency, since the magnitude of log(x- - A) is much smaller than that of X-
The synthetic data RRTOCCHP is simulated by MLBA taking gtrue input, where RT and CH are response time and choices. For implementing BOLFI and MOBOLFI, 100 samples from the prior are used as initialization, (0(9), (RT(), (1 where RT()) (RT(() 220)), and CHÚ (CH[",. ...., CH( are the simulated data given g00, For response times RT, write RT for the vector of order statistics of RT. Similarly, RT is the vector of order statistics of RT The training data is (A, RT)), (zCH(() )' where (RT, RTT) and 42(CH, CH²) are discrepancies for the response and choice data respectively. The first discrepancy is defined as
(RT RTT) (og((RT)) log( (RT)
(24)
14

---

Published in Transactions on Machine Learning Research (06/2025)
Closed form
Closed fomm NOSOLLF NOBOLT
NOROLFI MOROUII
1
0
-0.5
0.0
0.5
1.0
1.5
2
A1
logix A)
1.0
Closed
2.0
Cloeed
MOBOLI
BOUF
1.5
BOLTI
0.6
0.4
0.5
0.2
0.0
0.0
-25
-5
B2
Clised-fm
MOROLF MOBOLFI ROLF miee
0.6
MOBOLTI
ROLR
Donate
Deviity
0.4
0.4
0.2
0.2
0.0
0.0
62
8g
Figure 5: Approximate posteriors for MLBA example. The plots shows marginal posteriors for different methods and each parameter of interest. The green/purple curve represents the MOBOLFI approximate posterior/BO1 BOLFI approximate posterior. The red curve is the sample from the closed form posteriou caleulated using the closed form likelihood of MLBA. The orange/bil /blue curves are the MOBOLFT approximate posteriors calculated by marginal approximate likelihood of RT /CH data. The red dashed line is the value of girue
where |I denotes the Tn distance. For the second discrepancy, recall that the CH data represented by dimensional one-hot encoding of 3 alternatives, and we define
320
Ae CH*) 3 320
Σ ICH -C 1=1
(25)
The prior assumess independence between components of 0, with each marginal prior uniform on an interval of length 8 centred the true value, with the exception of the parameter À which has a prior uniform on [0, 1]. Differential emission MCMC (De-MCMC, see Turner et al. (2013) is used to sample from the approximate posterior. Further details of the experiment are given in Section C.2 of the Appendix.
Figure 5) shows approximate posterior distributions obtained using the BOLFI and MOBOLFI methods. We make number of observations. First, the variance of the MOBOLFI approximate posterior is smaller than that for BOLFI. In fact, MOBOLFI using only the response time data produces more accurate approx- imations of the true posterior than BOLFT given both data sources. This suggests that when it is natural to
15

---

Published in Transactions on Machine Learning Research (06/2025)
define discrepancies separately for different data sources, combining them linearly, even with weights, may result in information loss. Second, both MOBOLFI and BOLFI posterior variances are larger than for the true posterior, which is expected, due to the irreducible information loss from approximation. Third, when comparing the MOBOLFI posteriors calculated by likelihoods of different data sources, the posterior using both data sources is better than any posterior using only one data source. Fourth, for inference / and A1, the posterior using only CH has its marginal maximum a posteriori probability (MAP) value far from the value of (true; for inference of X the posterior using only CH has variance larger than that of the posterior using the joint likelihood, compared to using RT only or the joint likelihood. In psychology and economics, researchers usually focus on the choice data CH in parameter inference for SSMs like MLBA. The figure shows that using response time data in inference not only reduces the approximate posterior variance to make it closer to the true posterior variance, but also helps locate the area of the global maximum. Such conclusions are likely to extend to other variants of MLBA and LBA. Section C.3 of the Appendix details additional experiments on the factors affecting the performance of MOBOLFI in this example.
4.2.3 Misspecified synthetic data
We now extend the previous synthetic MLBA example to a misspecified scenario, illustrating the advantages of MOBOLFI for efficiently exploring the high likelihood region for both data source likeli- hoods when there is conflict. We define two sets of paramcters @ame (0.2, 25, -3., 6, 4, log(99) and Otrue (0.05, 24, -6.5, 3, 1.5, log(199)) for simulating response time and choice data respectively. We sim- ulate data (RT*, CHr1) by MLBA taking Olrue as input, and data (RTch, CHoh)) by MLBA taking fur input. The synthetic observed data (RT",CH*) = (RTT,, CHHh)) concatenates the response time data generated using R and the choice data simulated using Orr.. Therefore, the two data sources provide conflicting information about the parameter vector.
0.00
-0.05
-010
pooossal
-0.10
0.15
-0.15
init
BouFi
0.20
-0.20
BOLFI ch
025
-0.25
-0.55
-0.35
-040
-0.4
-0.2
Figure 6: Noisy Negative discrepancy values of Bayesian Optimization acquisition points. The horizontal axis denotes the negative RT discrepancy, the vertical axis denotes the negative CH discrepancy. The discrepancies have been scaled to a similar range. The yellow/orange point represents simulated discrepancies obtained using Oize// /0rzue as input. The green point is simulated vector of negative discrepancies obtained at the mode of the MOBOLFI approximate posterior.
We follow the same design as in See 4.2.2 with some minor changes. Firstly, we remove the rejection criterion for sampling the initial training set for BO discussed in Section C.2. With conflicting information, parameter samples are likely to obtain a low value discrepancy on one data source and a high value on the other. Secondly, the number of iterations for the BO algorithm set to 500.
Applying MOBOLFI such (RTO, CH*), we find that the multi-objective Bayesian Optimization explores range of local optimum points near both Oirue and gtrue. Although we cannot observe the expected discrep- ancy, Figure 6 (left) shows the noisy negative discrepancy values of MOBOLFI acquisitions. The plot shows that the algorithm explores a range of parameter values lcading to small discrepancy for one discrepancy or the other, even if both are not small simultaneously due to the conflict. Dots on the plot show simulated discrepancies for gtrue and orrue. The plot on the right in the figure shows simulated discprepancies (for both discrepancics) for BOLFI acquistions based on one data source. Green dots are for BOLFI acquistions based on choice data discrepancy, and dark blue dots are for BOLFI acquistions based on response data discrepancy. Due to the conflict we that the BOLFI acquisitions cannot explore in single run the
16

---

Published in Transactions on Machine Learning Research (06/2025)
full region where either one of the discrepancies is small. Hence cannot obtain reliable estimates of the likelihood for both data sources based on BOLFI acquistions for one data source only. A similar problem would arise for any composite discrepancy linearly combining the data-source specific discrepancies. Plots of approximate marginal posteriors are shown in Section C.3 of the appendices.
4.2.4 Real-world data
Next apply MOBOLFI and BOLFI to infer parameters of an MLBA model for real-world dataset from a consumer choice experiment. This experiment assesses the rental preferences of ride-hailing drivers in Singapore for electric vehicles (EVs) via a street-intercept survey (Ding et al. 2024 Before the stated preference experiment, the driver's basic information such as working days per week, and information about currently rented internal combustion engine vehicles (ICEVs) is collected. The drivers were asked to make a choice among three alternatives including ICEVs, Electric Vehicle Model A (EVA), and Electric Vehicle Model B (EVB) with three listed attributes, which are monthly rental cost (RC, in SGD), daily operating cost (OC, in SGD), and daily mileage (DR, in km). The time that elapses from the appearance of information to the confirmation of choice is recorded as RT. Note that EVA and EVB are assumed to be identical except for their values of three listed attributes. After data preprocessing, 149 participants with 584 valid observations are used in the paramcter ostimation. Since monthly rental cost (RC) and daily operating cost (OC) capture the monetary aspect, these two can be merged into one alternative specific altribule for the monthly total cost (TC, in SGD).
52
TCna RC + OCna x WF X 12
(26)
where n and a are the indices of observation and alternative respectively. WF, is the number of working days per week, hence TCno is the monthly total cost for alternative a in observation n. Lastly, due to the hcavy-tailed distribution of DR, it is transformed to log scale. In summary, the size of the attribute matrix X is 584 x 6, and the real output data for the MLBA simulator (ITT, CH) is 584 x 4.
The unknown parameter is 0 = (A1, STC, log Slog DR, SICEV, log (x AA), and the rest of the parame- ters are fixed at known valuos. We set constants A = -5, s = 1, 1,8 p = 0, and To = 1 to cnsure parameter identification. The allernative--peecfii constant ō is often interpreted as the initial preference of the cor- responding alternative before considering the attribute values in empirical studies. A mild and reasonable assumption made here is that Ō δεVB 0, since both are EVs. Furthermore, based on prospect theory that consumers pay more attention to loss than gain [Canerer] 2000), X is normally smaller than A1. Therefore, to simplify the estimation, we set A2 = 0 for the final empirical analysis. Finally, similarly to synthetic data estimation, STC, ADR and X A are log-transformed to increase sampling efficiency.
We use independent priors d1 ~ u[0,5], 5]. SICEV ~ L[-3,55 5), log 3TC ~ 2[- 5], log SDR ~ 31-3, 5]. and log (x- - A) UI 6]. For the implementation of MOBOLFI and BOLFI, 282 initial samples (06)) (KTO,) ( 282 are obtained using a space-filling design (Santner et al., 2018). The training data where (RT CH²) are observed empirical data. Unlike the synthetic data experiment where single simulated dataset is used for computing discrepancies at each parameter by equation 24 and oquation 25 (RT(,, CH()) = {(RT!!) (H.") is a colloction of multiple simulations independently generated by the MLBA simulator with 000 and attribute matrix X. The dis- crepancy for the RT data source (RT(()) RTT) is the log-average of 1-Wasserstein distance between the empirical distribution of log- transformedd RT and the = 50 replicates of (T.
is (0()) (RT() RT"), ,
(R RTO) log ES ||
log(RTT) log (RTT. M11
(27)
35
where RT is the vector of order statistics of RT.
The discrepancy AH(HH), CHo) for the choice data is the log-average of SSE of choice proportion over the S replicates for the synthetic choice data and the observed choice outcomes (i.e. the average squared
17

---

Published in Transactions on Machine Learning Research (06/2025)
difference between empirical choice proportions by simulations and observed outcomes):
S I2CHH! - -CH* CH
AHCH = log
(28)
9SN2
where NN- = 584 is the number of observations and IN E R is the vector whose each element is 1. In addition, BOLFT's one-dimensional discrepancy for joint data sources is the summation of AC(H )ad)and A5(RT(), RTT) in equation 28 and equation 27 Therefore, the training data for BOLFI is (069, (AI(RT()) RTT) + AHCHH( (H
2822
There are two reasons for adopting A instead of A in equation 24 and equation 25 for the empirical experiment. First, log transformation on both discrepancies encourages the MOBOLFI and BOLFI to explore more in the high-density area (when A is close to 0), and therefore, letting the approximate log-likelihood be more sensitive (larger gradient) around the high-density area in parameter space. Second, the simulation data from S replicates provide a more robust discrepancy for each data source by reducing randomness brought by the MLBA simulator. De-MCMC is used for the MOBOLFI/BOLFI approximate posterior and closed-form MLBA posterior estimation.
Figure shows approximate marginal MOBOLFI and BOLFI posteriors for selected parameters and the corresponding posteriors for the closed-form MLBA likelihood. We make two observations. First, the MOBOLFI approximate posteriors (green curves) are much closer to the closed-form posteriors (red curves) than those of BOLFI (purple curves) for all parameters with the same training iterations, highlighting MOBOLFT's advantages for multi-source data. Second, the figure demonstrates that the contributions of different data sources to the estimation of individual parameters can vary greatly according to the parameter, and MOBOLFI can reveal these relationships. The CHH data contributod more to the estimation of the front- end parameters of MLBA (ōrcev, A1) related to the preference across alternatives and attributes, while the RT data contributes more to estimation of the back-end parameters of MLBA (log (x A)). For example, the MOBOLFI marginal posterior of λ and log (BDR) is very diffuse when conditioning on only response- time data, showing that the MOBOLFI posterior conditioning on both data sources is informative mostly because of the choice data. In contrast, the log (x - AA) MOBOLFT marginal posterior conditioning on both choice and response time is more similar to its marginal posterior conditional on response time only. This is consistent with the interpretation that X reflects the decision difficulty and respondent's cautiousness, which are associated with RT. Lastly, the approximate MOBOLFI posterior conditional on both data sources is always less dispersed than the MOBOLFI posterior conditional on a single data source only. This observation is aligned with the finding (Li & Bansal 2024) that joint choice-RT data provides more elficient parameter estimate compared to an estimate applying marginal data sources (CH or RT only).
The posterior estimates of MOBOLFI are different from those from the closed-form likelihood in Figure 7 mainly on log STC and ÒiCEV. There are several sources of error arising in the application of the MOBOLFI method. One the choice of discrepancy, which part of the definition of the BOLFI or MOBOLFT approximate likelihoods, which differ from the true likelihood. A second source of error is the choice of tolerance parameter in forming the approximation. A third source of approximation is that the approximate likelihood is calculated from the BO Gaussian process surrogate. When the number of model simulations is small, the estimate of the expected discrepancy provided by the surrogate and estimates of its noise parameters may be poor, and the approximate likelihood obtained from the surrogate may be inaccurate even a good discrepancy has been used. This source of error will be reduced as the number of model simulations increases. A final major source of error can occur if the error model in the BO surrogate is inappropriate this can be mitigated by more sophisticated error models, or transformations. We have taken care in implementing the BOLFI and MOBOLFI approaches to use sufficient simulations, to choose the tolerance carefully and to ensure the adequacy of the error model. In this example find that the most important factor in the performance of BOLFI and MOBOLFI is the choice of discrepancy. Figure g examines how well the observed data are reproduced in the fitted model for point estimales from the posterior means and MAP estimates for the closed-from likelihood and MOBOLFI approximate posteriors. The left column in the figure shows that log discropancies for synthetic data simulated using the posterior mcan of MOBOLFI are generally larger than those for the closed-form likelihood posterior mean value on both CH and RT. The right column of the figure compares simulated discrepancies for MOBOLFF and close-form
18

---

Published in Transactions on Machine Learning Research (06/2025)
5
0.6
Densty
Densty
3
0.2
1
0.0
2
-2.5
A1
ÒICEV
2.5
1.5
2.0
Densty
1.0
Density
1.5
1.0
0.5
0.5
2
5.0
2.5
0.0
2.5
5.0
7.5
(BT()Tc)
(BgR))))
Closed-form
MOBOLFI with
1.5
MOBOLFI ch with A"
Density
MOBOLFI A"
BOLFI with A*
MOBOLFI mean with A
0.5
Closed-form mean
BOLFI with
0.0
-5
10
log(x- A)
Figure 7: Approximate and closcd-form posteriors for MLBA example. The plots show estimated marginal posteriors for different parameters. All posteriors are istimdded using De-MCMC sampling and kernel density estimation. The coloured lines and dots show the results for different methods. The methods compared are estimates obtained with closed-form likelihood (red), MOBOLFT with both data sources (green), MOBOLFI with only choice data (bluc), MOBOLFI with only response data (orange), and BOLFI with both data sources (purple).
MAP point estimates. The RT log discrepancies simulated using MOBOLFI MAP smaller, while CH log discrepancies simulated using closed-form MAP are smaller.
5
Discussion
We developed Multi-objective Bayesian Optimization for Likelihood Inference (MOBOLFI) to address Bayesian inference with multi-source data in complex models with intractable likelihood, such SSMs. MOBOLFI extends the classic BOLFI method, which acquires model simulations at the most beneficial parameter values by using Bayesian optimization (BO) applied to mininization of an expected discrepancy between synthetic and observed data. The surrogate model used in BOLFI is also used for approximating the likelihood. The MOBOLFI extension for multi-source data considers a discrepancy for each data source, and considers multi-objective BO for exploring the parameter values where any individual data source likelihood is high in simulation efficient manner. Major advantages of the approach include avoiding information
19

---

Published in Transactions on Machine Learning Research (06/2025)
closed-form mean
closed-form MAP
100
MOBOLFI mean
MOBOLFI MAP
80
frequency
frequency
60
40
20
20
-3.0
-2.8
-2.6
-2.4
3.4
3.3
3.2
3.1
3.0
2.9
RT discrepancy
RT discrepancy
120
closed form mean MOBOLL mean
120
closed-fone MAP
100
100
frequency
80
frequency
80
60
60
40
40
20
20
4.4
4.2
3.88
-3.6
-3.4
4.44
-4.3
4.2
-4.1 CH discrepancy
-3,9
CH discrepancy
Figure 8: The histogram of discrepancies with 1000 replicates based on different point estimates. The left column denotes the RT (top) and CH (bottom) discrepancy distributions for posterior means of MOBOLFI (green) and closed-form likelihood (red), respectively. The right column shows the analogous RT and CH discrepancy distributions of the posterior MAP based on MOBOLFI and closed-form.
loss from naive methods for combination of data-source specifie discrepancies into single discrepancy, and the ability to approximate both the joint likelihood and likelihood for individual data sources. The latter is useful for detecting conflict between data sources, and for understanding the importance of the data sources for estimation of individual parameters.
Although the MOBOLFI mothod is motivated by the parameter inference SSMs with multi-source data, can be applied Likelihood inference of other complex models. For example, Section D of the Appendix we consider an example on bacterial infection in day care centres where parameters of latent Markov process are inferred. We partition data summarics into two and apply MOBOLFI for inference. MOBOLFI is competitive with BOLFI even in this example with only an individual data source. For the case of SSMs, we use the ML.BA simulator as the simulator of interest in our experiments, and MOBOLFI could be used in a wide range of SSM variants, especially for those whose likelihood functions are intractable like DDM, MDFT, etc. An interesting direction for future work to consider different choices of the acquisition function in the MOBOLFT approach. Multi-objective BO is an active area of research Bayesian optimization, and new developments BO methodology may translate into improved performance in Likelihood-Pres applications. Stale-of-the-art computational methods for the MLBA with closed-form likelihood, including variants of the model including random effects, discussed in Gunawan ct al. 2020). Evans 2019) addresses the important issue of efficient simulation from various SSM models, which is particularly important for likelihood-free inference.
In applying MOBOLFI there are a number of challenges. One that the initialization the algorithm can be very important, and an adoquate initial covering of the space needed. This seems be especially crucial with acquisition functions based on expected hypervolume improvement, where the BO algorithm keep proposing points in small area of the space without. proper initialization. Another difficulty,
20

---

Published in Transactions on Machine Learning Research (06/2025)
common to many other LFI algorithms, is the choice of data summary statistics used to define discrepancies. In Section C.3 of the Appendix we consider alternative choices of summary statistics to the ones considered in section 4.2.2 based on auxiliary model approaches to summary statistic construction which are widely used in the ABC literature (Drovandi et al. 2015). The choice of informative summary statistics remains challenging when working with complex models like SSMs. Further exploration of these issues left to future work.
References
Ziwen An, David J Nott, and Christopher Drovandi. Robust Bayesian synthetic likelihood via a semi- parametric approach. Statislics and Computing, 30(3):543 557, 2020.
Alexander Aushev, Henri Pesonen, Markus Heinonen, Jukka Corander, and Samuel Kaski. Likelihood-free inference with deep Gaussian processes. Compulational Statistics and Dala Analysis, 174:107529, 2022.
Giwon Balig, Daniel Evans, Matthew Galdo, and Brandon M Turner. Gaussian process linking functions for mind, brain, and behavior. Proceedings of the National Acodemy of Sciences, 117(47):29398 29406, 2020.
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy. Botorch: A framework for efficient Monte-Carlo Bayesian optimization. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 21524-21538 Curran Associates, Inc., 2020.
Nikolay Bliznyuk, David Ruppert, Christine Shoemaker, Rommel Regis, Stefan Wild, and Pradeep Mugun- than. Bayesian calibration and uncertainty analysis for computationally expensive models using optimiza- tion and radial basis function approximation. Journal of Compulational and Graphical Statislics, 17(2): 270 100-2944220 2008.
Scott D Brown and Andrew Heatheote. The simplest complete model of choice response time: Linear ballistic accumulation. Cognitive psychology. 57(3):153 178, 2008.
Colin F. Camerer. Prospect theory in the wild: Evidence from the field. In Daniel Kahneman and Amos Tversky (eds.), Choices, Values, and Frames, pp. 288-300. Cambridge University Press, 2000.
Adam D. Cobb and Brian Jalaian. Scaling Hamiltonian Monte Carlo inference for Bayesian neural networks with symmetric splitting. In Cassio de Campos and Marloes H. Maathuis (eds.), Proceedings of the Teerth-Sewenth Conference on Uncertainty in Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pp. 675-685. PMLR, 2021.
Patrick R. Conrad, Youssef M. Marzouk, Natesh S. Pillai, and Aaron Smith. Accelerating asymptotically exact MCMC for compulationally intensive models via local approximations. Journal of the American Statistical Association, 15915)) 2016.
Dennis D Cox and Susan John. SDO: A statistical method for global optimization. In Mullidisccilinary Design Optimization: State of the Art, pp. 315 329, SIAM: Philadelphia, 1997.
Peter S Craig, Michael Goldstein, Allan H Schcult, and James A Smith. Pressure matching for hydrocarbon reservoirs: a case study in the use of Bayes linear strategies for large computer experiments. In Case Studies in Bayesian Statistics: Volume III, pp. 37-93. Springer, 1997.
Samuel Daulton, Maximilian Balandal, and Eytan Bakshy. Parallel Bayesian optimization of multiple noisy objectives with expected hypervolume improvement. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (ods.), Advances in Neural Information Processing Systems, volume 34, pp. 2187 2200. Curran Associates, Inc., 2021.
Jiaxuan Ding, Eui-Jin Kim, Vladimir Maksimenko, and Prateck Bansal. Can decoy effects mudge ride- hailhaiinn drivers' preferences for electric vehicles? Available at SSRN 4682413, 2024.
21

---

Published in Transactions on Machine Learning Research (06/2025)
Christopher C. Drovandi, Anthony N. Pettitt, and Anthony Lec. Bayesian Indirect Inference Using Para- metric Auxiliary Model. Statistical Science, 30(1):72 - 95, 2015.
Michael Emmerich. Single-and multiiobjecctive evolutionary design optimization assisted by Gaussion random field metamodels. PhD thesis, Fachbereich Informatik, University of Dortmund, 2005.
M. Evans and H. Moshonov. Checking for prior-data conflict. Bayesian Analysis, 1:893 911, 2006.
Nathan J. Evans. A method, framework, and tutorial for efficiently simulating models of decision-making Behavior Research Methods, 51(5):2390-2404 2019.
Mattoo Fasiolo, Simon N. Wood, Florian Hartig, and Mark V. Bravington. An extended empirical saddlepoint approximation for intractable likelihoods. Electronic Journal of Slalistics, 12:1544 1578, 2018.
Mark Fielding, David J. Nott, and Shie- Yui Liong, Efficient MCMC schemes for computationally expensive posterior distributions. Technometrics, 53(1):16-25 28, 2011.
David T. Frazier and Christopher Drovandi. Robust approximate Bayesian inference with synthetic likeli- hood. Journal of Computational and Graphical Statistics, 30(4):958 976, 2021.
David T. Frazier, Christian P. Robert, and Judith Rousscau. Model misspecification in approximate Bayesian computation: consequences and diagnostics. Journal of the Royal Statistical Society: Series B (Statistical Methodology), (2(2):421-444, 2020.
David T Frazier, David J Nott, Christopher Drovandi, and Robert Kohn. Bayesian inference using synthetic likelihood: asymptotics and adjustments. Journal of the American Statistical Association, 218(5 2832, 2022.
Roman Garnett. Bayesian Oplimization. Cambridge University Press, 2023.
David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke. Automatic posterior transformation for likelihood free inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA, volune 97 of Proceedings of Machine Learning Research, pp. 2404 2414. PMLR, 2019.
D. Gunawan, G.E. Hawkins, M.-N. Tran, R. Kohn, and S.D. Brown. New estimation approaches for the hierarchical Linear Ballistic Accumulator model. Journal of Mathematical Psychology, 96:102368, 2020.
Michael U. Gutmann and Jukka Corander. Bayesian optimization for likelihood-free inference of simulator- based statistical models. Journal of Machine Learning Research, 17(125):1 47, 2016.
Thomas O. Hancock, Stephane Hess, A.A.J. Marley, and Charisma F. Choudhury. An accumulation of preference: Two alternative dynamic models for understanding transport choices. Transportation Research Part B: Melhodoloyical, 149:250 282, jul 2021a.
Thomas O Hancock, Stephane Hess, Anthony AJ Marley, and Charisma F Choudhury. An accumulation of preference: two alternative dynamic models for understanding transport choices. Transportation Research Part B: Methodological, 149:250-282, 2021b.
Marko Järvenpää and Jukka Corander. Approximate Bayesian inference from noisy likelihoods with Gaussian process enulated MCMC, orXiv preprint arXiv:2104.039, 2021.
Marko Järvenpää, Michael U. Gutmann, Aki Vehtari, and Pekka Marttinen. Gaussian process modelling in approximate Bayesian computation to estimate horizontal gene transfer in bacteria. The Annals of Applied Statistics, 12(4):2228 2251, 2018.
Marko Järvenpää, Michael U. Gutmann, Arijus Pleska, Aki Vehtari, and Pekka Marttinen. Efficient Acqui- sition Rules for Model-Basee Approximate Bayesian Computation. Bayesian Analysis, 14(2):595 622, 2019.
22

---

Published in Transactions on Machine Learning Research (06/2025)
Marko Järvenpää, Michael U. Gutmann, Aki Vehtari, and Pekka Marttinen. Parallel Gaussian Process Surrogate Bayesian Inference with Noisy Likelihood Evaluations. Bayesian Analysis, 16(1):147 - 178, 2021.
K. Kandasamy, J. Schncider, and B. Poczos. Bayesian active learning for posterior cstimation. In Proceedings of 24th International Joint Conference on Artificial Intelligence (IJCAI '15), pp. 3605 3611, July 2015.
Marc C. Kennedy and Authony O'Hagan. Bayesian calibration of computer models. Journal of the Royal Statistical Sociely: Series B (Stalistical Methodology), 63(3):425-446) 464, 2001.
John R. Lewis, Steven N. MacEachern, and Yoonkyung Lee. Bayesian Restricted Likelihood Methods: Conditioning on Insufficient Statistics in Bayesian Regression (with Discussion). Bayesian Analysis, 16 (4):1393 2854, 2021.
Wentao Li and Paul Fearnhead. On the asymptotic efficiency of approximate Bayesian computation estima- tors. Biometrika, 105(2):285_ 299, 2018.
Xinwei Li and Prateek Bansal. The importance of response time in preference elicitation: Asymptotic results. Available al SSRN 4782582, 2024.
E. C. Marshall and D. J. Spiegelhalter. Identifying outliers in Bayesian hierarchical models: simmlation- based approach. Bayesian Analysis, 2:409-444, 2007.
Edward Meeds and Max Welling. GPS-ABC: gaussian process surrogate approximate Bayesian computation. Nevin L. Zhang and Jin Tian (ods.), Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI 2014, Quebee City, Queber, Canada, July 23-27, 2014, pp. 593 602. AUAI Press, 2014.
Elina Numminen, Lu Cheng, Mats Gyllenberg, and Jukka Corander. Estimating the transmission dynamics of streptococcus pneumoniae from strain prevalence data. Biometrics, 69(3):748-751) 2013.
George Papamakarios and lain Murray, Fast e-free inference of simulation models with Bayesian conditional density estimation. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (eds.), Advances in Newral Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
George Papamakarios, David Sterratt, and lain Murray. Sequential neural likelihood: fast likelihood free inference with autoregressive flows. In Kamalika Chaudhuri and Masashi Sugiyama (eds.), Proceedings of the Snenty-Seconel International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pp.
PMLR, 2019.
A. M. Presanis, D. Ohlssen, D. J. Spiegelhalter, and D. De Angelis. Conflict diagnosties in directed acyclic graphs, with applications in Bayesian evidence synthesis. Statistical Science, 28:376-397, 2013.
Leah F. Price, Christopher C. Drovandi, Authony C. Lee, and David J. Nott. Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics, (711)1111 2018
Stefan T. Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ull- rich Köthe, and Paul-Christian Bürkner. Bayesflow: Amortized bayesian workflows with neural networks.
Journal of Open Source Software, 8(89):5702, 2023.
Carl Edward Rasmusssn. Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian inte- grals. In Bayesian Statislics 7: Proceedings of the Seventh Valencia International Meeling, pp. Oxford University Press, 2003.
Thomas J. Santner, Brian J. Williams, and William I. Notz. The Design and Analysis of Computer Experi- ments. Springer New York, 2018
Marvin Schmitt, Stefan T Radev, and Saul-Christian Bürkner. Fuse it or lose it: Deep fusion for multimodal simulation-based inference. Xie preprint or 1rX6::2311.10671, 2023.
23

---

Published in Transactions on Machine Learning Research (06/2025)
S. A. Sisson, Y. Fan, and M. A. Beaumont (eds.). Handbook of Approzimate Bayesian Computation. Chapman & Hall/CRC, 2018.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXio:0912.3999 2009.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias W Seeger. Information-theoreti regret bounds for Gaussian process optimization in the bandit setting, IEEE transactions information theory, 58(5):3250 3265, 2012.
Andrew Terry, A.A.J. Marley, Avinash Barnwal, E.J. Wagenmakers, Andrew Hcathcote, and Scott D. Brown. Generalising the drift rate distribution for linear ballistic accumulators. Journal of Mathemalical Psychology, 68-69:49-58, Oct 2015.
Owen Thomas, Raquel Sá-Leão, Hermínia de Lencastre, Samuel Kaski, Jukka Corander, and Henri Pesonen. Misspeciffcattin--rbbss likelihood-free inference in high dimensions. arXiv preprint 0rXiv:2002.09877 2020.
Jennifer S Trueblood, Scott D Brown, and Andrew Heathcote. The multiattribute linear ballistic accumulator model of context effects in multialternative choice. Psychological review, 121(2):179, 2014.
Brandon M Turner and Per B Sederberg. A generalized, likelihood-free method for posterior estimation. Psychonomic bulletin & review, 21:227-250 2014.
Brandon M Turner, Per B Sederberg, Scott D Brown, and Mark Steyvers. method for efficiently sampling from distributions with correlated dinensions. Psychological methods, 18(3):368, 2013.
Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian Schmon. Robust neural posterior estimation and statistical model criticism. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Newral Information Processing Systems, volume 35, pp. 33845- 33859. Curran Associates, Inc., 2022.
Richard Wilkinson. Approximate Bayesian computation (ABC) gives exact results under the assumption of model error. Statistical Applications in Genetics and Molecular Biology, 12(2):129 141, 2013.
Richard Wilkinson. Accelcrating ABC methods using Gaussian procossos. In Samucl Kaski and Jukka Corander (eds.), Proceedings of the Sevenleenth International Conference Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pp. 1015-1023, Reykjavik, Iceland, 22-25 Apr 2014. PMLR.
Simon N Wood. Statistical inference for noisy nonlincar ecological dynamic systems. Nature, 466(7310): 1102 1104, 2010.
Richard Zhang and Daniel Golovin. Random hypervolune scalarizations for provable multi-objective black box optimization. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Con- ference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 11096-11105. PMLR, 2020.
A Derivation of Methods & details of implementation
Table glossary of the main notation used in this paper.
24

---

Published in Transactions on Machine Learning Research (06/2025)
Table 1: Glossary of notation
The model parameter taking a value in parameter space + R.. Data to be observed taking value in y The sampling density of given A The observed value of y A d-dimesional summary statistic to be observed
P(u10)
Wobs
S
p((10)
The sampling density of the summary statistic given The observed value of
Sobs
p(0 (ybbs))
The posterior distribution given the observed y The partial posterior distribution given the observed S
p(0) p (Sobs)
Pr (Sobal 10)
ABC likelihood with tolerance t > 0
(S, Sobs)
A discrepancy between a simulated summary S ~ p(S10) and Sobs
(0))
The expectation of (S Sohs) for S ~ p(S(0
(x (xT, (xT,wT) TT
For two-source data, y consists of component data sources x and w Tobs and Wabs are the observed values of x and w For two-source data, the vector of summaries concatenates summaries T for data source and U for data source w The observed values of the summary statistics T and U Likelihoods for T, U, UIT and T/U for two-source summary statistic data
Yobs = (Ioba Wohs)
S = (T"
Tobs: Uobs o(Tobs)0). (.0)).((Unn plUnraal0) 10)
p(Uobul (Tobs: 0), p(Tobxl PaTabalUabss
AT, Tobs)
A discrepancy between a simulated summary T ~ p(((())) and Tobs
wh (U, Uobs)
A discrepancy between a simulated summary ~ p(U]0) and Uobs
Di(0) D2
The expectation of (T, Tobs) for T TT~ppTT6
The expectation of (U, for U ~~ pPp()))
(S (., (Topx)))) 10), (Uobs
MOBOLFI approximate likelihoods for S, T, U, U and T/U with tolerance t (21,
Pe (Uob Tobs: 0), (Tobel Wobs,
A.1 Gaussian Processes and Bayesian optimization
Next we provide further background on Gaussian processes and Bayesian optimization (BO). To make the discussion self-contained, there is some repetition of definitions and concepts explained in the main text. BO with a Gaussian process surrogate is used in the BOLFI method which inspires the new MOBOLFI approach in our work.
Bayesian optimization attempts to find a global optimum of a function. We will consider Bayesian optimiza- tion for finding global minimum. We want minimize f(0), 0 E θ C R2, where derivatives of f()) are not available, and evaluations of f(-) may be corrupted by noise. BO is most suitable for problems where the dimension of 0 is not too large, although high-dimensional BO is an active area of current research. We model the noisy evaluations of f(0) with a "surrogate model" describing uncertainty about f(+) given the function ovaluations previously made. A common surrogate model in BO is a Gaussian process, and this is used the BOLFI method.
We establish some notation and definitions first. Suppose that θ is an n p matrix, with ith row 0, O For any function g θ E, write g(@) for the vector obtained by applying 9()) to the rows of O, Le. g()) (g("))....: TT..)) Suppose that θ an m x p matrix, where the ith row is 0, E O. For any function h O x θ + RR we write h(θ, O) for the n x m matrix with (i,j)th entry I(00, 0,)). 0,....,,, E θ, denote by 01m the n x p matrix with ith row 0,, random function f(-) defined on is a Gaussian process with mean function j θ + IR and positive definite covariance function C x θ - R if, for any n, and any f.....,o, O the random vector f(1m.) multivariate normally distributed, (01.a) N(r(01sa)) C(O1mn,
For
25

---

Published in Transactions on Machine Learning Research (06/2025)
Suppose we observe the Gaussian process f()) with noise at points . On O. are
The noisy observations
2; = ((0)) + 6i2 i == ,,...,n, n,
(29)
where E Fd (0, o²), for some variance o² > O. In many applications, the Gaussian noise assumption can be reasonable in the vicinity of the minimizer with an appropriate transformation. Write ²gn = (2),..., =n) We are interested in describing uncertainty about f((*) for some 6" E O, given the noisy observations zzn.
Since /(.) is a Gaussian process with mean function pa(-) and covariance function C(:,)), (0(0m)) [(0)) is
multivariate normal, and because of the independent normally distributed noise in equation 29 (z f(0)
is easily shown to be multivariate normal. It follows that [ZEn ~ N( 0 where
(0)) 1a(0)) C( 01,)) {C(01:n; 1mm) 022 (zgn (u(0m=))) (00 C(f", 00) C(8", .1m) { {C(01mm P1:)) r²13 C(01:n 8*).
(30)
(31)
The uncertainty quantification provided by equation 30 and equation 31 can be used to decide which 0* should be used to obtain a further noisy observation
x" = f(0) +e, e' N(,,0))
having the most benefit in the search for the minimizer. 0 is chosen through an optimization, of the "acquisition function". This acquisition function can be the expected loss for some formal decision problem where the expectation is taken with respect to the the Gaussian process uncertainty about f(.), or it might be chosen using more heuristic reasoning.
In the BOLFI method, Gutmann & Corander 2016 suggested using the lower confidence bound acquisition function [Cox John 1997).
An = (0) - o
(32)
where
n2 2 πt²
2 = 2 log
3en
with the default value e, = 0.1. With such design of the na the cumulative regret of a GP using Upper Confidence bound as the acquisition proved to be bounded with high probability (Srinivas et al., 2009). In the day care center example discussed in Section 4, we implement BOLFI using n However, in the two examples presented in the main text, the large number of parameters p results in a large na This increases the probability that the acquisition point is not close to the global optimum, necessitating more iterations and longer training times for the implementation of BOLFI. Additionally, this approach risks proposing points outside the range of the uniform prior, potentially leading to numerical problems. To ensure numerical stability and efficient training, for the two examples in main text, we replace v/ with
nn
Th = 2 log
36
We posit that using x specifically for these examples strikes an effective balance between exploration and exploitation.
Further discussion of the intuitive basis for this acquisition function is given in the main text. We have described the process of choosing algorithms. To completely specify BO algorithm, we need to describe a method of initialization (where some initial function values are obtained, perhaps using a space-filling design on θ), a stopping rule (most simply we might stop when a given budget of function evaluations is exhausted) and a method for estimating Gaussian process hyperparamcters, and updating estimates as the algorithm proceeds. For further details on these practical issues see (Garnett) 2023). For the two examples in the main text, we considered different initializations, with details given in Sections 2.1 and 3.2.
location for taking
noisy function evaluation in BO
26

---

Published in Transactions on Machine Learning Research (06/2025)
A.2 Related work on BOLFI
The BOLFI approach of Gutmann & Corander (016) has been very successful for simulation efficient estimation of posterior distributions with expensive simulators, and the basic method has been further de- veloped in a number of ways. Järvenpää et al. (2018) explores the importance of the Gaussian process formulation used in the BOLFI framework, including transformations of the discrepancy, heteroskedastic classifier Gaussian process formulations for likelihood approximations, and different utilities for Gaussian process model choice. Järvenpää et al. (2019) go beyond generic acquisition functions from the Bayesian op- timization literature and develop alternatives tailored to simulation-based inference problems targeting the reduction of posterior uncertainty. Järvenpää et al. (2021) consider Bayesian optimization for likelihood-free inference with noisy log- likeliid evaluations and batch sequential strategics amenable to parallel compu- tation. Aushev et al. (2022) consider replacing the Gaussian process used in BOLFI with a deep Gaussian process, and demonstrate that this can result in more accurate posterior approximations, particularly when the target posterior density is multi-modal. [Kadaasmy et al. (2015) uses Bayesian optinization with some novel acquisition functions to query an expensive to evaluate likelihood, before estimating the posterior using the cheap to evaluate surrogate. A generalized Bayesian version of BOLFI, which is suitable for misspecified models and when the summary statistic dimension is high, is described by Thomas et al. (2020).
The use of a Gaussian process or other "surrogate" models to obtain a likelihood approximation is not only used in the context of BO methods. Wilkinson (2014) considers a Gaussian process surrogate for a synthotic likelihood (Wood. 2010: Price et al. 2018) in so-called history matching algorithms (Craig et al., 1997). Meeds & Welling (014)) considers Gaussian process surrogate models for summary statistic means and variances for a synthetic likelihood, and adaptively acquire simulations in order to reduce uncertainty in sxcopt/reject decision for posterior simulation. A similar more sophisticated approach
a Metropolli-Haatting has recently been considered by Järvenpāä & Corander (2021), where the authors use a Gaussian process surrogate for the log-likelihood itself, and are more explicit about acquisition rules for augmenting the training set for the Gaussian process. A surrogate model can be useful too in applications where exact likelihood
caleulations can be made, but are expensive. See, for example, (Kennedy & O O'Hagan 2001 Rasmussen
2003 Bliznyuk et al. 2008 Fielding et al. 2011 Conrad et al. 2016) among many others.
A.3
Implementation of MOBOLFI
Here give details of multi-objective Bayesian optimization which are not given in the main text. Multi-objective BO is used in implementing the MOBOLFI method developed in the main text. Once again, there is some repetition with definitions and concepts in the main text to make the discussion self-contained
Let f(0) = components of f'(-). In general there is no value 0" E where all components are minimized simultaneously. Multi-objective optimization methods approximate the set of "non-dominated" solutions which are not ob- viously inferior to other solutions. A value 0 + dominates 0' E θ if f₂(0) S (0),), j ,,....K, K, with the inequality being strict for at least one j. The dominated solution is inferior in the sense that there is another point at which f'()) is strictly smaller along some dimensions and no larger for the other dimensions. Multi-objective optimization algorithms try to find the Pareto optimal set of non-dominated points in O.
(0)) be a multivariate function, for which we are interested in minimizing the
Numerical multi-objective optimization methods obtain finite approximations to the Pareto set. The Pareto frontier is the set of opai function values obtained by the points in the Pareto set. Multi-objective Bayesian optimization Garnett (2023) Section 11.7) uses surrogate models to implement multi-objective optimization for expensive to evaluate functions, possibly observed with noise. Similar to Bayesian optimization with a scalar objective, the representation of uncertainty given by the surrogate is used to efficiently decide where to perform the next function evaluation.
Multivariate Gaussian processes are a common choice of surrogate for multi-objective Bayesian optimiza- tion, and we give some background and notation now, extending the discussion of Section 2.3. Suppose that g : > K and that θ an n X p matrix with ith row 0, E e. We write g = (g(0,) ...., g(,)))(())))) e R Kn, Let h: O x + k be K x K matrix-valued function. Let 7 be an n x m matrix with ith row 0. We write h(0,0 0) for the partitioned matrix with n block rows and block columns where the (i,
27

---

Published in Transactions on Machine Learning Research (06/2025)
block entry is h(O,; B,) E R A random function f() (((() is a multivariate Gaussian process with mean function J - RRK and positive definite covariance function C: x + R if for any n and 01...., On, f01mm) is multivariate normal, N (m(O1:)) C(0mm, 0mm)).
Once again extending the discussion of Section 2.2, suppose we observe values of f()) with noise at 01.... 8n E e to obtain
2, = f(0) Es.
(33)
where now 2, E K and 6, add N(O, Σ) where E RR is some positive definite covariance matrix. As our discussion of the univariate case, for some 0* O, and writing ²gn = (2] ",,", (z ((0))T)" is multivariate Gaussian and the conditional density of ((0") is multivariate normal with mean vector and covariance matrix
Pn (0") = 1=(0)) + C(0", ,1I..) (C(O1:n 01:)) + Dr, (zga u((1m) E, (+*) C( 0)) 0(00000mm (C(61mm f1m) + Dn (E)) C(01m, 0*),
(34)
(35)
where (E) E RaKxnK is the block diagonal matrix with the K x K diagonal block entries equal to E. For reducing computation cost, our later numerical experiments we consider multivariate Gaussian processes where the components of f(-) are independent, but we consider correlated noise in equation 33 i.e. Σ is not diagonal.
The noise 6 can be observed by the variation of repeated simulation {=ري)ي= give input 0, Therefore, the dependent noise covariance matrix Σ is estimated by the covariance of a repeated finite simulation sample (An. where Air.j = D(6,.) + Cis; for some 6ig; is a simulated bivariate noisy discrepancy. For results Un, (0.)).
in this paper, we set ng= 100 and 0ig arg min (A, - (di-Mng(0 (0.) (2, (T.,,))TT.,,
Given the uncertainty quantification provided by cquation 34 and equation 35| if there is finite set of points, say 01...., O2, approximating the Parelo set, with corresponding approximation Jn of the Pareto frontier, expected hypervolume improvement (EHVT) measures the volume of the space dominated by the current approximation of the Pareto frontier and bounded below by a reference point, the so-called Pareto hypervolume. EHVI was firstly used as an acquisition function in multi-objective Bayesian optimization by |Emeerich| 2005), where an expectation of the hypervolume improvement is taken with respect to the surrogate model uncertainty to define the acquisition function. In the MOBOLFI method, multi-objective BO is applied to a vector of expected discrepancies, and the simulated discrepancies are noisy. For this reason, we use the noisy expected hypervolume improvement (NEHVI) (Daulton et al. (021)) for the acquisition function. NEHVI implements a Baycsian treatment when calculating EHVI, integrating uncertainty about the Pareto frontier. This makes NEHVI more suitable than EHVI as an acquisition function with noisy data.
Gutmann & Corander (016)) chose the tolerance t in the univariate BOLFI method as the q-quantile of Ai,...,Any where q 10, In the bivariate MOBOLFI method, extend the choice of tolerance (t1, to the 2-dimensional vector q-quantile of A ..., Any. where q 1)3 Given that the evaluation of ABC approximate likelihood involves f, we also did comparison study over different q-quantile tolerance for 3 different examples, by letting g performance of inference, varying by examples.
Our visuals present the affect of tolerance on
In multi-objective Bayesian Optimization, one practical difficulty is the scaling of objectives. For a fixed diagonal matrix V, instead of applying multi-objective BO a vector of noisy discrepancies A, we could apply it to V-1 A instead, and in general the results are not invariant to the choice of V. Figure 2 studies the effect of scaling on inference in toy example with V-1 diag(w, 1), where w is scalar weight. Since only two objectives involved, we follow the notation in Section 4.1, where A = (Ai(X (X((,, Xx), (X)))x))) denotes the joint objective over two data sources X and W. From the figure, we find that with different w, if we simply add the discrepancy to get a univariate discrepancy for the BOLFI posterior the approximate posterior differs markedly, while for MOBOLFI the approximate posterior is less sensitive to w. Not doing scaling equivalent to setting the V-1 Iz, i.c. the red curve in this figure, which not the choice with the best performance.
Az,(() WW))
28

---

Published in Transactions on Machine Learning Research (06/2025)
One classic scaling method in the machine learning literature is normalization by using Mean absolute deviation (MAD). Algorithm
| presents the detailed steps of doing scaling of joint noisy objective A for implementing BOLFI. They key idea of Algorithm g is to put each elements of A on a similar scale. The use of Algorithm I does not always result in the best performance. In the plot 2a the green curve with w = 0.7 is chosen by Algorithm | (rounding to decimal places). The scaling with best performance is example specific, depending on the information of each data source, the LFT method and the BO acquisition function. We leave further investigation of optimal scaling approaches to future work. For the results of this paper, by default choose the scaling the outcome (rounded to 1 decimal places) from algorithm | using n = 100. The auxiliary model example (orange curve in Figure 13) defines a discrepancy the score vector of an auxiliary model, and the scaling in this example is slightly different and will be discussed in Section C.3 of Appendix.
Algorithm 1 Scaling of discrepancies in (MO)BOLFT
Require: Prior πT)) of 0, target function = (A1, AK)" to scale, nunber of sample size n Sample 063 ~ (0), = ,,....,
For each 0,, evaluate the corresponding target function A(0(0)) (B(mn))
..... (A((0)))...,x(((())))))
For 1, -..... K, write e(ار = MAD(A)(0(), ( Write V = diag(u("), ,.. ,((K)
Scale the target function A = V-IA
(In BOLFI, A = A ... Ak, and the scaled discrepancy is Ascale Vi A1 + ... vk)
B Experiment setup and extra findings - Toy example
We give some further details of the implementation of MOBOLFI and some additional experiments that were not included in the main text due to space limitations. Code to implement all experiments can be obtained
git
B.1 Experiment setup toy example
In the toy example in the main text initial training data of 100 observations was used, and BOLFI/MOBOLFI was used to train I-dimensional/2-dime dimeesiioa surrogate GP models respectively with 200 BO acquisitions. Parameter samples are drawn from the approximate posterior distribution (Section 2.3) using Hamiltonian Monte Carlo (HMC) with tolerance set to be the 1% quantile of the training data dis- crepancies. This is done element- wise for each discrepancy to get the vector of tolerances for MOBOLFI. We implement HMC using the hamiltorch package (Cobb & Jalaian| 2021), running four chains for 8,000 itera- tions each, using step size 0.1 and one step within each proposal (which corresponds to a Metropolis Metropolis-Hastings adjusted Langevin algorithm).
For MOBOLFI, the objective function optimized is
D(0) =
where X is the synthetic data and X is the observed data for the first data source, A X²) is the discrep- ancy for the first data source, W is synthetic data and W is the observed data for the second data source, and A2 Wo) is the discrepancy for the second data source. The definition of the discrepancies is given in the main text. For implementing BOLFI, the training data is {0,, w- (x ( (W)()) W.) Wo) where 0.4 is the scaling with best performance from mumber of alternatives (scc figure 2a in Section A.3 of this Ap- pendix). We use the lower confidence bound acquisition function for BOLFI with = 0.1 in equation 32 For implementing MOBOLEL, the training data is (0,, (w (x X)),, (W Wo) We use the NEHVI acquisition function (Daulton al. 2021) with a reference point min,((A, (X(,, (x), (W((9, W))) 0.1. Each optimization in the multi-objective BO iterations approximates the optimum of the acquisition func- tion from 100 candidate samples in 10 restarts. To avoid observations with discrepancy values which are too
29

---

Published in Transactions on Machine Learning Research (06/2025)
1.0
post (50 iter)
(100 ster)
1.2
quantile
quantile
(150 iter)
01 quantile prior
0.8
(200
1.0
prror
Deruty
0.6
0.8
0.6
0.4
0.4
0.2
0.2
0.0
0.0
-2
0
1
2
01
(a) Different iterations for GP training
(b) Different choices tolerance
Figure 9: Approximate posterior for the toy example given different numbers of training iterations and tol- erances. The left column shows the MOBOLFI approximate posteriors given 50/ 1100/150/50 BO iterations. The right column presents the approximate posteriors for threshold set to 20% /1% quantiles of the training data discrepancies. All approximate posterior densities are kernel density estimates from HMC samples. The red dashed line is at ffrue = -0.7.
large (which affects GP training), we set a mild rejection criterion to filter out observations higher than a threshold, i.c. the 99%-quantile (rounded to 1 decimal places) obtained from simulated discrepancies that are drawn without any filtering criterion.
The botorch package maximizes the supplicd objective function by default, and so implementation of BOLFI we maxinize the negative expected discrepancy. Guttmann & Corander (2016) minimizes the lower confidence bound as (32) in their work, and when maximizing the negative expected discrepancy it is equivalent to use the upper confidence bound as acquisition function,
(0) = (0) + (0))
(36)
B.2 Extra findings toy example
In addition to the experiments in the main texl, Figure 9 shows MOBOLFI approximate posteriors of 0 for different choices of the number of iterations the BO algorithm and for different choices of tolerance. The approximate posteriors are kernel density estimates from HMC samples. The left-hand column of the figure demonstrates that the MOBOLFI approximate posterior variance decreases with more BO iterations. This is expected, since the surrogate model uncertainty contributes to the uncertainty in the MOBOLFT approximate posterior. The right-hand column of the figure demonstrates that the MOBOLFI approximate posterior is closer to the true posterior when smaller quantile of the training data is used for the tolerance. In Figure (b), 200 BO acquisitions were used. Gutmann & Corander (2016) suggested using the 5% quantile of the training data, but a 1% quantile attains comparable performance in this example. Both the selection of tolerance levels the likelihood approximation and the number of BO iterations used are crucial for achieving accurate approximate posterior distributions using MOBOLFI.
MOBOLFI differs from BOLFI by using multiple discrepancies in the BO algorithm. To explore the sensitivily of inference to including multiple data sources when the model for some data sources does not depend on some parameters, we modify the simulator as follows. The first 8 elements of Bure are shared by both X and Y. However, the 9th element of Ourue contributes to X only, and the 10th element of Ouroe contributes to Y only. Specifically,
0 1(
X, N(x 0,, 11, N
w = Ocl + od)(l)
(37)
0x (0true, .... Ocrue,8 8, (true_9)²
0w = (0true, Orruc, 0true,10) T
30

---

Published in Transactions on Machine Learning Research (06/2025)
We apply MOBOLFI with the same experiment setup in Section B.1 of this Appendix. The MOBOLFI and BOLFI approximate posteriors are presented in Figure TO We focus on the parameters 0g. 010, which are parameters influencing the distribution of only one of the data soruces. In plot 10a the MOBOLFI approximate posterior, leveraging correlated noise between data sources, still outperforms BOLFI inference of 0g. 010- In plot [ob. for 010 not depending on X, the MOBOLFI approximate posterior conditional on of X does not obtain smaller variance than the BOLFI approximate posterior. That not surprising see, since information from X viewed as useless and redundant for inference of 0 On the other hand, MOBOLFI still performs better than BOLFI in the inference of 0g. Plot [0c shows that approximate MOBOLFI posterior conditional on W only does not obtain smaller variance than the BOLFT approximate posterior, for 0g not depending on W.
1.0
1.0
BOLFI
BOLFI
0.8
post
0.8
post
true
Fur
Dentt
0.6
Demty
0.6
0.2
0.2
0.0
0.0
09
100
(a) MOBOLFI BOLFI
1.0
1.0
MOBOLFI
BOLFI
BOLFI
0.8
0.8
Fuue
Densty
Densty
0.2
0.2
0.0
0.0
09
010
(b) Information in X
1.0
MOBOLFI BOLFI
1.0
MOBOLFI
BOLFI
0.8
0.8
true
t've
0.6
Devott
0.6
0.2
0.2
0.0
0.0
010
(c) Information in W
Figure 10: Approximate posterior for the updated toy example with two parameters 0g 010 solely depend- ing on two data sources X, W, respectively. The row 10a shows the MOBOLFI and BOLFI approximate posteriors under the updated simulator, colored by blue/orad alue/eranne. The green curve is the real posterior, the red curve is the prior and the dash red line is the true value 6g, 010. In row 10b 10c the MOBOLFI and BOLFI approximate posteriors calculated by marginal likelihood of X/W are given by blue/orange lines. The green curve is the real posterior of X/W respectively.
31

---

Published in Transactions on Machine Learning Research (06/2025)
C
Experimental setup and extra findings: MLBA example
C.1 MLBA closed form likelihood function
Consider docision-making for a single individual first. For a drift rate % ~ ((a,, s²), a = 1..., M, the probability density function (pdf) of the time taken for the accumulator to reach the threshold X (see Brown & Heathcote (2008), Appendix A for a derivation):
ft(t) 1 do (* A A-1da)
ső (x (x-A-1de -dd + do (x tdn,
s x ((lda) Idg =
)}
(38)
ts
ts
ts
18
where <(.) and (() are the pdf and edf of the standard normal distribution. Its associated cumulative density function (odf) is:
Fa(t) = + X
Ido
X-A
A - Ida )
X ldo ( x Idn
A
ts
A
ts
(39)
X
tda
X - tdl
+ Φ
A
ts
A
ts
The joint pdf of CH = a and RT t + To
MLBAjoing joint (CH = a, RT = )- = fa(t)llogaa (1 Pb(t)),
(40)
and the marginal pdf of CH is
MLBAA (CH (CH= 1/ Sa(1) IT (1- Fo(t))dt,
(41)
beC
where C is the choice set.
We follow the adjustment suggested by Terry et al. 2015). where the drift rate follows a truncated normal distribution: Va ~ (d 8², 0,, oc), with and oo as lower and upper bounds of the support. This distribution helps correct the original pdf and edf in equation (38 and equation 39 Moreover, the truncated normal distribution resulted in superior performance in empirical experiments conducted by HHancock et al.) (2021a). Therefore, both ft) and (t) are additionally divided by factor ((())) in equation 38 and equation
39
C.2 Experimental setup MLBA
When working on synthetic data, we firstly construct training dataset (0(*) (RTOO,CH)) of 100 observations.
In BOLFI/ MMOLLL,,,,, univariate/bivariate GP is trained with 1000 iterations. Likelihood approxi- mations for samples were generated from the approximate posterior distributions using De-MCMC. These samples are compared to samples from the posterior with the closed form likelihood 40 above.
use tolerance given by the 1% quantile of the training discrepancies, and
For MOBOLFI, the objective function optimized is
D(0) = (RT,, RT))), FLACCHH,
where RT, CH are synthetic response time and choice data respectively, RT and CH are the correspond- ing observed data, (I(RT, RTT) is the discrepancy for the response time data and A2(CH, CHo) is the discrepancy for the choice data. The definition the discrepancies given in the main text. In BOLFI, the training objective w- (RT, RT*) A2(CH, CH²), where w 0.7. In MOBOLFI, the training objective is the 2-dimensional vector (w- (RT,, A2 CH*)). In BOLFI, the acquisiiin function is defined to be the Lower Confidence bound with variance weight 0.1 defined in equation
32

---

Published in Transactions on Machine Learning Research (06/2025)
In MOBOLFI the Bayesian Optimization acquisition function is chosen to be the qNEHVI (see attps://botordds wttp:/p.acag....siiaiaeeses http:///a.aaaiiii.oo.....ooooooo www.otorcch.qqussiiiiooo.. multi oww.ooiie...oonetee.ooo. carlo. htgppp/cceddtyeeooooommmmmrrrree
for details). Both MOBOLFI and BOLFI surrogate models are
trained with 1000 BO acquisitions.
For faster convergence, we set. different hyperparameterss for running De-MCMC to sample from MOBOLFI/BOLLE and the posterior with closed form likelihood. When sampling from the posterior for closed form likelihood, we set 9 chains, sample size 20000, burn-in size 18000 and migration rate 0.5. When sampling from BOLFI/MOBOLE 13000. The coefficient y in the Dc-MCMC algorithm is set to be the fixed constant 2.38/ /2 nog where no = 6 is the number of parameters we infer.
approximate posteriors, we set 9 chains, sample size 16000 and burn-in
Similar to the toy example, we set some rejection criteria in sampling the initial training set for BO. Points with values of any one of the discrepancics greater than a threshold, i.c. the 99%-quantilc obtained from simulated discrepancies that are drawn without. any filtering criterion, are discarded. In MLBA simulation we observe that there are some prior samples g that simulate HT = {KT .... 320) where RT, RT Vj E {1..... 320), where RT denotes the sample mean of RT" Ideally the prior should be constrained to avoid such degenerate regions of the parameter space, but these regions are not casily characterized analytically. Hence we also filter out 060) such that Var( < Var(ItT"). 0.7.
C.3 Extra findings MLBA example (Synthetic data)
Similar to the toy example, we build MOBOLFI approximate posteriors given different number of iterations for the BO algorithm. Figure II compares kernel density estimates obtained from samples from the approximate posteriors. Due to the complexity of the MLBA simulator, need many more iterations in the BO algorithm than for the toy example to obtain good approximations. The Figure shows that 100 iterations in the BO algorithm is too little, but the approximation to the marginal posteriors has mostly stabilized after 300 iterations for most parameters.
We also investigate the affect of tolerance on MOBOLFI performance in MLBA. In Figure 12 we compare the MOBOLFI approximate posteriors for tolerances specified as 10%/5%/11 /1% quantiles of the training discrepancies. Given enough Bо iterations, the 1% quantile is the best choice. Although the performance of MOBOLFI is sensitive to the choice of t, setting the tolerance to 1% or 5% quantile of the training data works well across many problems.
When a summary statistic is used in the definition of the discrepancy, its choice is another factor that affects the performance of MOBOLFI. For summary statistic based LFI methods, the choice of summary is important, because good summaries can reduce the computational burden with little loss of information. Tn Figure 13 we evaluate MOBOLFI approximate posteriors obtained using two approaches. The first approach uses the discrepancies discussed in the main text. The second approach changes the discrepancy used for the choice dala, by defining summary statistics using an auxiliary model (e.g.. Drovandi et al., 2015). Specifically, we adopt the score vector evaluated at the maximum likelihood estimate (MLE) for the observed data for a multinomial logit (MNL) model as the choice data discropancy.
To explain further, the MNL model is a discrete choice model based on random utility maximization theory (McFadden, 1974). In the MNL, the choice probability for the alternative is
exp (Va) (V.)
Po
(42)
where Va is systematic utility of the alternative a, represented by set of attributes X (X1, Xx) and corresponding parameters { = (C......k), where K is the total number of parameters. Since it has a closed-form likelihood function, the parameters can be estimated by maximum likelihood estimation. Write (C for the likelihood for the MNL auxiliary model for choice data CH and parameter g. The score function is (A(H::)) = Ve Iog E). Writing CH for the observed choice data, and & for the MLE for the observed choice data, we have SA(CH::)) R O, and for simulated choice data CH we use as summary
33

---

Published in Transactions on Machine Learning Research (06/2025)
300 heers
120 Rer
2.5
Ferr
300 Ker
320 bara
2.0
baney
1.5
1.0
0.0
-0.5
0.0
0.5
1.0
A1
log(x A1
0.5
1.50
Frre
300 BAD
0.4
1.25
1.00
0.3
0.75
0.1
0.25
0.0
-30.0
-27.5
-25.0
-22.5
-20.0
-175
0.00
81
82
0.6
0.5
0.5
2000
0.4
1200 Itera
0.4
Cenetr
0.3
0.1
0.0
10.0
0.0
7.5
02
5g
Figure 11: Approximate posteriors for MLBA example given different numbers of iterations for the BO algorithm. Each plot shows approximate marginal posteriors for one parameter of interest. The densities shown are kernel density estimates obtained from MCMC samples. The dashed red line shows the location of the corresponding true parameter.
statistios (CH; (i.e. we use the (K 1) score vector for the data CH evaluated at the MLE & for the observed choice data as the vector summary statistics).
Each component of the score vector has different scale from the others because the paramcters & correspond to different. scales of attributes (e.g., driving range and purchase price). Therefore, the score vector S (A should be converted into a scalar value be used choice data discrepancy A2(CH, CH*). We scale parameter-speeitt components of the score vector put all components on similar scale as follows. 1,..., For each i, evaluate the corresponding 5( (), where CH() is the choice data simulation given g(1) Then, we apply Algorithm [ to obtain the scaling weight V of the parameter-specifie K-components of the score vectors with sample size n 100 and target function Sx(=E). Finally, the choice data discrepancy is defined as (2(CH, CH²) V {)
Sample from the prior g(1) T0)),
vector
To implement the MOBOLFI approach, we nced to calculate the training objective by calculating two data source-specific discrepancies. The response data discrepancy is the one used in the main text, and the joint discrepancy is (Ai (RT, RT"), V A2(CH,CH*) where V₂-1 is obtained by Algorithm
34

---

Published in Transactions on Machine Learning Research (06/2025)
0.83 quamte
0.03 quandle
2.5
0.15 avantde
0.05 quaritte
quantll
0.3 quandle
2.0
Derety
1.5
1.0
1
0.0
0.0
1.0
3.75
4.00
4.25
4.50
4.75
5.00
A1
logix
A)
0.5
0.00 gimmle quamii
0.01 quamie auantic
0.4
1.25
1.00
0.75
0.1
0.25
0.0
0.00
P1
B2
0.01 guande
quaniie quaritle
0.5
0.5
0.4
Plue
0.4
Denatty
0.3
0.3
0.1
0.0
0.0
62
6g
Figure 12: Approximate posteriors of MLBA example given different levels of quantile of training datasct the tolerance 6. The blue/orange//ree grren curve represents the MOBOLFI approximate posterior (with 1000 iterations training) using 40%/5%/11 5%% quantile of the training dataset as the tolerance t respectively. The dashed red line shows the location of the true parameter.
In our experiment, generated n 100 MLBA datascts for K parameters of the MNL auxiliary model and obtained Vi- (128.4, 226.2, 1656.6, 519.9, 2269.1) and V/- 0.0007, respectively.
In Figure 13 we compare MOBOLFI using the choice discrepancy from the main text with the auxiliary model choice discrepancy (MOBOLFI AUX). The closed-form likelihood (MLBA) is used as benchmark model. Given enough training, the MOBOLFI and MOBOLFI AUX show different advantages in approx- imating the posterior with regard to point estimate and posterior variance. The lower posterior variance indicates more robust estimates, while the point estimate closer to the true parameter indicates accuracy for decision-making in practice (c.g., estimate of electric vehicle adoption rate in this MLBA case). Since the MLBA parameters behaviorally related to both reaction time and choice decision, the MOBOLFI AUX affects all paramcters. For \1, $1, and log(x AA, the MOBOLFI AUX outperforms the MOBOLFI by providing a slightly better point estimate while reducing posterior variance in parameler inference, in- dicating the validity of adopting summary statistics based on an auxiliary model. For B2, 82, and ő3, the approximate postcrior distributions from the MOBOLFI are superior to those from MOBOLFI AUX with similar posterior variance but much better point estimates. This result suggests that the performance of MOBOLFI and MOBOLFT AUX comparable.
35

---

Published in Transactions on Machine Learning Research (06/2025)
MOROLFI
MOSSOLM MILMA
MOGOLII
Daraity
2
1
-0.5
0.0
0.5
1.0
4.0
4.2
4.6
5.0
A
lng(x
A1
MOROOL AIX MROHO.
MORO F1 AUX
1.0
0.6
2
0.4
1
0.0
0
81
B2
MIGBOLFI AUX
MICIOLFI AUX
0.5
Bruit
0.5
0.4
0.3
0.2
0.2
0.0
0.0
D2
Figure 13: Approximate posteriors of MLBA example estimated by different data summary for response choice. The blue and orange curves represent the MOBOLFI approximate posterior (with 1000 iterations training) using the original distance (MOBOLFI) and score function an auxiliary model (MOBOLFI AUX) as the data summary, respectively. The red dash line the value of true parameter gtrue
C.3.1 Misspecified MLBA
MOBOLFT approximate posteriors of the parameters for the misspecified MLBA scenario of Section 4.2.3 are shown in Figure 14 For ÅT, both data sources contribute to the inference, and the posterior calculated using single data sources are similar to the posterior calculated by the joint likelihood. For log(x - AA,, the definition of determines that inference of it mainly depends on the response time data RT. It is unsurprising to see that the posterior calculated by the joint likelihood has similar location the posterior caleulated conditional on RT only. The posterior calculated conditional on CHH only provides much less information, having large variance. For 3,, the true values for the different data sources lie in the tail of the posterior calculated using the joint likelihood, while the posterior calculated conditioning data source are similar and consistent with the true values from both data sources. As discussed above, there is a set of local optimum points that could obtain similar performance in evaluating the discrepancy to gtrue and stmp. The marginal approximate posterior of 31 calculated by the joint likelihood puts weight on points other than fe and fhg consistent with the data. For B3, we found the posterior calculatod using the joint likelihood is closer the posterior calculated conditioning only on the CH data but with variance larger. This differs from the behaviour for 31. We believe this because given the fixed Bg -6,
only
36

---

Published in Transactions on Machine Learning Research (06/2025)
1.50
MODOLII
MOBOLF
MODOLII MOROLFIR IT
1.25
MIDBOUT
MOBOLII
2.0
1.00
1.5
0.75
0.50
0.25
0.00
0.0
-0.25
0.00
0.25
0.50
0.75
1.00
1.25
logix A)
MOROLA
MOROLE
1.0
MOBOLPIE IE
0.4
0.8
0.3
0.2
0.1
0.2
0.0
0.0
81
82
0.30
MORUH
MOROLI MOROLFI
0.25
0.25
0.20
TRUe
0.20
mane ch
Dontte
0.15
0.15
0.10
0.05
0.00
0.00
62
Figure 14: Approximate marginal posteriors for MLBA example under misspecification. The green curve represents the MOBOLFI approximate posterior calculated using the joint likelihood of two data sources. The orange/blut blue curve are the MOBOLFI approximate posteriors calculated by the approximate likelihood of RT/CH data only. The orange/bli /blue dash line shows the location of the true parameters f and gtrue
the values B2 in f and fa are close to the value of 33 but away from 31. This suggests that in the synthetic data example, attribute 1 is more important in evaluating evidence alternatives which directly affects the response time. It is not surprising then see that 8g measures influence of the less important attribute 2 and appears to be more sensitive to the choice data CH than B1. For 82 and 8g, the alternative specific constant is added into the drift rate mean calculation directly, unrelated to the attribute pairwise- comparison, which is the main source of variation of drift rate mean. Therefore, the posterior calculated using the joint likelihood appears similar to the posterior calculated conditional on CH only. Inference of 82 and 83 should mainly depend on the CH data. However, we do obscrve that the conflicting information from RT does brings more uncertainty to the approximate posterior calculated by joint likelihood, which shows larger variance than the approximate posterior calculated conditional on CHH only.
C.4 Extra findings MLBA example (empirical case)
Figure 15 shows the difference between the posterior mean and MAP estimates from the MOBOLFI ap- proximate posterior and closed-form true posterior. The largest diffferences for dICEV and log Brc.
37

---

Published in Transactions on Machine Learning Research (06/2025)
Point estimates from the MOBOLFI approximate posterior have smaller value of 8ICEV and larger value log, Brg, when compared to estimates from closed-form posterior.
Cloeed form
MOSOLI with 6
0.8
MOBOLI mean with
1.5
MENROL MAP with &
cloaed form
0.6
Density
1.0
0.4
0.5
0.0
0.0
1.0
1.5
2.0
2.5
3.0
3.5
SICEV
log(BTc)
Figure 15: The posterior density and point estimates of MOBOLFI and closed-form based CH and RT data sources.
D Day care center example
(2013) measuring infections of bacteria strains over individuals in number of day care centres. We refer the reader to Numminen al. (2013) for detailed description of the model. The data are indicators for whether individuals are infected with different strains at different times. We use the simplified model in Example 3, Gutmann & Corander (2016), and in their notation R, (t) is the rate of infection with strain 8 at time t. Rt(t) is informed as weighted sum of a probability for infection from within the day care centre (which may vary by strain and time) and a probability of infection from outside (which varies by strain). The weights the weighted sum are unknown parameters, denoted I and A. There of infection with strain between situations when an individual already infected with another strain, compared to when they not. The set of all unknowns is 0 (3,
Next we consider data from Numminen
further unknown parameter, which we denote by d, which controls the relative rate
consider synthetic data where the the number of day care centres, the number of attendee individuals sampled in each day care center and the number of strains are set to 29, 36, 33 respectively. We generate the observed data from the model with girue = (3.6, 0.6, 0.1), and implement MOBOLFI with 50/150/250 BO iterations and tolerances set as 20%/5%/1% quantiles of (S,, Sq. S4)
Following Gutmann & Corander 2016)
the training data discrepancies. We consider the same four summary statisties Numminen et al. (2013), which can be computed for each daycare centre and averaged. The prior is the Uniform distribution illustrated in Example Gutmann & Corander (2016). A discrepancy for BOLFI is then obtained by the Euclidean distance between simulated and observed summary statistics. For the MOBOLFI method, we consider two discrepancies: the first is the Euclidean distance between the first and fourth summary statistics for simulated and observed data, and the second is the Euclidean distance between the second and third summarics. The reason for this partitioning is that the first and fourth summaries represent diversity prevalence strains, and hence could be beneficial to consider discrepancy for these summaries separately. In this example, discrepancies are not scaled, yet MAD scaling (see algorithm T is applied to the four data summaries in simulation. We consider approximate posterior densities obtained kernel density estimates from an importance sampling with size 20000, and compare MOBOLFI approximate posterior densities with those obtained from BOLFI. The choice of auxiliary density for importance sampling has some flexibility, for results here we use a multivariate student-t distribution, with location middle points of the uniform prior and degree of freedom 3. No rejection criteria are used.
Figure 16 shows approximate BOLFI and MOBOLFI posteriors. In Figure 16a the MOBOLFI approxi- mate posterior exhibits similar variance compared to BOLFI. Figure 16b compares approximate posteriors trained with different numbers of BQ iterations. With more training iterations, the approximate marginal posterior densities stabilize. Figure 16c shows approximate posteriors with different tolerance levels, using 250 BO iterations. Inference with a 1% quantile tolerance appears superior.
38

---

Published in Transactions on Machine Learning Research (06/2025)
0.40
MOROLE
MOROUA
MOROLL
10
0.30
2.0
0.20
0.15
1.0
0.10
0.5
0.05
0.00
0.0
0.5
1.0
1.5
2.0
0.0
0.
0.6
0.6
0.8
1.0
(a) Joint Likelihood
0.36
0.30
200
2.5
250
150
10
250
190
Brue
Dus
0.5
2.0
1.5
0:15
1.0
0.10
0.5
0.05
0.00
0.0
(b) Varying the number of BO iterations
0:01 quarri
0:01 quaria
12
25
10
1.5
0.5
0.0
10
0.0
0.5
1.0
1.5
2.0
0.0
02
0.4
0.6
0.8
1.0
0
(c) Varying tolerances
Figure 16: Approximate posteriors for the day center example: (a) estimated posteriors using the joint likelihood, and the blue/orange curves show MOBOLFI/BOLFI approximate postcrior; (b) estimated posteriors using varying number of BO iterations, and the blue/ flue/oraage// grren curves show MOBOLFI ap- proximate posteriors with 250/150/55 BO iterations; (c) estimated posteriors for different tolerances and the blue/ /or//oraggggggaee green curves compare MOBOLFI approximate posteriors using 1%/5%/20 /20% tolerance levels. The dashed blue lines show the location of the true parameter values. The approximate posteriors kernel density estimates based importance sampling.
39


Reviews:
[{"review_id": "hQfxZH1T5B", "reviewer": "Reviewer_LEgz", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "S1. The core contribution is clear.\n\nS2. The paper is self-sufficient thanks to inclusion of thorough background about BOLFI. \n\nW1. The core contribution is somewhat incremental (but that doesn't preclude it from being accepted to TMLR).\n\nW2. It is not entirely clear to me yet how frequently do we expect to see the advantage of MOBOLFI over BOLFI in practice - see RC1 below.\n\nW3. Clarity of the writing could be improved here and there.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors extend BOLFI, an algorithm for computationally-efficient algorithm for Bayesian Optimization for Likelihood-Free-Inference, to use Multi-Objective optimization in the settings where the data-generating-process contains observed-random-variables which span different modalities or physical quantities, and hence would be difficult to capture model using a single likelihood function. The authors empirically showcase that MOBOLFI can outperform BOLFI on several data-generating-processes (as judged by divergence of the estimate from the true-posterior), while not requiring the tuning of \n\nThe secondary contribution is that authors show that the separate \"likelihood-approximations\" for individual modalities can be used to check the consistency between them and understand their individual impact on posterior estimation.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "RC1. As a reader, one question I was asking myself while reading the paper was: to what extent is the advantage of MOBOLFI vs BOLFI you're claiming dependent on choices of weights $w$? How sensitive is the performance of BOLFI to the choices of weights? How easy is it to choose pathologically bad weights? What's the performance difference (e.g. divergence between true-posterior and the BOLFI-estimate) between \"optimal\" and best-guess choice of weights? I think this is somewhat answered by Sec A.3 (lower half of page 25) and Fig 7, but it seems to me that these considerations are crucial to contextualize the advantage of MOBOLFI over BOLFI, and as such I would like to see these experiments moved to the main text, and the discussion extended. I think that's more important in the main paper than the details of the data-generating-processes, which could be moved to the Appendix instead.\n\nRC2. Re the MOBOLFI vs true-posterior discrepancy in Fig 5: could you comment more on this discrepancy, please? What is this discrepancy caused by? Would this discrepancy go away for higher computational budget (e.g. more evaluations)? Is there any guarantee that the LFI inference of this flavor converges under large computational budget?\n\nRC3. Are there any setting where you would expect BOLFI to perform better than MOBOLFI?\n\nRC4. Fig 5: would it be possible to include the performance of BOLFI on this problem as well?\n\nRC5. suggestion: Maybe try including the information from the first paragraph of Sec 3.3 into the Sec 1 Introduction? As a reader, I was missing the understanding of why we need the multiobjective optimization until I got to Sec 3.3.\n\nMinor:\n- Sec 1: \"the likelihood factorizes into terms for different data sources\". Given this paper treats about likelihood-free-inference, this sentence sounds imprecise.\n- typo: p(U_obs|T_obs|θ)\n- typo: \"when ∆ is closed to 0\"\n- Fig 2a: the color of \"prior\" and \"real post\" is inverted in this figure\n- Please use vector graphics. Overall, the legends on many of the figures could use more care.\n- duplicate citations: Daulton, 2021a & 2021b\n- Fig 6: more informative caption, in particular, indicate what is the difference between the columns. I think the caption should be standalone, without having to refer to the main text.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MhltidQw50", "reviewer": "Reviewer_S7en", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The work provides a fairly thorough analysis of the proposed method, showing both illustrated toy example and more applied examples.\n\nWhere the work struggles is on novelty and significance. It seems relatively straightforward to split apart a likelihood term and corresponding similarity metric. Where there might be challenges is trading off against these terms and deciding which to give more weight. But that isn't really discussed. It is mentioned that this can be used to figure out an unreliable data source, but an ablation study would do the same thing in outcome. It might even end up being the same method if likelihood terms are simply removed.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a generalisation of BOLFI to what is called the multi-objective setting. Namely one where the likelihood portion of the model can be decomposed into multiple independent terms. The paper shows that is possible to have multiple discrepancy measures and it's value when modelling domains with multiple data sources. Several illustrated examples are shown.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "It would be nice to really highlight the novelty and significance of this work. I don't understand what makes this different than an ablation study where we might remove portions of a complicated likelihood term to diagnose issues with it.", "broader_impact_concerns": "No ethical concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "hUAIoqLtWq", "reviewer": "Reviewer_ZGWA", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n* MOBOLFI addresses an important limitation of BOLFI methods, and seems to perform empirically better than BOLFI.\n* Checking for consistencies in the data source is a nice feature of MOBOLFI, that could prove useful in adversarial environments.\n* The experiments properly illustrate the contributions. \n* Detailed background on the considered frameworks is provided.\n\n**Weaknesses**\n* The paper introduces a lot of notations and this makes the reading difficult.\n* Using BO for minimizing a scalarization of discrepancies comes with asymptotic guarantees that the minimal argument of the scalarization will eventually be found. I'm not sure that the same kind of guarantees are found in MOBO, regarding the asymptotic recoverability of the Pareto set.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces MOBOLFI, a method that exploits Multi-Objective Bayesian Optimization for Likelihood-Free Inference. Broadly speaking, the method is developed for LFI tasks, based on synthetic data generated by a simulator (whose parameters need to be tuned). The simulator is assumed computationally expensive and difficult to model formally, which justifies the use of a sample-efficient black-box optimization framework such as Bayesian Optimization (BO) to tune its parameters. \n\nExisting methods (e.g., BOLFI) were designed to handle a unique data source. To use them on a problem with multiple data sources (each coming with its own measure of discrepancy), one usually needs to combine the multiple discrepancies into one before using BO to minimize this scalarization.\n\nThe main contribution of this paper is to design a method (MOBOLFI) that make use of Multi-Objective Bayesian Optimization (MOBO) able to account for multiple data sources and discrepancies. The advantages of this method are twofold: (i) an arbitrary decision is removed since MOBOLFI does not require any scalarization of the discrepancies and (ii) MOBOLFI can account for conflicting information provided by the data sources.\n\nBoth (i) and (ii) are illustrated with empirical results on synthetic and real-world experiments.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Here are my requested changes, labeled for future reference.\n\n**Critical for acceptance**\n\n1. Could you try to simplify your notations to make the reading easier? If that is not possible, could you introduce a table of notations that the reader can refer to?\n\n**Would strengthen the work**\n\n2. I find the references to equations a bit unusual. Could you refer to equations using one of the following, well-known formats: \"According to (2)\" or \"According to Equation (2)\"?\n\n3. Do you know if Multi-Objective Bayesian Optimization provides asymptotic guarantees for the recoverability of the Pareto set? If so, that is a nice additional feature of MOBOLFI. If not, that may be a drawback of using MOBOLFI if a good scalarization of discrepancies is otherwise available. Regardless of the answer to the previous point, can you add a few sentence to discuss this in the paper?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Two reviewers indicated that the work was incremental and contributions modest. However, the claims are supported and no flaws were identified. The reviewers were satisfied and recommended acceptance under the TMLR criteria. No open questions remained after the rebuttal phase, apart from further polishing the clarity of the paper.
Decision: Accept (poster)