Paper deep dive
Full Bayesian Reinforcement Learning via LF-IBIS
Stefano Masini, Cecilia Viscardi, Michela Baccini
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/5/2026, 11:14:43 AM
Summary
The paper proposes Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), a novel algorithm for Full Bayesian Reinforcement Learning (fBRL). Unlike traditional BRL methods that require an explicit or tractable likelihood function, LF-IBIS utilizes Approximate Bayesian Computation (ABC) to enable inference in environments where dynamics are described by complex simulators or black-box processes. The method provides approximate posterior distributions over both environment parameters and optimal policies, allowing for the quantification of policy uncertainty. This enables a principled approach to the exploration-exploitation trade-off, specifically through a Bayesian Mean Squared Error (BMSE) based policy updating rule. The algorithm is validated using a simulation study on response-adaptive randomization in clinical trials.
Entities (5)
Relation Signals (4)
LF-IBIS → combines → Approximate Bayesian Computation
confidence 100% · By combining Approximate Bayesian Computation with Iterated Batch Importance Sampling, LF-IBIS enables full Bayesian inference
LF-IBIS → isatypeof → Bayesian Reinforcement Learning
confidence 100% · We propose Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), a novel algorithm for BRL
LF-IBIS → addresses → Response-Adaptive Randomization
confidence 90% · We test the method on a simulation study in response-adaptive randomization in clinical trials
Maximum Entropy Reinforcement Learning → providesbasisfor → LF-IBIS
confidence 80% · the Bellman equation in (1) is not the classical RL formulation; rather, it derives from the Maximum Entropy (MaxEnt) RL framework
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Reinforcement Learning (RL) is a sequential decision-making framework in which an agent learns optimal policies through interaction with an environment by maximizing cumulative rewards. Among RL methods, Bayesian Reinforcement Learning (BRL) addresses common practical challenges related to data scarcity by leveraging prior knowledge about the environment and sequential belief updates. However, most BRL approaches require an explicit likelihood function, which is frequently inaccessible or intractable in real-world settings. We propose Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), a novel algorithm for BRL that updates the agent's beliefs online as new interactions become available. By combining Approximate Bayesian Computation with Iterated Batch Importance Sampling, LF-IBIS enables full Bayesian inference in settings where the environment dynamics are not described by an explicit or tractable likelihood. The method yields approximate posterior distributions over both environment parameters and optimal policies, providing a quantification of policy uncertainty useful for a Bayesian treatment of the exploration-exploitation trade-off. We test the method on a simulation study in response-adaptive randomization in clinical trials, where closed-form posteriors enable validation. Additional experiments address settings where the posterior has no closed form and illustrate online policy updating based on the posterior distribution of the optimal policy.
Tags
Links
- Source: https://arxiv.org/abs/2607.01741v1
- Canonical: https://arxiv.org/abs/2607.01741v1
Trouble viewing inline? Open PDF directly →
Full Text
115,689 characters extracted from source content.
Expand or collapse full text
Full Bayesian Reinforcement Learning via LF-IBIS Stefano Masini Department of Computer Science, University of Pisa, Italy. Email:stefano.masini@phd.unipi.it Cecilia Viscardi Department of Economics and Statistics, University of Salerno, Italy. Email: cviscardi@unisa.it Michela Baccini Department of Statistics, Computer Science, Applications ”G.Parenti”, University of Florence, Italy. Email: michela.baccini@unifi.it Abstract Reinforcement Learning (RL) is a sequential decision-making framework in which an agent learns optimal policies through interaction with an environment by maximizing cumulative rewards. Among RL methods, Bayesian Reinforcement Learning (BRL) addresses common practical challenges related to data scarcity by leveraging prior knowledge about the environment and sequential belief updates. However, most BRL approaches require an explicit likelihood function, which is frequently inaccessible or intractable in real-world settings. We propose Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), a novel algorithm for BRL that updates the agent’s beliefs online as new interactions become available. By combining Approximate Bayesian Computation with Iterated Batch Importance Sampling, LF-IBIS enables full Bayesian inference also in settings where the environment dynamics is not described by an explicit or tractable likelihood. The method yields approximate posterior distributions over both environment parameters and optimal policies, providing a quantification of policy uncertainty useful for a Bayesian treatment of the exploration–exploitation trade-off. We test the method on a simulation study in response-adaptive randomization in clinical trials, where closed-form posteriors enable validation. Additional experiments address settings where the posterior has no closed form and illustrate online policy updating based on the posterior distribution of the optimal policy. Keywords: Likelihood-free inference, Importance Sampling, Bayesian Reinforcement Learning, Sequential Learning, Response-Adaptive Randomization 1 Introduction Reinforcement Learning (RL) is a computational framework for sequential decision-making in which an agent interacts with an environment to learn how to choose actions that maximize cumulative rewards over time [33]. At each step, the agent observes the current state, selects an action, and receives a reward or penalty along with a successor state, gradually improving its behaviour through trial and error. RL has become widely used in several domains—including robotics, healthcare and online decision systems—thanks to its ability to autonomously discover effective strategies in complex and uncertain settings. RL methods are commonly classified into model-free and model-based approaches. Model-based RL learns a model of the environment behaviour —often formalized as a Markov Decision Process (MDP)— that specifies the state-transition dynamics [39]. The agent learns this model from interaction data and uses it to plan and compute its policies. In contrast, model-free methods learn policies (or value functions) directly from observed state–action–reward sequences. In this work, we focus on the model-based formulation. In classical RL, effective learning usually requires a large amount of data collected through many interactions with the environment [33]. In many real-world applications, however, acquiring such interaction data is expensive, limited, or risky. This exacerbates the well-known exploration–exploitation trade-off: the choice between exploring new actions and gather better information, or exploiting actions currently believed to be optimal to maximize immediate rewards [33, 1]. Bayesian Reinforcement Learning (BRL) addresses this challenge by integrating the principles of Bayesian inference with RL, enabling the incorporation of prior knowledge and a continual updating of the belief as new data become available. This leads to more data-efficient learning when interactions are limited [17, 8]. A major obstacle of any approach grounded in a statistical formulation of the RL problem, including BRL, is the requirement for an explicit likelihood function associated with the model that describes the environment’s behavior. In many real-world scenarios, such a likelihood is intractable or unavailable, for instance when the environment is modeled as a complex simulator or a black-box process [18]. To mitigate this limitation, we explore the use of Approximate Bayesian Computation (ABC), a class of likelihood-free inference methods that approximate the intractable likelihood through comparison between observed data and data simulated from the complex model [30]. Although ABC has been successfully applied in various domains, its use in sequential decision-making and BRL is still limited. A few studies have explored the connection between ABC and RL. Ritto et al. [25] combine RL-based model selection with ABC parameter estimation in a digital-twin setting, using Thompson sampling [36] (via a Beta-Bernoulli bandit) for model choice and ABC for updating the parameter posterior. Dimitrakakis and Tziortziotis [13] introduce ABC into RL as a method to perform inference in environments where the transition model is unknown or the likelihood is intractable. In their framework, the agent maintains a prior distribution over environment models and updates this distribution by sampling from it, simulating trajectories, and retaining only those samples whose simulated behavior closely matches observed data, according to a predefined distance metric. This likelihood-free approach enables the construction of approximate posterior distributions over environment models, which in turn supports principled exploration through model uncertainty. Their work represents a pioneering application of ABC to RL, demonstrating its potential to handle black-box environments without requiring access to differentiable or analytically defined transition dynamics. Our work extends this line of research in two key directions: (i) Online inference with progressive histories: unlike the offline-oriented setting of [13], our approach is designed for online RL and incrementally updates posterior beliefs as new interactions with the environment occur. (i) Joint posterior distributions over models and policies: while the original ABC-RL approach focuses on computing posterior distributions only over environment models, we propagate uncertainty over candidate stochastic policies, allowing the agent to exploit the whole posterior belief. We propose adapting policies by sampling from their posterior distributions, leading to a fully Bayesian analogue of Thompson sampling, in which exploration is driven by posterior uncertainty and exploitation naturally emerges as this uncertainty decreases. To this end, we introduce a framework for full Bayesian Reinforcement Learning (fBRL) based on a novel algorithm, Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), which combines Sequential Monte Carlo ABC (SMC-ABC) principles [30] with Iterated Batch Importance Sampling [6] in an RL setting. A similar approach is proposed by Roy et al. [28], who introduce a generalized Bayesian framework for deep reinforcement learning. Our approach shares with Roy et al.’s framework the fully Bayesian perspective, both the methods rely on SMC for posterior approximation, and exploit multiple posterior samples, extending classical Thompson sampling. However, their method is based on generalized Bayesian updating via prequential scoring rules, which act as a surrogate likelihood, whereas our approach adopts a fully likelihood-free framework based on ABC and direct simulation. Moreover, their framework leverages deep generative models and expected Thompson sampling, while our method exploits full posterior information and operates in a fully online, sequential setting as new data arrive. We complement our approach with an analysis of the asymptotic properties and validate it empirically in a Response-Adaptive Randomization problem for clinical trials [38], where a conjugate Bayesian solution allows direct comparison with LF-IBIS, standard rejection ABC, and the exact posterior. Structure of the paper Section 2 introduces RL, its statistical formulation and the extension to the full Bayesian setting. Section 3 presents the likelihood-free formulation and details our LF-IBIS algorithm, highlighting how adaptive ABC thresholding mechanism of SMC-ABC [10] is integrated with online inference in an IBIS paradigm[6]. Theoretical insights are provided in Section 3.3.4. In Section 4, we test the method empirically in the context of Response-Adaptive Randomization for clinical trials. Section 5 concludes the paper with a summary and discussion. Supplementary material contains proofs and additional details about algorithms and experimental settings. 2 Full Bayesian Reinforcement Learning In this work, we aim to define a procedure for providing a full Bayesian inference methodology in RL. To this end, we first cast RL into a statistical framework and define a likelihood function associated with the model assumed for the environment’s behavior. Then, we exploit the statistical formulation to provide Bayesian inference for the quantities of interest. Here, the term full Bayesian refers to an approach in which posterior uncertainty is not only represented at the level of the environment model, value functions, policy parameters, or rewards, but is propagated to the optimal policy. In contrast with Bayesian RL methods that mainly exploit posterior expected quantities, posterior summaries, or single posterior samples, our objective is to approximate the posterior distribution over the optimal policy itself. In what follows, capital letters denote random variables, while lowercase letters their realizations. Calligraphic letters denote sets and variables’ domains. For a finite set Y |||Y| denotes the cardinality and Δ() (Y) denotes the probability simplex over Y. ℕN, ℝR and ℝ+R^+ denote the sets of natural, real and non-negative real numbers, respectively. We use Pr(⋅) (·) to refer to probabilities of events, p(⋅)p(·) to generically denote probability density or mass functions, P to denote stochastic matrices. 2.1 Statistical Reinforcement Learning In RL, an agent learns a strategy to make optimal decisions while interacting with an environment and receiving rewards based on its actions. The interaction unfolds as a sequence of action–state–reward triples, the so-called observed history x=(an,sn,rn)n=1Nx=\(a_n,s_n,r_n)\_n=1^N, with N equal to the number of agent-environment interactions.111The history can alternatively be indexed by the time at which interactions occur. Here, we assume a single interaction per time step. From a statistical point of view, action ana_n and state sns_n can be seen as realizations of two random variables, An∈nA_n _n and Sn∈nS_n _n, while rnr_n, shorthand for r(an,sn)r(a_n,s_n), is the reward resulting from a (either deterministic or stochastic) function rn:n×n→ℛ⊆ℝ.r_n:A_n×S_n . Throughout, the state space nS_n, the action space nA_n, and the reward space ℛR are assumed to be finite sets. A common assumption is that the dynamic of the system can be described by a homogeneous Markov Decision Process (MDP) [3, 39], meaning that the following two assumptions hold: ASS.1 Sn+1⟂aj,sjj=1n−1∣(an,sn)S_n+1 \!\!\! \a_j,s_j\_j=1^n-1 (a_n,s_n) for each n∈1,…N−1n∈\1,...N-1\. This establishes that the sequence of states forms a first order Markov chain and that the environment’s behavior is fully described through |n||A_n| transition matrices denoted by Pn(a)P_n(a), with entries Pr(s′∣s,a):=Pr(Sn+1=s′∣Sn=s,An=a) (s s,a):= (S_n+1=s S_n=s,A_n=a) for (s,s′)∈n2(s,s ) _n^2 and a∈na _n. ASS.2 n≡,n≡,Pn(a)≡P(a)S_n ,A_n ,P_n(a)≡ P(a) for each n∈1,…,Nn∈\1,...,N\ and a∈a , where S and A are two finite sets, each P(a)P(a) is a (kS×kS)(k_S× k_S)-stochastic matrix, with kS:=||k_S:=|S|. This establishes the homogeneity of the Markov process. Note that the observed history is a realization of a multivariate discrete random variable X that takes values in NX^N, where :=×ℛX:=A×S×R. Its joint probability distribution depends on 1) kA:=||k_A:=|A| stochastic matrices, P(a)P(a), with a∈a ; 2) a stochastic policy π:→Δ()π:S→ (A) where π(a∣s)π(a s) is the probability of choosing action a from state s. 222A stochastic policy is a mapping π:→Δ()π:S→ (A) , whereas a deterministic policy is a mapping π:→π:S . We restrict attention to stochastic policies. Given a policy π(⋅∣⋅)π(· ·), the value function Vπ(s)V_π(s) quantifies the expected future reward that can be gained from state s. Assuming the transition matrices P(a)P(a) as known, it can be computed recursively by the Bellman equation [3] Vπ(s)=∑a∈π(a∣s)Qπ(s,a)+λℋ(π(⋅∣s))∀s∈V_π(s)= _a π(a s)Q_π(s,a)+ (π(· s)) ∀ s (1) where ℋ(π(⋅∣s))H(π(· s)) is the Shannon entropy of the stochastic policy ℋ(π(⋅∣s)):=−∑a∈π(a∣s)logπ(a∣s)H(π(· s)):=- _a π(a s) π(a s) and Qπ(s,a)Q_π(s,a) denotes state–action value function Qπ(s,a):=∑s′∈Pr(s′∣s,a)[r(a,s′)+γVπ(s′)].Q_π(s,a):= _s (s s,a)[r(a,s )+γ V_π(s )]. Here γ∈[0,1]γ∈[0,1] is the discount factor, which determines how future rewards are weighted relative to immediate ones. Note that the Bellman equation in (1) is not the classical RL formulation; rather, it derives from the Maximum Entropy (MaxEnt) RL framework [40]. Unlike standard RL—where optimizing the expected total reward often leads to deterministic policies—the MaxEnt objective includes an entropy term that explicitly promotes stochasticity. The parameter λ∈ℝ+λ ^+ controls the relative weight of the entropy contribution, thereby regulating the degree of randomness in the resulting policy. Larger values of λ promote higher-entropy policies and thus greater randomness. As noted in [20], the temperature parameter (typically denoted by α in the reinforcement learning literature) can be absorbed into the reward function through a simple rescaling by λ−1λ^-1. Accordingly, it does not affect the structure of the optimal policy, and we omit it in the following by fixing λ=1λ=1 without loss of generality. An optimal policy π∗(a∣s)π^*(a s) is any policy that allows the agent to reach the optimal value function V∗(s)=maxπVπ(s)V^*(s)= _πV_π(s) for every state. Under the MaxEnt objective, several classical RL algorithms—such as Policy Iteration, Value Iteration, and Q-Iteration—can be used to compute π∗π^* [33]. All these methods rely on access to the transition model P, or on accurate estimates of either its entries or the corresponding value functions. In many cases, empirical estimates of transition probabilities can be obtained through Monte Carlo sampling [33]. In what follows, we cast RL within a statistical inference framework that enables learning of both the transition dynamics and the optimal policy. Likelihood function in RL Recall that under the MDP assumptions, the dynamic of the environment is fully described by kAk_A matrices denoted by P(a)P(a) with a∈a . These are (kS×kS)(k_S× k_S)-stochastic matrices, where each row corresponds to a state s and contains the conditional probability distribution over S, denoted by Pr(⋅∣s,a) (· s,a). Let us denote by μ∈ℳμ , with ℳ⊆ℝdM ^d a vector of parameters of size d for some d∈ℕd . We assume that the rows of the matrix P(a)P(a) are unknown elements of a family of conditional probability distributions parametrized by μ: ℱ:=Pr(⋅∣s,a,μ):μ∈ℳ.F:=\ (· s,a,μ):μ \. Under this assumption, given a policy π, the random history X has a distribution represented by the following probability function: Pr(X=x∣μ,π)=ρ(s1)∏n=1Nπ(an∣sn)Pr(Sn+1=sn+1∣sn,an,μ)x∈N, (X=x μ,π)=ρ(s_1)Π _n=1^Nπ(a_n s_n) (S_n+1=s_n+1 s_n,a_n,μ) x ^N, (2) where ρ(⋅)ρ(·) is the initial state distribution. This probability has key relevance when addressing inferential purposes since it is proportional to the likelihood function for the parameters μ that should be inferred to learn the environment transition matrix. In particular, let us define f(x∣μ)≔∏n=1NPr(Sn+1=sn+1∣sn,an,μ).f(x μ) Π _n=1^N (S_n+1=s_n+1 s_n,a_n,μ). Given the initial state s1s_1 and a policy π(⋅∣⋅)π(· ·), from Eq (2) follows that ℒ(μ;x,π,s1)∝f(x∣μ)L(μ;x,π,s_1) f(x μ) (3) where ℒ(μ;x,π,s1)L(μ;x,π,s_1) denotes the likelihood function, which will henceforth be written as ℒ(μ;x)L(μ;x) for brevity. This statistical formulation lays the foundation for inferential procedures both in a frequentist and Bayesian framework where an agent wants to learn the behavior of an environment governed by μ, based on the sequence of its states s1,…,sNs_1,…,s_N in an observed history x. In this work, we focus on a Bayesian treatment of the RL problem, which is presented in the following section. 2.2 Bayesian Reinforcement Learning Bayesian approaches to RL may represent posterior uncertainty over different quantities, including environment models, value functions, policy parameters, policy gradients, or reward functions. In the framework proposed in this thesis, inference is performed over the parameters μ governing the environment dynamics, in line with model-based BRL approaches such as BAMDP-based formulations [15] and posterior sampling methods [23, 31]. In this context, the parameters μ are treated as random variables. The agent starts with a prior belief expressed by the prior distribution p(μ)p(μ) over ℳM. Then, it uses the likelihood function in equation (3) to update it and compute the posterior p(μ∣x)p(μ x) through the Bayes formula: p(μ∣x)∝p(μ)ℒ(μ;x).p(μ x) p(μ)L(μ;x). (4) In this work, we adopt a full Bayesian perspective. Unlike the model-based BRL approaches the posterior is not used solely to identify a specific parameter vector—e.g., the maximum a posteriori estimate or posterior expected quantities—to evaluate the value function and determine an optimal policy, but all the information contained in the posterior distribution in Eq (4) is exploited. In fact, the optimal policy in turn depends on the vector of parameters μ through the transition probabilities. We make this dependence explicit in the value functions and the optimal policy. The optimal state–action value function Qμ∗(s,a)Q^*_μ(s,a) for a given μ is Qμ∗(s,a)=∑s′∈Pr(s′∣s,a,μ)[r(a,s′)+γVμ∗(s′)],Q^*_μ(s,a)= _s (s s,a,μ) [r(a,s )+γ\,V^*_μ(s ) ], where Vμ∗(s)V^*_μ(s) is the optimal (entropy-regularized) value function retrieved by solving the soft Bellman fixed-point Vμ∗(s)=λlog(∑a∈exp(1λQμ∗(s,a))).V^*_μ(s)=λ \! ( _a \! ( 1λQ^*_μ(s,a) ) ). This value function corresponds to an optimal policy πμ∗(a∣s)=exp(1λQμ∗(s,a))∑a′∈exp(1λQμ∗(s,a′)).π^*_μ(a s)= \! ( 1λQ^*_μ(s,a) )Σ _a \! ( 1λQ^*_μ(s,a ) ). (5) The soft Bellman fixed–point equation and the corresponding Boltzmann optimal policy are standard results in maximum–entropy reinforcement learning; an explicit derivation can be found in [19], building on the maximum-entropy formulation of [41]. We exploit their results to provide a posterior distribution for the optimal policy, rather than a point estimate. In fact, if we denote by πμ∗π^*_μ the stochastic matrix of probabilities πμ∗(a∣s):(a,s)∈×\π^*_μ(a s):(a,s) ×S\, then it follows from equation (5) that πμ∗π^*_μ is a function of μ. To compute πμ∗π^*_μ given a specific μ value, a standard optimal-policy computation method can be employed; in our case, we adopt Policy Iteration. Since each μ value is associated with a posterior density (or probability), the uncertainty about μ can be propagated to πμ∗π^*_μ, thus computing p(π∗∣x)p(π^* x). Exploiting the posterior distribution of the optimal policy requires a simulation method that allows us to repeatedly sample from it. To this end, Section 3 introduces the proposed algorithm for approximating the posterior distribution of π⋆π . Once this distribution has been approximated, it can be used to address the exploration–exploitation trade-off in an uncertainty-aware way. A possible decision strategy based on this posterior information is then presented in Section 4.3. 2.2.1 Posterior-Based Policy Updating The posterior distribution over policies provides valuable information for deciding when it is advantageous to update the current policy. To this end, we adopt a sampling-based decision rule grounded in the Bayesian Mean Squared Error (BMSE). Let π(ℓ)ℓ=1L\π^( )\_ =1^L denote a collection of policies sampled from the posterior distribution of the optimal policy. The BMSE with respect to a given policy π is defined as follows: BMSE(π)=1L∑ℓ=1L|π(ℓ)−π|2.BMSE(π)= 1L _ =1^L|π^( )-π|^2. (6) In the implementation of an RL procedure, this quantity can be used to assess how well a given policy aligns with the posterior distribution of optimal policies. In particular, at each update step, we compute BMSE(π)BMSE(π) for the current policy and BMSE(π~)BMSE( π) for a new candidate policy π~ π drawn from the posterior distribution. The policy with the lower BMSE is then selected. If the selected policy coincides with the current policy, we interpret this as evidence that the posterior mass is already concentrated around it, and thus further exploration of the environment may be beneficial. Conversely, if the sampled policy yields a lower BMSE, we update the current policy accordingly, exploiting the information contained in the posterior. This procedure naturally balances exploration and exploitation: the posterior distribution encodes uncertainty over optimal policies, while the BMSE-based criterion provides a principled mechanism to decide whether to retain or update the current policy. 3 Likelihood-free Iterated Batch Importance Sampling (LF-IBIS) The posterior distribution in Eq. (4) is often not available in closed form. As a consequence, the posterior distribution over the optimal policy is not available analytically either. As in Bayesian analysis in general, posterior computations often require methods for getting samples from the posterior distribution such as Monte Carlo or Markov Chain Monte Carlo algorithms [26]. In the presented RL framework, two issues make this task particularly challenging: 1. Intractable or unavailable likelihood. The likelihood ℒ(μ;x)L(μ;x) depends on the transition matrices parameterized by μ, whose entries may be complex or unknown functions of the parameters. In other words, the family ℱF may be unknown, or it may consist of non-computable elements. This occurs because environment dynamics are often described by computational models that can be simulated but lack analytical form. 2. Sequential data. The history x is observed gradually over time, and recomputing posterior distributions as new interactions agent-environment occur can be computationally prohibitive. Likelihood–based Monte Carlo methods are therefore not directly applicable. To address these issues, we propose a hybrid algorithm, Likelihood-free Iterated Batch Importance Sampling, which combines Approximate Bayesian Computation (ABC) with Iterated Batch Importance Sampling (IBIS). 3.1 ABC in Reinforcement Learning Approximate Bayesian Computation (ABC) provides a likelihood-free approach to Bayesian inference when likelihood evaluation is infeasible or analytically intractable [29, 32]. Instead of computing the exact likelihood, ABC relies on a Monte Carlo approximation based on the simulation of pseudo-data, i.e. simulated histories, from the generative model associated with the environment – see [30] for details about ABC and [13] for a first application to RL. The simplest ABC sampling scheme [35, 24], i.e., the standard Rejection ABC, proceeds by sampling parameters from the prior, using these parameters to generate pseudo-data from a simulator that reproduces the generative process underlying the intractable likelihood, computing distances between simulated and observed data, and accepting parameters for which the distance is within a tolerance threshold ϵ>0ε>0. The output is a sample from an approximate posterior distribution, with the quality of approximation depending on the choice of ϵε. Let y∈Ny ^N denote a simulated history, sampled from a simulator that reproduces the environment behaviour. The output of the ABC procedure is a sample from the following approximate joint posterior distribution: pϵ(μ,y∣x)=p(μ)ℒ(μ;y)Aϵ,x(y)∫ℳ∑y∈Aϵ,xp(μ)ℒ(μ;y)dμ∝p(μ)f(y∣μ)Aϵ,x(y)p_ε(μ,y x)= p(μ)L(μ;y)I_A_ε,x(y) _MΣ _y∈ A_ε,xp(μ)L(μ;y)dμ p(μ)f(y μ)I_A_ε,x(y) . Here, Z(⋅)I_Z(·) is the indicator function of a set Z, the set Aϵ,x=y∈N:d(y,x)≤ϵA_ε,x=\y ^N:d(y,x)≤ε\ identifies the acceptance region, and d:N×N→ℝ+d:X^N×X^N ^+ is a distance function. By marginalizing out y, we obtain the approximate posterior over the parameters μ: pϵ(μ∣x)∝p(μ)fϵ(x∣μ)p_ε(μ x) p(μ)f_ε(x μ) where: fϵ(x∣μ)≔Pr(d(y,x)≤ϵ∣μ)=∑y∈Nf(y∣μ)Aϵ,x(y).f_ε(x μ) (d(y,x)≤ε μ )=Σ _y ^Nf(y μ)I_A_ε,x(y). (7) This quantity represents the ABC likelihood, that can be estimated via Monte Carlo sampling. Specifically, one can draw M i.i.d. pseudo-dataset y1,…,yM∼f(⋅∣μ)\y_1,…,y_M\ f(· μ) and approximate the probability in Eq (7) by fϵ,M(x∣μ)≔1M∑m=1MAϵ,x(ym).f_ε,M(x μ) 1MΣ _m=1^MI_A_ε,x(y_m). (8) Note that the probability of drawing pseudo-data close to the observed data, Pr(d(x,y)≤ϵ∣μ) (d(x,y)≤ε μ ), depends also on N. Thus, to make the algorithm more scalable and get a good Monte Carlo approximation without requiring huge M, a common practice is to summarize both observed and simulated data through lower-dimensional summary statistics η(⋅)η(·). This implies a new definition of the acceptance region as Aϵ,η(x)=y∈N:d(η(y),η(x))≤ϵA_ε,η(x)=\y ^N:d(η(y),η(x))≤ε\, and the introduction of a new source of approximation. Moreover, the indicator function is often replaced by a smooth kernel function κ(⋅)κ(·), which allows for a continuous scaling from 1 (when x=yx=y) to 0 (when d(η(y),η(x))>ϵd(η(y),η(x))>ε), leading to the following posterior approximation: pϵ(μ∣x)∝p(μ)fϵ(x∣μ)≈p(μ)⋅1M∑m=1Mκ(d(ym,x);ϵ)p_ε(μ x) p(μ)f_ε(x μ)≈ p(μ)· 1MΣ _m=1^Mκ (d(y_m,x);ε ) (9) Henceforth, to simplify notation, x and y will denote, according to the context, either the complete data or their corresponding summary statistics η(x)η(x) and η(y)η(y). In Eq (9) the ABC kernel κ is typically a smoothing function (e.g., Epanechnikov or Gaussian) that uses ϵε as scale parameter and is often defined over a compact support. Note that the standard Rejection ABC algorithm uses M=1M=1 and relies on a crude Monte Carlo estimate of the ABC likelihood which, for the sake of this work, is instead a key quantity, as will become clear in the remainder of the paper. In fact, more advanced alternative to Rejection ABC, such as the Sequential Monte Carlo Approximate Bayesian Computation (SMC-ABC) algorithm proposed by Del Moral et al. [10], uses the Monte Carlo approximation in Eq (8). This method adapts classical SMC methodology to the ABC framework by constructing a sequence of intermediate approximate distributions pϵi(μ∣x)i=1n\p_ _i(μ x)\_i=1^n, each associated with one of the decreasing tolerance levels ϵ1>⋯>ϵn _1>…> _n. The algorithm iterates importance sampling steps as displayed in Alg 1. At each iteration i, the algorithm gets a weighted sample of size L, μi(ℓ),ωi(ℓ)(ϵi)ℓ=1L\ _i^( ), _i^( )( _i)\_ =1^L, through three main steps: (i) reweighting particles from the previous iteration – i.e. the weighted sample from pϵi−1p_ _i-1 – according to the ABC likelihood that uses the current threshold ϵi _i, (i) resampling L particles using new weights, (i) perturbing particles through a transition kernel – e.g., MCMC kernel of invariant distribution pϵip_ _i that uses a proposal distribution q(⋅,⋅)q(·,·). In the version proposed by Del Moral et al. [10], the sequence ϵii=1n\ _i\_i=1^n is not fixed in advance but instead chosen adaptively to gradually improving the quality of the approximation while controlling the accuracy of the posterior estimate. The key idea is starting with a large tolerance ϵ1 _1, which allows more simulated dataset to be accepted, and then progressively reduce it to focus the inference on regions of higher posterior density. To adapt the threshold, the algorithm monitors the Effective Sample Size (ESS): the next threshold ϵ∗ε^* is chosen by solving the following equation through a bisection algorithm: ESS(ωi(ϵ∗))=α⋅ESS(ωi−1(ϵi−1)),ESS ( _i(ε^*) )=α·ESS ( _i-1( _i-1) ), (10) where ωi(⋅) _i(·) denote the corresponding (unnormalized) importance weights, and α∈(0,1)α∈(0,1) is a tuning parameter (e.g., α=0.9α=0.9) that controls how aggressively the tolerance is reduced. Algorithm 1 SMC-ABC Algorithm (Del Moral et al. 2012) 1:Input: target threshold ϵε, starting threshold ϵ1 _1, α∈(0,1)α∈(0,1), nr of simulations L, number of pseudo-data M, observed data x; proposal distribution q(⋅,⋅)q(·,·) 2:Draw μ1(ℓ)∼p(⋅) _1^( ) p(·) for ℓ∈1,…,L ∈\1,...,L\ 3:For each ℓ , simulate y1,m(ℓ)∼f(⋅∣μ1(ℓ))y_1,m^( ) f(· _1^( )) for m∈1,…,Mm∈\1,...,M\ 4:Compute initial weights ω1(ϵ1)=ω1(ℓ)(ϵ1)ℓ=1L _1( _1)=\ _1^( )( _1)\_ =1^L where ω1(ℓ)(ϵ1)=∑m=1Mκ(d(y1,m(ℓ),x);ϵ1)M _1^( )( _1)= _m=1^Mκ (d(y_1,m^( ),x); _1 )M 5:i=1i=1 6:while ϵi>ϵ _i>ε do 7: Set i=i+1i=i+1 8: Reweighting:: Compute new ϵi _i by solving with bisection method: ESS(ωi(ϵi))=αESS(ωi−1(ϵi−1))ESS ( _i( _i) )= ( _i-1( _i-1) ) 9: Where ωi(ℓ)(ϵi)=ωi−1(ℓ)(ϵi−1)∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi)∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi−1) _i^( )( _i)=ω^( )_i-1( _i-1) _m=1^Mκ (d(y_i-1,m^( ),x); _i ) _m=1^Mκ (d(y_i-1,m^( ),x); _i-1 ) 10: Resampling: Sample L pairs from (μi−1(ℓ),yi−1,m(ℓ))ℓ=1L\( _i-1^( ),y_i-1,m^( ))\_ =1^L using their weights ωi(ℓ)(ϵi)∑ℓ=1Lωi(ℓ)(ϵi) _i^( )( _i) _ =1^L _i^( )( _i) as probabilities. Denote new sampled pairs, with a slight abuse of notation, by (μi−1(ℓ),yi−1,m(ℓ))ℓ=1L\( _i-1^( ),y_i-1,m^( ))\_ =1^L. 11: Set new weights: ωi(ℓ)(ϵi)=1ω^( )_i( _i)=1 12: Moving: Draw μ∗∼q(μi−1(ℓ),⋅)μ^* q( _i-1^( ),·) for ℓ∈1,…,L ∈\1,...,L\ 13: For each ℓ , simulate ym∗∼f(⋅∣μ∗)y_m^* f(· μ^*) for m∈1,…,Mm∈\1,...,M\ 14: Compute the acceptance ratio: ar(μi(ℓ)→μ∗)=p(μ∗)∑m=1Mκ(d(ym∗,x);ϵi)q(μ∗,μi(ℓ))p(μi(ℓ))∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi)q(μi(ℓ),μ∗)ar( _i^( )→μ^*)= p(μ^*) _m=1^Mκ (d(y_m^*,x); _i )q(μ^*, _i^( ))p( _i^( )) _m=1^Mκ (d(y_i-1,m^( ),x); _i )q( _i^( ),μ^*) for ℓ∈1,…,L ∈\1,...,L\ 15: Draw u(ℓ)∼Unif(0,1)u^( ) (0,1) for ℓ∈1,…,L ∈\1,...,L\ 16: if u(ℓ)<min1,ar(μi(ℓ)→μ∗)u^( )< \1,ar( _i^( )→μ^*)\ then 17: Set μi(ℓ)=μ∗ _i^( )=μ^* and yi,m(ℓ)=ym∗y^( )_i,m=y^*_m for ℓ∈1,…,L ∈\1,...,L\ 18: else 19: Set μi(ℓ)=μi−1(ℓ) _i^( )= _i-1^( ) and yi,m(ℓ)=yi−1,m(ℓ)y^( )_i,m=y_i-1,m^( ) for ℓ∈1,…,L ∈\1,...,L\ 20: end if 21:end while Note that in [10], the ABC kernel is specified as an indicator function, whereas Algorithm 1 is formulated more generally to accommodate arbitrary ABC kernels. Using multiple simulations (M≥1M≥ 1) per particle and/or alternative kernel choices may improve the stability of the algorithm by reducing variability among the importance weights. Further details on this aspect are provided in Appendix A.1. More generally, the kernel function κ, the summary statistics η(⋅)η(·), the distance function, the metrics used to evaluate the efficiency of posterior estimates (e.g., ESS) etc, can be replaced by more sophisticated alternatives to better capture the relevant features of the data and improve inference quality. These aspects will be discussed specifically with reference to the case at hand in the following sections of the paper. 3.2 Iterated Batch Importance Sampling (IBIS) The Iterated Batch Importance Sampling (IBIS) algorithm, introduced by Chopin [6], is a Sequential Monte Carlo (SMC) method designed for static scenarios, where a single posterior distribution is targeted. The central idea is to build a sequential data setting to address high-dimensional problems in which standard Monte Carlo methods are impractical, due to the computational burden associated with data size and/or the lack of well-informed proposal distributions. The algorithm iterates importance sampling steps as displayed in Alg 2. At iteration i, it draws a weighted sample μi(ℓ),ωi(ℓ)ℓ=1L\ _i^( ), _i^( )\_ =1^L from the partial posterior p(μ∣x1:i)p(μ x_1:i), conditioned on observations x1:i=xjj=1ix_1:i=\x_j\_j=1^i. The algorithm proceeds through a reweighting, a resampling, and a moving step. During the reweighting phase, the weighted sample from the previous iteration, i.e., the particles from p(μ∣x1:i−1)p(μ x_1:i-1), is adapted to account for the new observation xix_i, yielding an updated approximation of p(μ∣x1:i)p(μ x_1:i).333For simplicity, we assume that only one observation is added per iteration, although the algorithm is generally formulated to handle batches larger than one. Starting from the prior, the algorithm iteratively refines the posterior by using the approximate distribution from the previous iteration as the proposal for the next iteration. Therefore, at the i-th iteration, the unnormalized importance weight for the ℓ -th particle is ωi(ℓ)=ωi−1(ℓ)p(μi−1(ℓ)∣x1:i)p(μi−1(ℓ)∣x1:i−1)=ωi−1(ℓ)p(μi−1(ℓ))f(x1:i∣μi−1(ℓ))p(μi−1(ℓ))f(x1:i−1∣μi−1(ℓ))=ωi−1(ℓ)⋅f(xi∣x1:i−1,μi−1(ℓ))ω^( )_i=ω^( )_i-1 p(μ^( )_i-1 x_1:i)p(μ^( )_i-1 x_1:i-1)=ω^( )_i-1 p(μ^( )_i-1)f(x_1:i μ^( )_i-1)p(μ^( )_i-1)f(x_1:i-1 μ^( )_i-1)=ω^( )_i-1· f(x_i x_1:i-1, _i-1^( )) (11) where f(xi∣x1:i−1,μi−1(ℓ))f(x_i x_1:i-1, _i-1^( )) is the predictive likelihood and serves as reweighting factor to account for the new data point xix_i. Algorithm 2 IBIS Algorithm (Chopin, 2002) 1:Input Prior p(μ)p(μ), number of particles L, observed data x, proposal distribution q(⋅,⋅)q(·,·) 2:Draw μ1(ℓ)∼p(μ) _1^( ) p(μ) for ℓ∈1,…,L ∈\1,...,L\ 3:Compute initial weights ω1=ω1(ℓ)ℓ=1L _1=\ _1^( )\_ =1^L where ω1(ℓ)=1 _1^( )=1 for ℓ∈1,…,L ∈\1,...,L\ 4:for each time step i=2,…,ni=2,…,n do 5: Reweighting: ωi(ℓ)=ωi−1(ℓ)⋅f(xi∣x1:i−1,μi−1(ℓ)) _i^( )= _i-1^( )· f(x_i x_1:i-1, _i-1^( )) for ℓ∈1,…,L ∈\1,...,L\ 6: Resampling: Sample L particles μi−1(ℓ),ωi(ℓ)ℓ=1L\ _i-1^( ), _i^( )\_ =1^L using weights ωi(ℓ)∑ℓ=1Lωi(ℓ) _i^( ) _ =1^L _i^( ) as probabilities 7: Set new weights: ωi(ℓ)=1ω^( )_i=1 8: Moving: Draw μ∗∼q(μi−1(ℓ),⋅)μ^* q( _i-1^( ),·) for ℓ∈1,…,L ∈\1,...,L\ 9: Compute the acceptance ratio: ar(μi−1(ℓ)→μ∗)=p(μ∗)q(μ∗,μi−1(ℓ))p(μi−1(ℓ))q(μi−1(ℓ),μ∗)ar( _i-1^( )→μ^*)= p(μ^*)q(μ^*, _i-1^( ))p( _i-1^( ))q( _i-1^( ),μ^*) for ℓ∈1,…,L ∈\1,...,L\ 10: Draw u(ℓ)∼Unif(0,1)u^( ) (0,1) for ℓ∈1,…,L ∈\1,...,L\ 11: if u(ℓ)<min1,ar(μi−1(ℓ)→μ∗)u^( )< \1,ar( _i-1^( )→μ^*)\ then 12: Set μi(ℓ)=μ∗ _i^( )=μ^* for ℓ∈1,…,L ∈\1,...,L\ 13: else 14: Set μi(ℓ)=μi−1(ℓ) _i^( )= _i-1^( ) for ℓ∈1,…,L ∈\1,...,L\ 15: end if 16:end for This sampling approach is particularly well adapted to situations in which data become gradually available rather than being accessible from the beginning, like in RL. However, in our framework, pointwise evaluations of the likelihood function are infeasible; therefore, it cannot be directly applied. The next section presents a new algorithm that extends its applicability to reinforcement learning. 3.3 LF-IBIS In RL, at each time step, the agent observes the outcome of its action and uses this information to improve its understanding of the environment. This means that we an algorithm that updates the posterior for μ, adapting a sample from p(μ∣x1:i)p(μ x_1:i) to a sample from p(μ∣x1:i+1)p(μ x_1:i+1), as new data xi+1x_i+1 become available. To this end, our proposal is an algorithm that alternates two key phases: • IBIS steps, where the observed history is extended by including the new interaction; • SMC-ABC steps, where the discrepancy threshold ϵi _i is adaptively reduced. The proposed algorithm is presented in Alg. 3, and its main computational flow is summarized in the flowchart shown in Fig. 7, reported in Appendix A.3. Note that, while the resampling and the moving step are coherent with the original algorithms described in Sect 3.1 and Sect 3.2, both phases require specifically tailored reweighting steps to work in a setting where the likelihood function is not available in closed form and inference is performed in a sequential data context. 3.3.1 The IBIS step In the standard IBIS algorithm, the reweighting factor turns out to be equal to the predictive likelihood – see Eq (11). This quantity, in the RL framework, has not a closed form. For this reason we propose the following ABC-like estimate for the reweighting factor based on ABC approximations of the two quantities, f(x1:i−1∣μi(ℓ))f(x_1:i-1 μ^( )_i) and f(x1:i∣μi(ℓ))f(x_1:i μ^( )_i) as in Eq (9): f(xi∣x1:i−1,μi(ℓ))≈f~ϵ,M(xi∣x1:i−1,μi(ℓ))=∑m=1Mκ(d(y1:i,m(ℓ),x1:i);ϵ)∑m=1Mκ(d(y1:i−1,m(ℓ),x1:i−1);ϵ).f(x_i x_1:i-1, _i^( ))≈ f_ε,M(x_i x_1:i-1, _i^( ))= Σ _m=1^Mκ (d(y_1:i,m^( ),x_1:i);ε )Σ _m=1^Mκ (d(y_1:i-1,m^( ),x_1:i-1);ε ). A crucial detail lies in the simulation process during the IBIS step – see step 13 in Alg 3. We generate M pseudo-histories for each μi(ℓ) _i^( ) by appending one new agent/environment interaction to the already existing simulated history. Formally, this process is denoted as: y1:i,m(ℓ)∼f(⋅∣μi(ℓ),y1:i−1,m(ℓ))m=1M,\y^( )_1:i,m f(· _i^( ),y^( )_1:i-1,m)\_m=1^M, where each y1:i,m(ℓ)y^( )_1:i,m represents a full trajectory of length i, and where y1:i−1,m(ℓ)y^( )_1:i-1,m denotes the simulated history at the previous step from that specific particle parameter μi(ℓ) _i^( ). 3.3.2 SMC-ABC step Despite appearing as a standard iteration of Algorithm 1, the SMC-ABC steps present some specific features. In particular, the size of the observed data increases throughout the algorithm, which reduces acceptance probabilities and may lead to sample degeneracy [21, 6]. This behavior may have two main drawbacks: (1) the effective sample size can decrease even for a fixed ϵε value; (2) the Metropolis–Hastings acceptance ratio in the moving step may become unstable, even for large values of M. For these reasons, it is necessary to define summary statistics, a distance function, a kernel function, and an ϵε-reduction criterion specifically tailored to this framework. Acceptance-rejection criteria We will use two alternative summary statistics η(⋅)η(·) with different distance functions: • Observation-based criterion: We define an empirical distribution PxP^x, represented as a third-order tensor of dimensions kS×kS×kAk_S× k_S× k_A, where kS:=||k_S:=|S| and kA:=||k_A:=|A|. The entry Pj,k,zxP^x_j,k,z gives the following empirical conditional probability (relative frequency) of transition, estimated from the observed history x1:i=(a1,s1,r1),…,(ai,si,ri)x_1:i=\(a_1,s_1,r_1),…,(a_i,s_i,r_i)\: Pj,k,zx=∑t=1i−1at=z,st=j,st+1=ki−1,P^x_j,k,z\;=\; _t=1^i-1I \a_t=z,s_t=j,s_t+1=k \i-1, where ⋅I\·\ denotes the indicator function. Specifically, the entry in row j and column k corresponds to the empirical frequency of transitions from state j to state k when action z was taken. If a state–action pair is not observed in the history, the corresponding row is conventionally set to zero. Similarly, for each simulated history y1:i,m(ℓ)y^( )_1:i,m, we define a corresponding empirical distribution PyP^y whose entries Pj,k,zyP^y_j,k,z are computed in the same way as Pj,k,zxP^x_j,k,z. We then evaluate the Hellinger distance between these two empirical distributions as d(y1:i,m(ℓ),x1:i)=12∑z∈∑j∈∑k∈(Pj,k,zy−Pj,k,zx)2.d (y^( )_1:i,m,x_1:i )= 12 _z _j _k ( P^y_j,k,z- P_j,k,z^x )^2. (12) • Utility-based criterion: Let U denote the utility function, defined as: U=∑t=1∞γt−1rtU= _t=1^∞γ^t-1r_t, where γ∈[0,1]γ∈[0,1] is the discount factor commonly used in RL, and rtr_t the reward received at time t. Given the observed history x1:ix_1:i and a simulated history y at step i, we compute the corresponding utilities: Ux1:i=∑t=1iγt−1rtx,Uy1:i=∑t=1iγt−1rtyU_x_1:i= _t=1^iγ^t-1r_t^x, U_y_1:i= _t=1^iγ^t-1r^y_t where with rtxr^x_t and rtyr^y_t we indicate the reward at time t in the observed and simulated histories respectively. We then define the discrepancy measure as the Euclidean distance between the two utilities: d(y1:i,m(ℓ),x1:i)=|Ux1:i−Uy1:i,m(ℓ)|d (y^( )_1:i,m,x_1:i )= |U_x_1:i-U_y^( )_1:i,m | (13) This distance captures the difference in cumulative discounted rewards between the observed and simulated histories, emphasizing similarity in overall performance rather than specific sequence details. Kernel function We consider the following hybrid ABC kernel: κ(d(y1:i,m(ℓ),x1:i);ϵi)=1if d(y1:i,m(ℓ),x1:i)≤ϵie−d(y1:i,m(ℓ),x1:i)ϵi2if d(y1:i,m(ℓ),x1:i)>ϵiκ (d(y_1:i,m^( ),x_1:i); _i )= cases1&if d(y_1:i,m^( ),x_1:i)≤ _i\\ e^- d(y_1:i,m^( ),x_1:i) _i^2&if d(y_1:i,m^( ),x_1:i)> _i cases (14) This kernel combines a uniform behavior for distances below the threshold ϵi _i with an exponential decay for distances above it. It enables sharp acceptance for accurate simulations while gradually penalizing less accurate ones. By avoiding a hard cutoff, it assigns non-zero weights even to simulations that slightly exceed the threshold. Since M pseudo-observations are simulated per particle, this kernel allows each of them to contribute to the approximation of the likelihood function in the particle’s weight, even if their distances exceed ϵi _i. Although this contribution is exponentially small, it helps stabilize the acceptance ratio computed in the moving step of Alg 3. Moreover, based on arguments similar to those introduced in [37], we conjecture that it mitigates particle degeneracy issues –typically affecting algorithms based on importance sampling– which, in this setting, are further exacerbated by the fact that the size of the observed data increases throughout the algorithm, thereby raising the probability of generating pseudo-data lying at a distance larger than ϵε from the observed history. However, the use of a kernel defined on a non-compact support introduces a bias; nonetheless, this bias tends to zero as ϵ→0ε→ 0 see Proposition 1. Adapting the threshold To adapt the threshold ϵε, as an alternative to the ESS criterion described in Sect. 3.1, we consider a strategy based on the number of unique particles surviving after the resampling step performed with ϵ∗ε^* [4]. Specifically, let μi(ℓ)ℓ=1L\ _i^( )\_ =1^L be L particles resampled at the i-th iteration according to weights ωi(ℓ)(ϵ∗) _i^( )(ε^*). We define: UP(ωi(ϵ∗)):=μ∈ℳ:∃μi(ℓ)=μ.UP( _i(ε^*)):=\μ :∃ _i^( )=μ\. The number of unique particles corresponds to |UP(ωi(ϵ∗))| |UP( _i(ε^*)) |. The next threshold ϵ∗ε^* is then computed by solving the following equation via the bisection method: |UP(ωi(ϵ∗))|=α⋅|UP(ωi(ϵi))|. |UP ( _i(ε^*) ) |=α· |UP ( _i( _i) ) |. (15) Unlike the ESS, whose evaluation depends on the variability among the particle weights [16], the number of unique particles provides a measure of the structural diversity of the population after the resampling step [4]. Algorithm 3 LF-IBIS. The corresponding flowchart is reported in Appendix A.3, Fig. 7. 1:Input: target threshold ϵε, starting threshold ϵ1 _1, α∈(0,1)α∈(0,1), number of simulations L, number of pseudo-data M, final length of the history n, initial observed data x1x_1; proposal distribution q(⋅,⋅)q(·,·); ϵε adaptation criterion (ESS or UP) 2:Draw μ1(ℓ)∼p(⋅) _1^( ) p(·) for ℓ∈1,…,L ∈\1,...,L\ 3:For each ℓ , simulate y1,m(ℓ)∼f(⋅∣μ1(ℓ))y_1,m^( ) f(· _1^( )) for m∈1,…,Mm∈\1,...,M\ 4:Compute initial weights ω1(ϵ1)=ω1(ℓ)(ϵ1)ℓ=1L _1( _1)=\ _1^( )( _1)\_ =1^L where ω1(ℓ)(ϵ1)=∑m=1Mκ(d(y1,m(ℓ),x1);ϵ1)M _1^( )( _1)= _m=1^Mκ (d(y_1,m^( ),x_1); _1 )M 5:i=1i=1 6:while i≤ni≤ n or ϵi>ϵ _i>ε do 7: if no new agent-env interactions are introduced then perfom a SMC-ABC step then 8: Reweighting: Compute new ϵ∗ε^* by solving with bisection method: ESS(ωi(ϵ∗))=αESS(ωi(ϵi))ESS ( _i(ε^*) )= ( _i( _i) ) or |UP(ωi(ϵ∗))|=α|UP(ωi(ϵi))| |UP ( _i(ε^*) ) |=α |UP ( _i( _i) ) | 9: Where ωi(ℓ)(ϵ∗)=ωi(ℓ)(ϵi)∑m=1Mκ(d(y1:i,m(ℓ),x1:i);ϵ∗)∑m=1Mκ(d(y1:i,m(ℓ),x1:i);ϵi) _i^( )(ε^*)=ω^( )_i( _i) _m=1^Mκ (d(y_1:i,m^( ),x_1:i);ε^* ) _m=1^Mκ (d(y_1:i,m^( ),x_1:i); _i ) 10: Set ϵi=ϵ∗ _i=ε^* 11: else if new data is introduced perfom a IBIS step 12: Set i=i+1i=i+1 and the new ϵi _i equal to the previous one 13: For each ℓ , simulate a new state and append it to y1:i−1,m(ℓ)y^( )_1:i-1,m, for m∈1,…,Mm∈\1,...,M\ : y1:i,m(ℓ)∼f(⋅∣μi−1(ℓ),y1:i−1,m(ℓ))y^( )_1:i,m f(· _i-1^( ),y^( )_1:i-1,m) 14: Set μi(ℓ)ℓ=1L=μi−1(ℓ)ℓ=1L\ _i^( )\_ =1^L=\ _i-1^( )\_ =1^L 15: Compute new weights: ωi(ℓ)(ϵi)=ωi−1(ℓ)(ϵi)∑m=1Mκ(d(y1:i,m(ℓ),x1:i);ϵi)∑m=1Mκ(d(y1:i−1,m(ℓ),x1:i−1);ϵi) _i^( )( _i)=ω^( )_i-1( _i) _m=1^Mκ (d(y_1:i,m^( ),x_1:i); _i ) _m=1^Mκ (d(y_1:i-1,m^( ),x_1:i-1); _i ) 16: end if 17: Resampling: Sample L pairs from (μi(ℓ),y1:i,m(ℓ))ℓ=1L\( _i^( ),y_1:i,m^( ))\_ =1^L using ωi(ℓ)(ϵi)∑ℓ=1Lωi(ℓ)(ϵi) _i^( )( _i) _ =1^L _i^( )( _i) as probabilities. Denote new sampled pairs, with a slight abuse of notation, by (μi(ℓ),y1:i,m(ℓ))ℓ=1L\( _i^( ),y_1:i,m^( ))\_ =1^L. 18: Set new weights: ωi(ℓ)(ϵi)=1ω^( )_i( _i)=1 19: Moving: Draw μ∗∼q(μi(ℓ),⋅)μ^* q( _i^( ),·) for ℓ∈1,…,L ∈\1,...,L\ 20: For each ℓ , simulate y1:i,m∗∼f(⋅∣μ∗)y_1:i,m^* f(· μ^*) for m∈1,…,Mm∈\1,...,M\ 21: Compute the acceptance ratio: ar(μi(ℓ)→μ∗)=p(μ∗)∑m=1Mκ(d(y1:i,m∗,x1:i);ϵi)q(μ∗,μi(ℓ))p(μi(ℓ))∑m=1Mκ(d(y1:i,m(ℓ),x1:i);ϵi)q(μi(ℓ),μ∗)ar( _i^( )→μ^*)= p(μ^*) _m=1^Mκ (d(y_1:i,m^*,x_1:i); _i )q(μ^*, _i^( ))p( _i^( )) _m=1^Mκ (d(y_1:i,m^( ),x_1:i); _i )q( _i^( ),μ^*) for ℓ∈1,…,L ∈\1,...,L\ 22: Draw u(ℓ)∼Unif(0,1)u^( ) (0,1) for ℓ∈1,…,L ∈\1,...,L\ 23: if u(ℓ)<min1,ar(μi(ℓ)→μ∗)u^( )< \1,ar( _i^( )→μ^*)\ then 24: Set μi(ℓ)=μ∗ _i^( )=μ^* and y1:i,m(ℓ)=y1:i,m∗y^( )_1:i,m=y^*_1:i,m for ℓ∈1,…,L ∈\1,...,L\ 25: else 26: Set μi(ℓ)=μi−1(ℓ) _i^( )= _i-1^( ) and y1:i,m(ℓ)=y1:i−1,m(ℓ)y^( )_1:i,m=y_1:i-1,m^( ) for ℓ∈1,…,L ∈\1,...,L\ 27: end if 28: For each μi(ℓ) _i^( ), compute the correspondent optimal policy πμiℓ∗π^*_ _i by solving Eq. (5) 29:end while 3.3.3 Output and stopping criterion LF-IBIS terminates when the maximum number of iterations is reached (i.e., i=ni=n), or when ϵi<ϵ _i<ε. At each iteration i, the algorithm produces a set of L particles μi(ℓ)ℓ=1L\ _i^( )\_ =1^L, which approximate the posterior distribution over the model parameters. For each particle μi(ℓ) _i^( ), the corresponding optimal policy πμi(ℓ)∗π^*_ _i^( ) is computed by solving Eq. (5) via policy iteration, using the transition probabilities Pr(s′∣s,a,μi(ℓ)) (s s,a, _i^( )). As a result, at each iteration we obtain a sample πμi(ℓ)∗ℓ=1L\π^*_ _i^( )\_ =1^L from an approximation of the posterior distribution over optimal policies. This sequential availability of posterior samples over policies enables the application of the decision criterion introduced in Subsection 2.2.1. In particular, these samples are used to evaluate the Bayesian Mean Squared Error (BMSE) and determine whether the current policy should be updated or retained, thereby supporting an adaptive exploration–exploitation strategy. The reliability of these estimates depends on the convergence of the particle system generated by Algorithm 3 to the true posterior distribution. As discussed in the next section, the resulting estimators are expected to be biased but consistent. 3.3.4 Asymptotic properties Algorithm 3 builds upon the IBIS algorithm proposed by Chopin [6] and the SMC-ABC algorithm by Del Moral et al. [10]. As regards SMC-ABC, the algorithm is in turn grounded in the likelihood-based SMC algorithm in Del Moral et al.[9]. While in [10] they only discuss convergence results in their adaptive framework – i.e. when the tolerance threshold is adaptively determined– formal proofs can be found in Beskos et al. [5]. Regarding IBIS, since it is effectively a particle filter algorithm, estimators’ consistency and their asymptotic normality follow from theory about IS and particle filter algorithms –see [7, 14, 6]. Morover, in [6] Theorem 1 shows that, under regularity conditions (see Appendix to [6] ), relative precision of batch importance sampling remains asymptotically stable as the number of observations increases, provided the proportion of new data is constant, which ensures that the number of particles required to maintain accuracy does not need to grow with the total number of observations. However, in our framework, consistency does not directly follow from these results since the reweighting factor is replaced by ABC-like estimates. Given an observed history x∈ix ^i, for same i≥1i≥ 1, an i.i.d. sample of simulated histories y1,…,yMy_1,…,y_M from f(⋅∣μ)f(· μ), and a predefined ABC acceptance region Aϵ,xA_ε,x, we can compute the following two alternative approximate likelihoods: • fϵ,M(x∣μ)f_ε,M(x μ) : the standard ABC estimator defined as in Eq. (8) • f~ϵ,M(x∣μ):=1M∑m=1Mκ(d(ym,x);ϵ) f_ε,M(x μ):= 1M _m=1^Mκ\! (d(y_m,x);\,ε ), our ABC estimator based on the kernel κ(⋅;⋅)κ(·;·) introduced in Eq. (14). In what follows Proposition 1 states that our approximate likelihood f~ϵ,M(x∣μ) f_ε,M(x μ) is a biased estimate for fϵ,M(x∣μ)f_ε,M(x μ) and provides an upper bound for that bias. Proposition 2 guarantees the convergence of the corresponding approximate IBIS reweighting factor to the exact one. Moreover, empirical evidence of convergence is provided through the experiments in Sect 4.2. Proposition 1. Let β(μ;ϵ,M):=|f~ϵ,M(x∣μ)−fϵ,M(x∣μ)|β(μ;ε,M):= | f_ε,M(x μ)-f_ε,M(x μ) | be the bias induced by the proposed kernel function. For each μ∈ℳμ , the bias is bounded from above by the following quantity that depends on μ and ϵε: β(μ;ϵ,M)≤L(μ;ϵ,M)=1M∑m=1Mexp[−d(x,ym)ϵ2].β(μ;ε,M)\;≤\;L(μ;ε,M)= 1M _m=1^M \! [- d(x,y_m)ε^2 ]. Proof of Proposition 1 is provided in the Supplementary Material (Appendix A.2). Note that this upper bound corresponds to the extreme case in which all histories simulated under μ are far from the observed history and would have been rejected under the standard uniform kernel. From Proposition 1, it also follows that as ϵε goes to zero, β(μ;ϵ,M)⟶0β(μ;ε,M) 0. Proposition 2. Let x1:ix_1:i and x1:i−1x_1:i-1 denote the observed histories up to iteration i and i−1i-1, respectively. For a given parameter vector μ, let y1:i−1,mm=1M\y_1:i-1,m\_m=1^M be i.i.d. simulated histories generated from f(⋅∣μ)f(· μ) and let each yi,my_i,m be independently sampled from f(⋅∣y1:i−1,m,μ)f(· y_1:i-1,m,μ). Define the approximate predictive likelihood used in the IBIS reweighting step as f~ϵ,M(xi∣x1:i−1,μ)=f~ϵ,M(x1:i∣μ)f~ϵ,M(x1:i−1∣μ)=∑m=1Mκ(d(y1:i,m,x1:i);ϵ)∑m=1Mκ(d(y1:i−1,m,x1:i−1);ϵ), f_ε,M(x_i x_1:i-1,μ)= f_ε,M(x_1:i μ) f_ε,M(x_1:i-1 μ)= _m=1^Mκ(d(y_1:i,m,x_1:i);ε) _m=1^Mκ(d(y_1:i-1,m,x_1:i-1);ε), where κ(⋅,⋅)κ(·,·) is defined in Eq. (14) and d(⋅,⋅)d(·,·) denotes the Euclidean distance. Then, limM→∞limϵ→0f~ϵ,M(xi∣x1:i−1,μ)=f(xi∣x1:i−1,μ)a.s. _M→∞ _ε→ 0 f_ε,M(x_i x_1:i-1,μ)=f(x_i x_1:i-1,μ) .s. provided that f(x1:i−1∣μ)>0f(x_1:i-1 μ)>0. Proof. Consider the numerator f~ϵ,M(x1:i∣μ) f_ε,M(x_1:i μ). We first take the limit ϵ→0ε→ 0. From Proposition 1, it follows that limϵ→0f~ϵ,M(x1:i∣μ) _ε→ 0 f_ε,M(x_1:i μ) =limϵ→0[fϵ,M(x1:i∣μ)+β(μ;ϵ,M)] = _ε→ 0 [f_ε,M(x_1:i μ)+β(μ;ε,M) ] =limϵ→0fϵ,M(x1:i∣μ) = _ε→ 0f_ε,M(x_1:i μ) =limϵ→01M∑m=1MAϵ,x1:i(y1:i,m) = _ε→ 0 1MΣ _m=1^MI_A_ε,x_1:i(y_1:i,m) =1M∑m=1Mx1:i(y1:i,m) = 1MΣ _m=1^MI_x_1:i(y_1:i,m) where x1:i(y1:i,m)I_x_1:i(y_1:i,m) denotes the indicator function that takes the value 1 when y1:i,m=x1:iy_1:i,m=x_1:i and 0 otherwise. This function is a Bernoulli random variable whose expected value is the probability f(x1:i∣μ)f(x_1:i μ). Hence, from the Strong Law of Large Numbers it follows that limM→∞limϵ→0f~ϵ,M(x1:i∣μ)=limM→∞1M∑m=1Mx1:i(y1:i,m)=f(x1:i∣μ)a.s. _M→∞ _ε→ 0 f_ε,M(x_1:i μ)= _M→∞ 1MΣ _m=1^MI_x_1:i(y_1:i,m)=f(x_1:i μ) a.s. The same arguments apply to the denominator: limM→∞limϵ→0f~ϵ,M(x1:i−1∣μ)=f(x1:i−1∣μ)a.s _M→∞ _ε→ 0 f_ε,M(x_1:i-1 μ)=f(x_1:i-1 μ) a.s Since in Algorithm 3, the condition f(x1:i−1∣μ)>0f(x_1:i-1 μ)>0 is always satisfied because the observed history has strictly positive probability under the model, we have that limM→∞limϵ→0f~ϵ,M(x1:i∣μ)f~ϵ,M(x1:i−1∣μ)=f(x1:i∣μ)f(x1:i−1∣μ)=f(xi∣x1:i−1,μ)a.s. _M→∞ _ε→ 0 f_ε,M(x_1:i μ) f_ε,M(x_1:i-1 μ)= f(x_1:i μ)f(x_1:i-1 μ)=f(x_i x_1:i-1,μ) a.s. ∎ Corollary 3. Under the assumptions of Proposition 2, let ω~i(ℓ)=ω~i−1(ℓ)f~ϵ,M(xi∣x1:i−1,μ(ℓ))andωi(ℓ)=ωi−1(ℓ)f(xi∣x1:i−1,μ(ℓ)) ω_i^( )= ω_i-1^( )\, f_ε,M(x_i x_1:i-1,μ^( ))\;\;and\;\; _i^( )= _i-1^( )\,f(x_i x_1:i-1,μ^( )) denote the approximate importance weight computed in Algorithm 3, and corresponding exact IBIS importance weight, respectively. Then, for every particle μ(ℓ)μ^( ), limM→∞limϵ→0ω~i(ℓ)=ωi(ℓ)a.s. _M→∞ _ε→ 0 ω_i^( )= _i^( ) .s. Proof. The result follows immediately from Proposition 2 by substituting the approximate predictive likelihood with its limiting value in the recursive weight update. ∎ Corollary 3 provides a theoretical justification for the reweighting step employed by Algorithm 3. In the asymptotic regime M→∞M→∞ and ϵ→0ε→ 0, the proposed likelihood-free reweighting scheme recovers the exact IBIS reweighting mechanism. 4 Experiments We evaluate our algorithm in a controlled setting –where exact Bayesian inference is tractable– and in scenarios where no closed-form posterior is available. This allows us to assess three complementary aspects. First, by comparing LF-IBIS with the exact posterior, we verify that the approximation induced by the ABC tolerance ϵε remains negligible and that the method provides a good approximation to the true posterior distribution. Second, we compare LF-IBIS against acceptance–rejection ABC (AR-ABC), which represents the state of the art in likelihood-free inference in RL, as used in [13], and is executed in an offline manner using the complete observed history available at the final time. Third, for both summary statistics and distance metric considered—the Hellinger distance on observation-based statistics and the Euclidean distance on utility-based statistics—we report results for LF-IBIS under the two adaptive thresholding criteria: the one based on the ESS and the one based on the number of UP. The experiments are conducted in an online setting in which data arrive sequentially through real-time interactions with the environment, while keeping the policy fixed throughout the entire execution. In a separate experiment, we show how to exploit the full posterior distribution of the optimal policy for adaptation, demonstrating that LF-IBIS naturally supports the exploration–exploitation trade-off as posterior uncertainty is progressively reduced. Once no further data are available, the algorithm continues to decrease the tolerance ϵε until the target value is reached. For comparability across methods, the observed history is kept identical across all experiments. Tuning parameters The parameter M is fixed to M=50M=50 across all experiments, while the number of particles L is chosen according to the complexity of the problem and, in particular, to the informativeness of the summary statistics: larger values of L are used when the statistics are less informative, in order to ensure a more accurate approximation of the posterior distribution. Finally, the parameter α is calibrated empirically in a preliminary phase, once M and L are fixed. In particular, α is selected so as to control the proportion of particles discarded at each iteration, ensuring that it is appropriate for the complexity of the experiment and for the chosen values of M and L. Table 1 reports the tuning parameters used in each experimental setting. Table 1: Tuning parameters used in Experiments 1 and 2. Experiment Method Parameter L M α Exp. 1 (Observation-based) LF-IBIS (ESS) μ 50,000 50 0.99 LF-IBIS (UP) μ 50,000 50 0.98 LF-IBIS (ESS) β 50,000 50 0.99 LF-IBIS (UP) β 50,000 50 0.98 Exp. 2 (Utility-based) LF-IBIS (ESS) μ 100,000 50 0.92 LF-IBIS (UP) μ 100,000 50 0.95 LF-IBIS (ESS) β 100,000 50 0.93 LF-IBIS (UP) β 100,000 50 0.95 4.1 Experimental Setting Response-Adaptive Randomization (RAR) is a special adaptive design used in clinical trials to update the allocation to treatments of incoming study participants to favour the most promising and beneficial treatments as indicated by interim results [38, 27]. The application of RL to RAR is promising for the potentiality of RL to allow dynamic updates to randomization probabilities based on real‐time predictions of treatment success [11, 22]. In this setting, the environment represents patients, with states corresponding to their treatment responses. We assume there are N patients entering the trial sequentially, and that each patient’s response is observed before the allocation of the next patient. The agent’s goal is to assign the n-th patient (1≤n≤N1≤ n≤ N) to treatment or placebo in order to maximize the number of disease remissions. Here, the action space is =0,1A=\0,1\, and each action AnA_n is a Bernoulli random variable representing the assignment of the n-th patient to treatment (An=1A_n=1) or control (An=0A_n=0). The policy π∈[0,1]π∈[0,1] corresponds to the Bernoulli parameter governing the i.i.d. sequence of random variables Ann=1N\A_n\_n=1^N. The environment is characterized by two possible states: =0,1S=\0,1\. If the n-th patient achieves remission, the state is Sn=1S_n=1; otherwise, Sn=0S_n=0. The responses of the environment are two Bernoulli random variables: Z0∼Ber(μ0)Z_0 ( _0) and Z1∼Ber(μ1)Z_1 ( _1), where Z0Z_0 denotes Sn∣An=0S_n A_n=0 and Z1Z_1 denotes Sn∣An=1S_n A_n=1, while μ0 _0 and μ1 _1 are the probabilities of remission in the control and treatment group, respectively. The reward depends on both the observed state and the treatment decision. Let us suppose that in our experiments a collateral effect may occur with probability 0.70.7 when treatment is administered. In the absence of treatment, it is assumed that the reward coincides with the state, namely Rn(an=0,sn=s)=sfor s∈0,1.R_n(a_n=0,s_n=s)=s s∈\0,1\. When treatment is administered, conditional on treatment, the reward is assumed to be: Rn(an=1,sn=0)=−0.2,if a collateral effect occurs,0,if no collateral effect occurs.R_n(a_n=1,s_n=0)= cases-0.2,&if a collateral effect occurs,\\ 0,&if no collateral effect occurs. cases Rn(an=1,sn=1)=0.8,if a collateral effect occurs,1,if no collateral effect occurs.R_n(a_n=1,s_n=1)= cases0.8,&if a collateral effect occurs,\\ 1,&if no collateral effect occurs. cases Note that, while in the absence of treatment the reward is deterministic, under treatment the reward is a random variable that depends on both the observed state sns_n and the treatment decision ana_n, as well as on the independently sampled collateral effect. The observed history at time n (corresponding to the n-th patient) is defined as x1:n=(a1,s1,r1),…,(an,sn,rn)∈nx_1:n=\(a_1,s_1,r_1),…,(a_n,s_n,r_n)\ ^n, where rnr_n denotes the realization for Rn(an,sn)R_n(a_n,s_n). Our objective is to perform inference on μ0 _0 and μ1 _1 in order to learn the optimal policy, conditional on the history up to time n: p(π∗,μ0,μ1∣x1:n)p(π^*, _0, _1 x_1:n). We assume uninformative Beta priors for the parameters, i.e., μ0,μ1,π∼Beta(1,1) _0, _1,π (1,1). To move the particles, we construct proposal distributions of the form q(μ0,i−1(ℓ),⋅)q( _0,i-1^( ),·), which at iteration i are defined as Beta distributions whose parameters are chosen via the method of moments. This ensures that the proposal is centered at the previously accepted particle μ0,i−1(ℓ) _0,i-1^( ) and that its variance matches twice the empirical variance of the set of previously accepted particles μ0,i−1(ℓ)ℓ=1L\ _0,i-1^( )\_ =1^L. An analogous construction is used for q(μ1,i−1(ℓ),⋅)q( _1,i-1^( ),·). Although such Beta proposal distributions are not standard, their design is inspired by the adaptive proposal strategy widely used for multivariate normal proposals, as popularized in the adaptive ABC framework of Beaumont et al. [2]. Full details of the proposal definition are provided in the Appendix B.1. The clinical application described so far can also be faced by expressing μ0 _0 and μ1 _1 as functions of alternative parameters: μ0=11+e−β0,μ1=11+e−(β0+β1), _0= 11+e^- _0, _1= 11+e^-( _0+ _1), with β0,β1∈ℝ _0, _1 . This inverse-logit (expit) transformation allows us to work with parameters defined on the entire real line. Under this parametrization, one can adopt normal priors for β0 _0 and β1 _1; however, these priors are not conjugate, and as a result the posterior distributions do not have closed-form expressions. An additional advantage of this logistic representation is that it naturally accommodates the inclusion of covariates (e.g., patient-specific characteristics) [12]. Under this parametrization, the target distribution becomes the posterior p(π,β0,β1∣x1:n)p(π, _0, _1 x_1:n) from which p(π,μ0,μ1∣x1:n)p(π, _0, _1 x_1:n) can be recovered through the expit transformation. The prior distributions used in this setting are β0∼(0,16),β1∼(2,16). _0 (0,16), _1 (2,16). The proposal distribution at iteration i is defined as a bivariate Normal distribution centered at the empirical mean of the particles alive at the previous iteration m^i−1 m_i-1, with covariance matrix Σ^i−1 _i-1 equal to twice their empirical covariance [2]: q((β0,i(ℓ),β1,i(ℓ))∣(β0,i−1(ℓ),β1,i−1(ℓ))ℓ=1L)=(m^i−1, 2Σ^i−1).q (( _0,i^( ), _1,i^( )) \( _0,i-1^( ), _1,i-1^( ))\_ =1^L )=N ( m_i-1,\;2 _i-1 ). 4.2 Experimental Results In the following, we present results for the described application in both its original formulation and the version based on logistic transformation. The observed history x1:Nx_1:N was generated using the true values μ0true=0.3 _0^true=0.3 and μ1true=0.8 _1^true=0.8, that, under the logistic formulation, correspond to the following values for β0true _0^true and β1true _1^true: β0true=logit(μ0true)≈−0.85,β1true=logit(μ1true)−β0true≈2.24. _0^true=logit( _0^true)≈-0.85, _1^true=logit( _1^true)- _0^true≈ 2.24. In the first experiment with observation-based statistics and Hellinger distance, the final threshold ϵfinal _final was determined from the AR-ABC simulations. Specifically, ϵfinal _final was chosen as the empirical 0.15%0.15\% quantile of the simulated distance distribution. In more challenging settings, such as the utility-based experiment, AR-ABC becomes substantially less efficient. In such cases, ϵfinal _final is calibrated empirically by considering candidate tolerance values and assessing the quality of the resulting posterior approximations. Whenever the exact Bayesian posterior is available, it is used as a benchmark to evaluate the agreement between the approximate and exact posteriors. 4.2.1 Experiment 1: Hellinger distance on observation-based statistics Figure 1 presents the posterior distributions for μ0 _0 and μ1 _1 obtained using LF-IBIS with ESS and LF-IBIS with UP (left panels), along with the corresponding learned policies (right panels). Both variants of the LF-IBIS algorithm yield posterior approximations that are consistent with exact Bayesian inference, similarly AR-ABC approach. The final threshold value was ϵfinal=0.02165 _final=0.02165. With this threshold, AR-ABC was run with 3M simulated samples, yielding 4,936 accepted particles. For the LF-IBIS experiments (both ESS- and UP-based), the initial history length was set to 3, the initial threshold to ϵ=1ε=1, and the final history length to 48. The tuning parameters M, L, and α were set as reported in Table 1. The average execution time was approximately 23 minutes for LF-IBIS with ESS and 29 minutes for LF-IBIS with UP. The total number of iterations required to reach the target tolerance was 94 for the ESS-based version and 89 for the UP-based version. Figure 1: Posterior distributions of μ (left panel) and the corresponding policies (right panel), when using Hellinger distance on observation-based statistics. Light blue: LF-IBIS with ESS. Dark blue: LF-IBIS with UP. Orange: exact Bayesian inference. Green: AR-ABC. The dashed vertical lines indicate μ0true _0^true and μ1true _1^true in the left panel, and the corresponding true posterior optimal policy π∗π^* in the right panel, computed via policy iteration using μ0true _0^true and μ1true _1^true. It is worth noting that the posterior distributions obtained with LF-IBIS appear to be slightly more concentrated than those produced by both exact Bayesian inference and the AR-ABC method. This difference is expected, as both exact Bayesian inference and AR-ABC are performed in an offline setting using the complete observed history, whereas LF-IBIS processes observations sequentially as they become available. Consequently, the posterior distributions obtained by these approaches cannot be expected to coincide exactly. Nevertheless LF-IBIS yields posterior distributions that are very close to those of the offline methods and remain highly consistent in shape. Figure 2 reports the results obtained when the problem is reparameterized using the logistic transformation. In this setting, the reparameterization makes exact Bayesian inference intractable. The parameter settings for both LF-IBIS and AR-ABC are identical to those used in Figure 1, and the same considerations regarding the quality of the approximations apply here as well. With the same threshold value ϵ=0.02165ε=0.02165, the AR-ABC procedure accepted 1,355 samples out of 3M simulated draws. The average execution time was approximately 35 minutes with 125 iterations for LF-IBIS with ESS, and approximately 36 minutes with 108 iterations for LF-IBIS with UP. Another noteworthy point is that the UP strategy seems to require fewer iterations to reach the target tolerance ϵε, in both the current example reparameterized via the logistic transformation and in the previous one. Figure 2: Posterior distributions of β0 _0 and β1 _1 (left panel) and the corresponding policies (right panel), when using Hellinger distance on observation-based statistics. Light blue: LF-IBIS with ESS. Dark blue: LF-IBIS with UP. Green: AR-ABC. The dashed red vertical line indicates β0true _0^true and β1true _1^true in the left panel, and the corresponding true posterior optimal policy π∗π^* in the right panel, computed via policy iteration using β0true _0^true and β1true _1^true. 4.2.2 Experiment 2: Euclidean distance on utility-based statistics Figure 3 shows the results obtained using utility-based statistics. In this case, all results were produced with a final threshold value of ϵ=2.1×10−5ε=2.1× 10^-5. For LF-IBIS with ESS, the initial history length was set to 3, the initial threshold to ϵ=0.5ε=0.5, and the final history length to 48. The tuning parameters M, L, and α were set as reported in Table 1. The average execution time was approximately 16 minutes for LF-IBIS with ESS (with 103 iterations) and approximately 27 minutes for LF-IBIS with UP (with 110 iterations). Figure 3: Posterior distributions of μ (left panel) and the corresponding policies (right panel), when using Euclidean distance on utility-based statistics. Light blue: LF-IBIS with ESS. Dark blue: LF-IBIS with UP. Orange: exact Bayesian inference. Green: AR-ABC. The dashed red vertical line indicates μ0true _0^true and μ1true _1^true in the left panel, and the corresponding true posterior optimal policy π∗π^* in the right panel, computed via policy iteration using μ0true _0^true and μ1true _1^true. In this experiment, unlike in the previous one, the posterior approximation obtained with LF-IBIS appears to be more accurate than that produced by AR-ABC. In the previous experiment, observation-based statistics were employed; these are constructed from empirical distributions and therefore tend to preserve more information about the underlying parameters μ. In contrast, the utility-based statistics used here constitute a suboptimal summary of the data, as they aggregate information into coarse measures and consequently convey less information about μ. In this setting, the sequential (online) nature of LF-IBIS, as opposed to the fully offline approach adopted by AR-ABC—which processes the entire observed history at once—appears to effectively compensate for the reduced informativeness of the utility-based statistics, resulting in a more coherent posterior approximation. Despite being less informative than observation-based summaries, utility-based statistics can be considerably more practical in applied settings. Our experimental framework is intentionally simplified, as we consider an MDP for which informative statistics of the process are known. In more realistic applications, however, such knowledge may be unavailable. For instance, in model-free contexts or in problems involving complex or partially observed dynamics, identifying sufficient statistics can be extremely challenging. In such cases, utility-based statistics offer a flexible alternative. An additional consequence of the reduced informativeness of the utility-based statistics is the need for a larger particle population to achieve stable inference. Specifically, we employ L=100,000L=100,000 particles, compared to L=50,000L=50,000 in the observation-based setting. Nevertheless, utility-based statistics are computationally less demanding, as they do not require the construction of empirical distributions at each iteration. This results in shorter execution times: LF-IBIS with utility-based statistics completes in approximately 16–27 minutes (for both the ESS- and UP-based variants), compared to roughly 23–29 minutes for the observation-based setting. Finally, the tuning parameter α was reduced relative to the previous experiments. In the current study, we set α=0.92α=0.92 for the ESS-based variant and α=0.95α=0.95 for the UP-based variant. Smaller values of α allow particles with negligible importance weights to be pruned and replaced earlier in the procedure, thereby accelerating adaptation of the particle population to the observed data. Figure 4 reports the results obtained with utility-based statistics when the problem is reparameterized through the logistic transformation. The parameter values for LF-IBIS and AR-ABC are the same as in Figure 3, and the same considerations regarding the results apply here as well. For the AR-ABC, 693 samples are accepted. The execution time is ≈18≈ 18 minutes with 109 iterations for the ESS version and ≈29≈ 29 minutes with 116 iterations for the UP version. Figure 4: Posterior distributions of β0 _0 and β1 _1 (left panel) and the corresponding policies (right panel), when using Euclidean distance on utility-based statistics. Light blue: LF-IBIS with ESS. Dark blue: LF-IBIS with UP. Green: AR-ABC. The dashed red vertical line indicates β0true _0^true and β1true _1^true in the left panel, and the corresponding true posterior optimal policy π∗π^* in the right panel, computed via policy iteration using β0true _0^true and β1true _1^true. Table 2 summarizes the experimental settings adopted in the two experiments. The table includes both algorithmic hyperparameters (e.g., particle population size and tuning parameters) and execution-related measures, such as the number of LF-IBIS iterations and runtime, as well as the final ESS for LF-IBIS and the total number of simulated samples and accepted particles for AR-ABC. It is important to note that in both the observation-based and utility-based experiments, although the UP variant often converges in the same or fewer iterations compared to the ESS-based approach, it is generally more time-consuming overall. This is because each application of the bisection method requires an additional resampling step. Table 2: Summary of experimental settings for Experiments 1 and 2. Experiment Method Param. L M α ϵfinal _final Iterations / Time Final ESS AR-ABC Simulated / Accepted Exp. 1 (Obs.-based) LF-IBIS (ESS) μ 50,000 50 0.99 0.02165 94 / ∼ 23 min 19576 LF-IBIS (UP) μ 50,000 50 0.98 0.02165 89 / ∼ 29 min 18701 AR-ABC μ – – – 0.02165 ∼ 13 sec 3M / 4,936 LF-IBIS (ESS) β 50,000 50 0.99 0.02165 125 / ∼ 35 min 5732 LF-IBIS (UP) β 50,000 50 0.98 0.02165 108 / ∼ 36 min 5758 AR-ABC β – – – 0.02165 ∼ 18 sec 3M / 1,355 Exp. 2 (Utility-based) LF-IBIS (ESS) μ 100,000 50 0.92 2.1×10−52.1× 10^-5 103 / ∼ 16 min 260 LF-IBIS (UP) μ 100,000 50 0.95 2.1×10−52.1× 10^-5 110 / ∼ 27 min 449 AR-ABC μ – – – 2.1×10−52.1× 10^-5 ∼ 4 min 150M / 725 LF-IBIS (ESS) β 100,000 50 0.93 2.1×10−52.1× 10^-5 109 / ∼ 18 min 381 LF-IBIS (UP) β 100,000 50 0.95 2.1×10−52.1× 10^-5 116 / ∼ 29 min 497 AR-ABC β – – – 2.1×10−52.1× 10^-5 ∼ 6 min 150M / 693 As a final assessment, to quantify the discrepancy between the approximate posterior distributions on μ produced by AR-ABC and by the two variants of LF-IBIS (ESS- and UP-based) and the exact Bayesian posterior, we use the Energy Distance [34] for both experiments. The energy distance is a metric that measures the distance between probability distributions based on pairwise distances between samples. For each method, we report the mean and standard deviation of the Energy distance computed over ten independent runs of the algorithm, all performed using the same observed history as in the experiments above. The Energy distance is reported only for the non-reparameterized versions of the experiments, since under the logistic reparameterization the exact Bayesian posterior is not available. Overall, the UP-based variant tends to produce posterior distributions that are closer to the exact Bayesian posterior. Moreover, the approximations obtained with UP exhibit lower variability across independent runs, and this difference is particularly pronounced in Experiment 2. In Appendix B.2, we also report the KL divergences between the posteriors of μ and π. Table 3: Mean (standard deviation) of the Energy distance computed over ten independent runs of the algorithms with respect to the exact Bayesian posterior. Results are reported for inference on the parameter μ and the policy π. Method Exp. 1 (Obs.-based) Exp. 2 (Utility-based) μ π μ π LF-IBIS (ESS) 0.0033 (0.0005) 0.0342 (0.0040) 0.0178 (0.0172) 0.0676 (0.0402) LF-IBIS (UP) 0.0026 (0.0003) 0.0265 (0.0029) 0.0031 (0.0025) 0.0273 (0.0131) AR-ABC 0.0001 (0) 0.0040 (0.0011) 0.1335 (0.0076) 0.2199 (0.0064) 4.3 Updating policy from posterior distribution The posterior distribution over policies provides valuable information for deciding when it is advantageous to update the current policy. In this section, policy updates are performed according to the criterion introduced in Section 2.2.1. At each iteration, samples from the posterior distribution of optimal policies are used to evaluate the Bayesian Mean Squared Error (BMSE) for both the current policy and a candidate policy drawn from the posterior. The policy with the lower BMSE is selected, as it is, on average, closer to the set of plausible optimal policies. The Bayesian Mean Squared Error for a policy π is defined in Eq (6) Intuitively, when the posterior distribution assigns high probability mass to the current policy, exploration of the environment is encouraged; otherwise, the policy is updated to exploit the information contained in the posterior. In Figure 5, we illustrate the value function as a function of the policy parameter π for Experiment 1, described in Section 4.2.1, when LF-IBIS is run using the ESS-based criterion for decreasing the tolerance threshold ϵε. The blue curve represents the posterior mean of the value function, while the shaded region corresponds to the 90%90\% posterior credibility bands. The dashed red vertical line indicates the true posterior optimal policy, computed via policy iteration using the true environment parameters μ0true _0^true and μ1true _1^true. The posterior distribution of the value function is concentrated around a mean function that assigns higher values to policies in the neighborhood of the true optimal policy. This highlights a novel aspect of our approach, namely the ability to obtain a posterior distribution over the value function by propagating uncertainty in the parameters governing the environment dynamics. Figure 5: Value function when the policy π varies from 0 to 1. Blue curve: posterior mean of the value function with 90%90\% posterior credibility bands. In Figure 6, we compare the value function obtained when the policy is kept fixed throughout the experiment with that obtained after switching to the updated policy. The results refer to Experiment 1 (Section 4.2.1) ESS version. We also report the corresponding number of patients assigned to the treatment group, with and without policy switching. Credible bands are computed using posterior samples of the environment parameters μ. The figure shows that, from a practical perspective, switching the policy is relevant. Whenever an improvement in the value function is observed, a corresponding increase in the number of patients assigned to the treatment group is also observed. The variability around the value function reflects the posterior uncertainty of the environment parameters μ: when the history length increases, this variability progressively decreases. In this setting, it is possible to plot the value function because the state space consists of only two states, 0,1\0,1\, and—by construction—the value function does not differ between them. When the blue curve dominates the red one, a clear separation between the two value functions is observed, whereas the converse is not true. Moreover, the blue curve almost never falls below the red one, except for a negligible deviation at a single point. Figure 6: Value function (top panel) and number of patients assigned to treatment (bottom panel) with and without policy switch (see Figure 1, ESS version). Here posterior sampling is performed after every 2 new interactions with the environment. Vertical lines indicate decision points: grey lines represent the time steps (i.e., lengths of the observed history) at which a policy switch is evaluated, while blue lines denote the time steps at which a policy switch is actually performed. 5 Discussion and Conclusion In this work, we introduced a fully Bayesian Reinforcement Learning (fBRL) procedure based on Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), designed to address data scarcity and to support online learning in RL. The proposed framework enables posterior inference over both environment parameters and policies without requiring access to an explicit likelihood. We evaluated LF-IBIS in a controlled experimental setting where exact Bayesian inference is available, thereby allowing a direct assessment of the quality of the approximation. Despite the typical sources of approximation introduced by ABC—most notably the choice of the tolerance ϵε, distance functions, and summary statistics [30]—our results show that the approximate posteriors produced by LF-IBIS closely match the exact posterior distributions. Moreover, we demonstrated that processing data incrementally yields posterior distributions that are comparable to those obtained when the entire observed history is processed at once, thus confirming the validity of the online formulation. Beyond the empirical validation, we also established asymptotic properties of the proposed likelihood-free reweighting scheme. Specifically, we showed that the ABC approximation provides a consistent approximation of the predictive likelihood and that the resulting importance weights converge almost surely to those of the exact IBIS algorithm as the ABC tolerance tends to zero and the number of simulated histories increases. A further contribution of this work is the comparison between observation-based and utility-based summary statistics. Observation-based statistics, constructed from empirical outcome distributions, are more informative and lead to highly accurate posterior approximations. However, they rely on knowledge of sufficient statistics, which may not be available in realistic or complex environments. Utility-based statistics, while less informative due to their aggregated nature, provide a practical and flexible alternative that is computationally more efficient and does not require explicit modeling of the underlying process. Our experimental results indicate that the incremental, sequential nature of LF-IBIS partially compensates for the reduced informativeness of utility-based statistics, yielding shorter execution times while maintaining competitive posterior accuracy. In addition, we investigated the ability of LF-IBIS to support adaptive decision-making through posterior-driven policy updates. When policies are sampled from their posterior distributions, the algorithm implements a Bayesian exploration–exploitation strategy: early posterior uncertainty promotes exploratory behavior, while posterior concentration gradually shifts the agent toward exploitation. This demonstrates that LF-IBIS not only performs inference over environment parameters, but also provides a principled mechanism for policy adaptation in online settings. More generally, the proposed framework is not limited to MDP with analytically tractable dynamics. Whenever a simulator of the environment is available—even if the underlying process does not admit an explicit MDP formulation—the LF-IBIS methodology can be naturally extended to this model-free RL setting. Although in our experimental setting the likelihood can be written analytically—due to the assumption of a known probability of side effects in the reward function—the proposed approach does not fundamentally rely on this property. In more general scenarios, LF-IBIS can be applied in exactly the same way without requiring explicit specification of the transition dynamics or the side-effect generation mechanism, provided that a simulator is available. This highlights the fully likelihood-free nature of the method and its suitability for settings in which analytical modeling is infeasible or undesirable. Overall, our results indicate that LF-IBIS constitutes a flexible and robust framework for likelihood-free Bayesian Reinforcement Learning, capable of supporting online inference, policy adaptation, and principled uncertainty quantification in complex decision-making problems. References [1] P. Auer, N. Cesa-Bianchi, and P. Fischer (2002) Finite-time analysis of the multiarmed bandit problem. Machine Learning 47 (2–3), p. 235–256. External Links: Document Cited by: §1. [2] M. A. Beaumont, J.-M. Cornuet, J.-M. Marin, and C. P. Robert (2009) Adaptive approximate bayesian computation. Biometrika 96 (4), p. 983–990. External Links: Document Cited by: §4.1, §4.1. [3] R. Bellman (1957) A markovian decision process. Journal of Mathematics and Mechanics 6 (5), p. 679–684. External Links: Link Cited by: §2.1, §2.1. [4] E. Bernton, P. E. Jacob, M. Gerber, and C. P. Robert (2019) Approximate bayesian computation with the wasserstein distance. Journal of the Royal Statistical Society: Series B 81 (2), p. 235–269. External Links: Document Cited by: §3.3.2, §3.3.2. [5] A. Beskos, A. Jasra, N. Kantas, and A. Thiery (2016) On the convergence of adaptive sequential monte carlo methods. The Annals of Applied Probability 26 (2), p. 1111–1146. External Links: Document Cited by: §3.3.4. [6] N. Chopin (2002) A sequential particle filter method for static models. Biometrika 89 (3), p. 539–551. External Links: Document Cited by: §1, §1, §3.2, §3.3.2, §3.3.4. [7] D. Crisan and A. Doucet (2000) Convergence of sequential monte carlo methods. Technical Report Technical Report CUED/F-INFENG/TR381, University of Cambridge. Cited by: §3.3.4. [8] R. Dearden, N. Friedman, and D. Andre (1999) Model-based bayesian exploration. In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, p. 150–159. External Links: Document Cited by: §1. [9] P. Del Moral, A. Doucet, and A. Jasra (2006) Sequential monte carlo samplers. Journal of the Royal Statistical Society: Series B 68 (3), p. 411–436. External Links: Document Cited by: §3.3.4. [10] P. Del Moral, A. Doucet, and A. Jasra (2012) An adaptive sequential monte carlo method for approximate bayesian computation. Statistics and Computing 22 (5), p. 1009–1020. External Links: Document Cited by: §A.1, §A.1, §1, §3.1, §3.1, §3.1, §3.3.4. [11] N. Deliu and S. S. Villar (2025) On the finite-sample and asymptotic error control of a randomization-probability test for response-adaptive clinical trials. Biometrics 81 (2), p. ujaf069. External Links: Document Cited by: §4.1. [12] N. Deliu (2021) Reinforcement learning in modern biostatistics: benefits, challenges and new proposals. Ph.D. Thesis, Sapienza University of Rome. External Links: Link Cited by: §4.1. [13] C. Dimitrakakis and N. Tziortziotis (2013) ABC reinforcement learning. In Proceedings of the 30th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 28, p. 684–692. External Links: Link Cited by: §1, §1, §3.1, §4. [14] A. Doucet, N. de Freitas, and N. Gordon (Eds.) (2001) Sequential monte carlo methods in practice. Statistics for Engineering and Information Science, Springer, New York, NY. External Links: Document Cited by: §3.3.4. [15] M. O. Duff (2002) Optimal learning: computational procedures for bayes-adaptive markov decision processes. Ph.D. Thesis, University of Massachusetts Amherst. Cited by: §2.2. [16] V. Elvira, L. Martino, and C. P. Robert (2022) Rethinking the effective sample size. International Statistical Review 90 (3), p. 525–550. External Links: Document Cited by: §3.3.2. [17] Z. Ghahramani (2015) Probabilistic machine learning and artificial intelligence. Nature 521 (7553), p. 452–459. External Links: Document Cited by: §1. [18] M. U. Gutmann and J. Corander (2016) Bayesian optimization for likelihood-free inference. Journal of Machine Learning Research 17 (125), p. 1–47. Cited by: §1. [19] T. Haarnoja, H. Tang, P. Abbeel, and S. Levine (2017) Reinforcement learning with deep energy-based policies. In Proceedings of the 34th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 70, p. 1352–1361. External Links: Link Cited by: §2.2. [20] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine (2018) Soft actor-critic: off-policy maximum entropy deep reinforcement learning with a stochastic actor. In Proceedings of the 35th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 80, p. 1861–1870. External Links: Link Cited by: §2.1. [21] J. S. Liu and R. Chen (1998) Sequential monte carlo methods for dynamic systems. Journal of the American Statistical Association 93 (443), p. 1032–1044. External Links: Document Cited by: §3.3.2. [22] D. Merrell, T. Chandereng, and Y. Park (2023) A markov decision process for response-adaptive randomization in clinical trials. Computational Statistics & Data Analysis 178, p. 107599. External Links: Document Cited by: §4.1. [23] I. Osband, D. Russo, and B. Van Roy (2013) (More) efficient reinforcement learning via posterior sampling. In Advances in Neural Information Processing Systems, Vol. 26. External Links: Link Cited by: §2.2. [24] J. K. Pritchard, M. T. Seielstad, A. Pérez-Lezaun, and M. W. Feldman (1999) Population growth of human y chromosomes: a study of y chromosome microsatellites. Molecular Biology and Evolution 16 (12), p. 1791–1798. External Links: Document Cited by: §3.1. [25] T. G. Ritto, S. Beregi, and D. A. W. Barton (2022) Reinforcement learning and approximate bayesian computation for model selection and parameter calibration applied to a nonlinear dynamical system. Mechanical Systems and Signal Processing 181, p. 109485. External Links: ISSN 0888-3270, Document Cited by: §1. [26] C. P. Robert and G. Casella (2004) Monte carlo statistical methods. 2 edition, Springer Texts in Statistics, Springer, New York, NY. External Links: Document Cited by: §3. [27] D. S. Robertson, K. M. Lee, B. C. López-Kolkovska, and S. S. Villar (2023) Response-adaptive randomization in clinical trials: from myths to practical considerations. Statistical Science 38 (2), p. 185–208. External Links: Document Cited by: §4.1. [28] S. S. Roy, R. G. Everitt, C. P. Robert, and R. Dutta (2024) Generalized bayesian deep reinforcement learning. External Links: 2412.11743, Link Cited by: §1. [29] D. B. Rubin (1984) Bayesianly justifiable and relevant frequency calculations for the applied statistician. The Annals of Statistics 12 (4), p. 1151–1172. External Links: Document Cited by: §3.1. [30] S. A. Sisson, Y. Fan, and M. A. Beaumont (Eds.) (2018) Handbook of approximate bayesian computation. Chapman and Hall/CRC, Boca Raton, FL. External Links: Document Cited by: §1, §1, §3.1, §5. [31] M. Strens (2000) A bayesian framework for reinforcement learning. In Proceedings of the 17th International Conference on Machine Learning, p. 943–950. Cited by: §2.2. [32] M. Sunnåker, A. G. Busetto, E. Numminen, J. Corander, M. Foll, and C. Dessimoz (2013) Approximate bayesian computation. PLoS Computational Biology 9 (1), p. e1002803. External Links: Document Cited by: §3.1. [33] R. S. Sutton and A. G. Barto (2018) Reinforcement learning: an introduction. 2 edition, MIT Press, Cambridge, MA. Cited by: §1, §1, §2.1. [34] G. J. Székely and M. L. Rizzo (2004) Testing for equal distributions in high dimension. InterStat 5, p. 1249–1272. External Links: Link Cited by: §4.2.2. [35] S. Tavaré, D. J. Balding, R. C. Griffiths, and P. Donnelly (1997) Inferring coalescence times from dna sequence data. Genetics 145 (2), p. 505–518. External Links: Document Cited by: §3.1. [36] W. R. Thompson (1933) On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika 25 (3–4), p. 285–294. External Links: Document Cited by: §1. [37] C. Viscardi, M. Boreale, and F. Corradi (2021) Weighted approximate bayesian computation via sanov’s theorem. Computational Statistics 36 (4), p. 2719–2753. External Links: Document Cited by: §3.3.2. [38] I. Wilson, S. Julious, C. Yap, S. Todd, and M. Dimairo (2025) Response adaptive randomisation in clinical trials: current practice, gaps and future directions. Statistical Methods in Medical Research 34 (9), p. 1851–1874. External Links: Document Cited by: §1, §4.1. [39] A. Wong, T. Bäck, A. V. Kononova, and A. Plaat (2023) Deep multiagent reinforcement learning: challenges and directions. Artificial Intelligence Review 56, p. 5023–5056. External Links: Document Cited by: §1, §2.1. [40] B. D. Ziebart, A. Maas, J. A. Bagnell, and A. K. Dey (2008) Maximum entropy inverse reinforcement learning. In Proceedings of the AAAI Conference on Artificial Intelligence, p. 1433–1438. External Links: Document Cited by: §2.1. [41] B. D. Ziebart (2010) Modeling purposeful adaptive behavior with the principle of maximum causal entropy. Ph.D. Thesis, Carnegie Mellon University. External Links: Link Cited by: §2.2. Appendix A Supplementary Material for Section 3 A.1 Technical Details for ABC–SMC with M Simulations, Kernel Construction and Weight Simplification In [10], the authors show the advantage of considering M>1M>1 pseudo-datasets for each parameter proposal. In this case, the target distribution becomes p~ϵi(μ,∣x)∝p(μ)(∏m=1Mf(ym∣μ))(1M∑m=1Mκ(d(ym,x);ϵi)), p_ _i(μ, y x) p(μ) ( _m=1^Mf(y_m μ) ) ( 1M _m=1^Mκ (d(y_m,x); _i ) ), (16) where =(y1,…,yM) y=(y_1,…,y_M) denotes a vector of M independent realizations of the random variable Y. The corresponding approximated marginal posterior is p~ϵi(μ∣x) p_ _i(μ x) ∝∫Mp(μ)(∏m=1Mf(ym∣μ))(1M∑m=1Mκ(d(ym,x);ϵi)) _X^Mp(μ) ( _m=1^Mf(y_m μ) ) ( 1M _m=1^Mκ (d(y_m,x); _i ) )\,d y =p(μ)∫Mf(∣μ)(1M∑m=1Mκ(d(ym,x);ϵi)) =p(μ) _X^Mf( y μ) ( 1M _m=1^Mκ (d(y_m,x); _i ) )\,d y =p(μ)[κ(d(Y,x);ϵi)] =p(μ)\,E\! [κ (d(Y,x); _i ) ] =p(μ)Pr(d(Y,x)≤ϵi). =p(μ)\, \! (d(Y,x)≤ _i ). Hence, the likelihood function is still replaced by an estimate of Pr(d(Y,x)≤ϵi) \! (d(Y,x)≤ _i ). This estimator is less variable than the crude single-simulation version, although it is computationally more expensive, since each evaluation requires M pseudo-datasets. [10] also propose a specific choice for the forward and backward kernels. • The forward kernel Ki((μi−1,i−1),(μi,i))K_i (( _i-1, y_i-1),( _i, y_i) ) is a Markov kernel with invariant distribution p~ϵi(μ,∣x) p_ _i(μ, y x) and proposal distribution q((μi−1,i−1),(μi,i))=q(μi−1,μi)f(i∣μi).q (( _i-1, y_i-1),( _i, y_i) )=q( _i-1, _i)\,f( y_i _i). It is defined as Ki((μi−1,i−1),(μi,i)) K_i (( _i-1, y_i-1),( _i, y_i) ) =q(μi−1,μi)∏m=1Mf(yi,m∣μi)r((μi−1,μi),(i−1,i)),if p~ϵi(μi−1,i−1∣x)>p~ϵi(μi,i∣x),q(μi−1,μi)∏m=1Mf(yi,m∣μi),otherwise. = casesq( _i-1, _i) _m=1^Mf(y_i,m _i)\,r (( _i-1, _i),( y_i-1, y_i) ),\\[11.99998pt] p_ _i( _i-1, y_i-1 x)> p_ _i( _i, y_i x),\\[8.00003pt] q( _i-1, _i) _m=1^Mf(y_i,m _i), . cases where the acceptance rate is r((μi−1,μi),(i−1,i)) r (( _i-1, _i),( y_i-1, y_i) ) =p(μi)f(i∣μi)(1M∑m=1Mκ(d(yi,m,x);ϵi))q(μi,μi−1)f(i−1∣μi−1)p(μi−1)f(i−1∣μi−1)(1M∑m=1Mκ(d(yi−1,m,x);ϵi))q(μi−1,μi)f(i∣μi) = p( _i)\,f( y_i _i)\, ( 1M _m=1^Mκ(d(y_i,m,x); _i) )q( _i, _i-1)\,f( y_i-1 _i-1)p( _i-1)\,f( y_i-1 _i-1)\, ( 1M _m=1^Mκ(d(y_i-1,m,x); _i) )q( _i-1, _i)\,f( y_i _i) =p(μi)(1M∑m=1Mκ(d(yi,m,x);ϵi))q(μi,μi−1)p(μi−1)(1M∑m=1Mκ(d(yi−1,m,x);ϵi))q(μi−1,μi). = p( _i)\, ( 1M _m=1^Mκ(d(y_i,m,x); _i) )q( _i, _i-1)p( _i-1)\, ( 1M _m=1^Mκ(d(y_i-1,m,x); _i) )q( _i-1, _i). • The backward kernel Li−1((μi,i),(μi−1,i−1))L_i-1 (( _i, y_i),( _i-1, y_i-1) ) is defined as Li−1((μi,i),(μi−1,i−1)) L_i-1 (( _i, y_i),( _i-1, y_i-1) ) =p~ϵi(μi−1,i−1∣x)p~ϵi(μi,i∣x)Ki((μi−1,i−1),(μi,i)) = p_ _i( _i-1, y_i-1 x) p_ _i( _i, y_i x)\,K_i (( _i-1, y_i-1),( _i, y_i) ) =p(μi−1)∏m=1Mf(yi−1,m∣μi−1)(1M∑m=1Mκ(d(yi−1,m,x);ϵi))p(μi)∏m=1Mf(yi,m∣μi)(1M∑m=1Mκ(d(yi,m,x);ϵi)) = p( _i-1) _m=1^Mf(y_i-1,m _i-1) ( 1M _m=1^Mκ(d(y_i-1,m,x); _i) )p( _i) _m=1^Mf(y_i,m _i) ( 1M _m=1^Mκ(d(y_i,m,x); _i) ) ⋅qi(μi−1,μi)∏m=1Mf(yi,m∣μi)p(μi)(1M∑m=1Mκ(d(yi,m,x);ϵi))qi(μi,μi−1)p(μi−1)(1M∑m=1Mκ(d(yi−1,m,x);ϵi))qi(μi−1,μi) · q_i( _i-1, _i) _m=1^Mf(y_i,m _i)p( _i) ( 1M _m=1^Mκ(d(y_i,m,x); _i) )q_i( _i, _i-1)p( _i-1) ( 1M _m=1^Mκ(d(y_i-1,m,x); _i) )q_i( _i-1, _i) =(∏m=1Mf(yi−1,m∣μi−1))qi(μi,μi−1). = ( _m=1^Mf(y_i-1,m _i-1) )q_i( _i, _i-1). The importance weights become ωi(ℓ) _i^( ) =ωi−1(ℓ)⋅i(ℓ)i(ℓ), = _i-1^( )· N_i^( )D_i^( ), (17) where i(ℓ) _i^( ) =p(μi(ℓ))f(i(ℓ)∣μi(ℓ))(1M∑m=1Mκ(d(yi,m(ℓ),x);ϵi)) =p( _i^( ))\,f( y_i^( ) _i^( )) ( 1M _m=1^Mκ (d(y_i,m^( ),x); _i ) ) ×Li−1((μi(ℓ),i(ℓ)),(μi−1(ℓ),i−1(ℓ))), × L_i-1 (( _i^( ), y_i^( )),( _i-1^( ), y_i-1^( )) ), i(ℓ) _i^( ) =p(μi−1(ℓ))f(i−1(ℓ)∣μi−1(ℓ))(1M∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi−1)) =p( _i-1^( ))\,f( y_i-1^( ) _i-1^( )) ( 1M _m=1^Mκ (d(y_i-1,m^( ),x); _i-1 ) ) ×Ki((μi−1(ℓ),i−1(ℓ)),(μi(ℓ),i(ℓ))). × K_i (( _i-1^( ), y_i-1^( )),( _i^( ), y_i^( )) ). Since Li−1((μi(ℓ),i(ℓ)),(μi−1(ℓ),i−1(ℓ)))Ki((μi−1(ℓ),i−1(ℓ)),(μi(ℓ),i(ℓ))) L_i-1 (( _i^( ), y_i^( )),( _i-1^( ), y_i-1^( )) )K_i (( _i-1^( ), y_i-1^( )),( _i^( ), y_i^( )) ) =p~ϵi(μi−1(ℓ),i−1(ℓ)∣x)p~ϵi(μi(ℓ),i(ℓ)∣x) = p_ _i( _i-1^( ), y_i-1^( ) x) p_ _i( _i^( ), y_i^( ) x) =p(μi−1(ℓ))f(i−1(ℓ)∣μi−1(ℓ))(1M∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi))p(μi(ℓ))f(i(ℓ)∣μi(ℓ))(1M∑m=1Mκ(d(yi,m(ℓ),x);ϵi)), = p( _i-1^( ))f( y_i-1^( ) _i-1^( )) ( 1M _m=1^Mκ (d(y_i-1,m^( ),x); _i ) )p( _i^( ))f( y_i^( ) _i^( )) ( 1M _m=1^Mκ (d(y_i,m^( ),x); _i ) ), Eq. (17) simplifies to ωi(ℓ)=ωi−1(ℓ)⋅∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi)∑m=1Mκ(d(yi−1,m(ℓ),x);ϵi−1). _i^( )= _i-1^( )· _m=1^Mκ (d(y_i-1,m^( ),x); _i ) _m=1^Mκ (d(y_i-1,m^( ),x); _i-1 ). (18) A.2 Proofs Proof of Proposition 1. β(μ;ϵ,M) β(μ;ε,M) :=|f~ϵ,M(x∣μ)−fϵ,M(x∣μ)| := | f_ε,M(x μ)-f_ε,M(x μ) | =1M|∑m=1Mκ(d(ym,x);ϵ)−∑m=1MAϵ,x(ym)| = 1M |Σ _m=1^Mκ(d(y_m,x);ε)-Σ _m=1^MI_A_ε,x(y_m) | =1M|∑m=1MAϵ,x(ym)+∑m=1MAϵ,xc(ym)exp[−d(x,ym)ϵ2]−∑m=1MAϵ,x(ym)| = 1M |Σ _m=1^MI_A_ε,x(y_m)+Σ _m=1^MI_A^c_ε,x(y_m) [- d(x,y_m)ε^2 ]-Σ _m=1^MI_A_ε,x(y_m) | =1M∑m=1MAϵ,xc(ym)exp[−d(x,ym)ϵ2] = 1MΣ _m=1^MI_A^c_ε,x(y_m) [- d(x,y_m)ε^2 ] ≤∑m=1Mexp[−d(x,ym)ϵ2] ≤Σ _m=1^M [- d(x,y_m)ε^2 ] where Aϵ,xc:=y∈i:d(y,x)>ϵA^c_ε,x:=\y ^i:d(y,x)>ε\. ∎ A.3 LF-IBIS flowchart Figure 7: Overview of the proposed LF-IBIS algorithm. Starting from an initial approximation of the posterior distribution, the algorithm iteratively performs either an SMC-ABC step, in which the ABC tolerance is adaptively reduced, or an IBIS step, in which newly observed data are incorporated through likelihood-free importance reweighting. When necessary, particles are resampled and rejuvenated via an MCMC move. The resulting particles provide an approximation of the posterior distribution over model parameters and optimal policies. Appendix B Supplementary Material for Section 4 B.1 Beta proposal distribution In the clinical trial experiment, at iteration i, each new particle is sampled according to the Beta proposal distribution q(μ0,i−1(ℓ),⋅)q( _0,i-1^( ),·) defined as follow: μ0,i(ℓ)∼q(μ0,i−1(ℓ),⋅)=Beta(ai(ℓ),bi(ℓ)), _0,i^( ) q( _0,i-1^( ),·)=Beta (a_i^( ),b_i^( ) ), where ai(ℓ)a_i^( ) and bi(ℓ)b_i^( ) are computed via the method of moments as ai(ℓ) a_i^( ) =μ0,i−1(ℓ)(μ0,i−1(ℓ)(1−μ0,i−1(ℓ))2Var^i−1(μ0)−1), = _0,i-1^( ) ( _0,i-1^( )(1- _0,i-1^( ))2 Var_i-1( _0)-1 ), bi(ℓ) b_i^( ) =(1−μ0,i−1(ℓ))(μ0,i−1(ℓ)(1−μ0,i−1(ℓ))2Var^i−1(μ0)−1). = (1- _0,i-1^( ) ) ( _0,i-1^( )(1- _0,i-1^( ))2 Var_i-1( _0)-1 ). By construction, this proposal distribution satisfies [μ0,i(ℓ)] [ _0,i^( )] =ai(ℓ)ai(ℓ)+bi(ℓ)=μ0,i−1(ℓ), = a_i^( )a_i^( )+b_i^( )= _0,i-1^( ), Var[μ0,i(ℓ)] [ _0,i^( )] =ai(ℓ)bi(ℓ)(ai(ℓ)+bi(ℓ))2(ai(ℓ)+bi(ℓ)+1)=2Var^i−1(μ0). = a_i^( )b_i^( )(a_i^( )+b_i^( ))^2(a_i^( )+b_i^( )+1)=2 Var_i-1( _0). Where the empirical variance of the particles at iteration i−1i-1 is defined as Var^i−1(μ0)=1L−1∑ℓ=1L(μ0,i−1(ℓ)−μ¯0,i−1)2, Var_i-1( _0)= 1L-1 _ =1^L ( _0,i-1^( )- μ_0,i-1 )^2, where μ¯0,i−1=1L∑ℓ=1Lμ0,i−1(ℓ) μ_0,i-1= 1L _ =1^L _0,i-1^( ). Let Vi−1≔2Var^i−1(μ0).V_i-1 2\, Var_i-1( _0). When the empirical variance is too large, i.e., Vi−1≥μ0,i−1(ℓ)(1−μ0,i−1(ℓ))V_i-1≥ _0,i-1^( )(1- _0,i-1^( )), no unimodal Beta distribution exists with the desired moments. In this case, the proposal parameters are set such that ai(ℓ)+bi(ℓ)=1a_i^( )+b_i^( )=1, yielding a non-unimodal Beta distribution with substantial mass near the boundaries of the unit interval. To avoid this behavior, we enforce the following condition for the method-of-moments parameters: for each particle ℓ , ai(ℓ)=μ0,i−1(ℓ),if Vi−1≥μ0,i−1(ℓ)(1−μ0,i−1(ℓ)),μ0,i−1(ℓ)(μ0,i−1(ℓ)(1−μ0,i−1(ℓ))Vi−1−1),otherwise,a_i^( )= cases _0,i-1^( ),&if V_i-1≥ _0,i-1^( ) (1- _0,i-1^( ) ),\\[6.0pt] _0,i-1^( )\! ( _0,i-1^( ) (1- _0,i-1^( ) )V_i-1-1 ),&otherwise, cases bi(ℓ)=1−μ0,i−1(ℓ),if Vi−1≥μ0,i−1(ℓ)(1−μ0,i−1(ℓ)),(1−μ0,i−1(ℓ))(μ0,i−1(ℓ)(1−μ0,i−1(ℓ))Vi−1−1),otherwise.b_i^( )= cases1- _0,i-1^( ),&if V_i-1≥ _0,i-1^( ) (1- _0,i-1^( ) ),\\[6.0pt] (1- _0,i-1^( ))\! ( _0,i-1^( ) (1- _0,i-1^( ) )V_i-1-1 ),&otherwise. cases In the logistic reparameterized setting of the clinical application, we do not face the constraint of keeping the variance small in order to avoid negative parameters, as previously required for the Beta proposal on μ0 _0 and μ1 _1. Since β0,β1∈ℝ _0, _1 , the bivariate Normal proposal is always well defined, and scaling the covariance matrix does not introduce any feasibility issues. B.2 Additional Results: KL Divergence Analysis We compare the posterior distributions obtained in the experiments under the different inference strategies by computing the Kullback–Leibler (KL) divergence between the exact Bayesian posterior and its approximations produced by AR-ABC and by the two variants of LF-IBIS (ESS- and UP-based), for both experiments. For each method, we report the mean and standard deviation of the KL divergence computed over ten independent runs of the algorithm, all performed using the same observed history as in the experiments above. The KL divergence is reported only for the non-reparameterized versions of the experiments, since under the logistic reparameterization the exact Bayesian posterior is not available. All KL divergences are estimated from posterior samples using kernel density approximations. Table 4: Mean and standard deviation of the Kullback–Leibler divergence computed over ten independent runs of the algorithms with respect to the exact Bayesian posterior. Results are reported for inference on the parameters μ and on the policy π. Method Exp. 1 (Observation-based) Exp. 2 (Utility-based) μ π μ π LF-IBIS (ESS) 0.0954 (0.0132) 0.0632 (0.0101) 1.1589 (0.3567) 0.3034 (0.2010) LF-IBIS (UP) 0.1239 (0.0211) 0.0640 (0.0138) 0.6992 (0.2246) 0.1207 (0.0338) AR-ABC 0.0085 (0.0035) 0.0052 (0.0027) 1.1671 (0.0441) 0.9696 (0.0431)