domain: stat math

text:
Guillaume Kon

Kam King

Andrea Pandolfi

BrandDelta Marco Piretto

Matteo Ruggiero
email: mat-teo.ruggiero@unito.it

Approximate filtering via discrete dual processes

Keywords: Duality, Filtering, Bayesian inference, Diffusion, Particle filtering, Smoothing

We consider the task of filtering a dynamic parameter evolving as a diffusion process, given data collected at discrete times from a likelihood which is conjugate to the marginal law of the diffusion, when a generic dual process on a discrete state space is available. Recently, it was shown that duality with respect to a death-like process implies that the filtering distributions are finite mixtures, making exact filtering and smoothing feasible through recursive algorithms with polynomial complexity in the number of observations. Here we provide general results for the case of duality between the diffusion and a regular jump continuous-time Markov chain on a discrete state space, which typically leads to filtering distribution given by countable mixtures indexed by the dual process state space. We investigate the performance of several approximation strategies on two hidden Markov models driven by Cox-Ingersoll-Ross and Wright-Fisher diffusions, which admit duals of birth-and-death type, and compare them with the available exact strategies based on death-type duals and with bootstrap particle filtering on the diffusion state space as a general benchmark.

Introduction

Hidden Markov models are widely used statistical models for time series that assume an unobserved Markov process (X t ) t≥0 , or hidden signal, driving the process that generates the observations (Y t i ) i=0,...,n , e.g., by specifying the dynamics of one or more parameters of the observation density f Xt (y), called emission distribution. See [START_REF] Cappé | Inference in hidden Markov models[END_REF] for a general treatment of hidden Markov models. In this framework, the first task is to estimate the trajectory of the signal given noisy observations collected at discrete times 0 = t 0 < t 1 < • • • < t n = T , which amounts to performing sequential Bayesian inference by computing the so-called filtering distributions p(x t i |y t 0 , . . . , y t i-1 ), i.e., the marginal distributions of the signal at time t i conditional on observations collected up to time t i-1 . Originally motivated by real-time tracking and navigation systems and pioneered by [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF][START_REF] Kalman | New results in linear filtering and prediction theory[END_REF], classical and widely known explicit results for this problem include: the Kalman-Bucy filter, when both the signal and the observation process are formulated in a gaussian linear system; the Baum-Welch filter, when the signal has a finite state-space as the observations are categorical; the Wonham filter, when the signal has a finite state-space and the observations are Gaussian. These scenarios allow the derivation of so-called finite-dimensional filters, i.e., a sequence of filtering distributions whose explicit identification is obtained through a parameter update based on the collected observations and on the time separation between the collection times, such that the resulting computational cost increases linearly with the number of observation times. Other explicit results include [START_REF] Comte | Multiplicative Kalman filtering[END_REF][START_REF] Ferrante | On necessary conditions for the existence of finite-dimensional filters in discrete time[END_REF][START_REF] Ferrante | Finite dimensional filters for nonlinear stochastic difference equations with multiplicative noises[END_REF][START_REF] Genon-Catalot | A non-linear explicit filter[END_REF][START_REF] Genon-Catalot | Random scale perturbation of an AR(1) process and its properties as a nonlinear explicit filter[END_REF][START_REF] Runggaldier | Sufficient conditions for finite dimensionality of filters in discrete time: a Laplace transform-based approach[END_REF][START_REF] Sawitzki | Finite-dimensional filter systems in discrete time[END_REF]. Outside these classes, explicit solutions are difficult to obtain, and their derivation typically relies on ad hoc computations. This is especially true when the map x → f x is non-linear and when the signal transition kernel is known up to a series expansion, often intractable, as is the case for many widely used stochastic models. When exact solutions are not available, one must typically make use of approximate strategies, whose state of the art is most prominently based on extensions of the Kalman and particle filters. See, for example, [START_REF] Bain | Fundamentals of stochastic filtering[END_REF][START_REF] Chopin | An introduction to sequential Monte Carlo[END_REF].

A somewhat weaker but useful notion with respect to that of a finite-dimensional filter was formulated in [START_REF] Chaleyat-Maurel | Computable infinite-dimensional filters with applications to discretized diffusion processes[END_REF], who introduced the concept of computable filter. This extends the former class to a larger class of filters whose marginal distributions are finite mixtures of elementary kernels rather than single kernels. Unlike the former case, such a scenario entails a higher computational cost, usually polynomial in the number of observation times, but avoids the infinite-dimensionality typically implied by series expansion of the signal transition kernel. See [START_REF] Chaleyat-Maurel | Filtering the Wright-Fisher diffusion[END_REF].

Recently, [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] derived sufficient conditions for computable filtering based on duality. A dual process is a process D t which enjoys the identity

E[h(X t , d)|X 0 = x] = E[h(x, D t )|D 0 = d]. (1)
Here the expectation on the left-hand side is taken with respect to the transition law of the signal X t , and that on the right hand side with respect to that of D t , while the class of functions h(x, d) which satisfy the above identity are called duality functions. See [START_REF] Jansen | On the notion(s) of duality for Markov processes[END_REF] for a review and for the technical details we have overlooked here. Duality has received an increasing amount of attention recently, and has found applications in Markov process theory, interacting particle systems, statistical physics and population genetics, among other fields. See, e.g., [START_REF] Barbour | A transition function expansion for a diffusion model with selection[END_REF][START_REF] Carinci | Dualities in population genetics: A fresh look with new dualities[END_REF][START_REF] Ethier | Markov processes: characterization and convergence[END_REF][START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection[END_REF][START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection and the lambda coalescent limit[END_REF][START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF][START_REF] Giardinà | Duality and exact correlations for a model of heat conduction[END_REF][START_REF] Giardinà | Duality and hidden symmetries in interacting particle systems[END_REF][START_REF] Giardinà | Correlation inequalities for interacting particle systems with duality[END_REF][START_REF] Griffiths | Dual process in the two-parameter Poisson-Dirichlet Petrov diffusion[END_REF][START_REF] Möhle | The concept of duality and applications to Markov processes arising in neutral population genetics models[END_REF][START_REF] Ohkubo | Duality in interacting particle systems and boson representation[END_REF].

In the framework of duality, [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] showed that for a reversible signal whose marginal distributions are conjugate to the emission distribution (i.e., the Bayesian update at a fixed t can be computed in closed-form), computable filtering is guaranteed if the stochastic part of the dual process evolves on a finite state space. This in turn allows one to derive recursive formulae for the filtering distributions for non-linear hidden Markov models involving signals driven by Cox-Ingersoll-Ross (CIR) processes and K-dimensional Wright-Fisher (WF) processes. Along similar lines, duality was exploited for computable smoothing, whereby one conditions also on future data points, in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] and for nonparametric hidden Markov models driven by Fleming-Viot and Dawson-Watanabe diffusions in [START_REF] Ascolani | Predictive inference with Fleming-Viotdriven dependent Dirichlet processes[END_REF][START_REF] Ascolani | Smoothing distributions for conditional Fleming-Viot and Dawson-Watanabe diffusions[END_REF][START_REF] Papaspiliopoulos | Conjugacy properties of time-evolving Dirichlet and gamma random measures[END_REF].

In this paper, we investigate the impact on filtering problems for hidden Markov models when the dual process takes the more general form of a continuous-time Markov chain on a discrete state space. This is of interest for example in some population genetic models with selection [START_REF] Barbour | A transition function expansion for a diffusion model with selection[END_REF] or interaction [START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF] whose known dual processes are of birth-and-death (B&D) type, whose specific filtering problems are currently under investigation by some of the authors. When the dual process evolves in a countable state space, the filtering distribution can in general be expected to be countably infinite mixtures. This leads to inferential procedures which are not computable in the sense specified above, since the computation of the filtering distribution can no longer be exact. However, it is natural to wonder how the inferential procedures obtained in such a scenario perform, possibly aided by some suitable approximation strategies.

The paper is organized as follows. In Section 2 we identify sufficient conditions for filtering based on discrete duals and provide a general description of the filtering operations in this setting. In Section 3, we apply these results to devise practical algorithms which allow to evaluate in recursive form filtering and smoothing distributions under this formulation. Section 4 and 5 investigate hidden Markov models driven by Cox-Ingersoll-Ross and K-dimensional Wright-Fisher diffusions, and show that the first admits a dual given by a one-dimensional B&D process, and the second a dual given by a K-dimensional Moran model. Section 6 discusses several approximation strategies used to implement the above algorithms with these dual processes, and compares their performance with exact filtering based on the results in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] and with a bootstrap particle filter as a general benchmark. Finally, we conclude with some brief remarks.

Filtering via discrete dual processes

Let the hidden signal (X t ) t≥0 , which here takes the role of a temporally evolving target of estimation, be given by a diffusion process on

X ⊂ R K , for K ≥ 1. Observations Y t i ∈ Y ⊂ R D , D ≥ 1 are collected at discrete times 0 = t 0 < t 1 < • • • < t n = T with distribution Y t i ind ∼ f x (•), given X t i = x.
Given an observation Y = y, define the update operator ϕ y acting on measures ξ on X by

ϕ y (ξ)(x) := f x (y)ξ(x) µ ξ (y) , µ ξ (y) := X f x (y)ξ(x). (2)
Here the probability measure ξ acts as prior distribution which encodes the current knowledge on the signal X t , whereas µ ξ (y) is the marginal likelihood of a data point y when X t has distribution ξ. The update operator thus amounts to an application of Bayes' theorem for conditioning the probability measure ξ on a new observation y, leading to the updated measure ϕ y (ξ). Define also the propagation operator ψ t by

ψ t (ξ)(dx ′ ) := ξP t (dx ′ ) = X ξ(x)P t (x, dx ′ ). ( 3
)
where P t is the transition kernel of the signal. Hence ψ t (ξ) is the probability measure at time t obtained by propagating forward the law ξ of the signal at time 0 by means of the signal semigroup.

We will make three assumptions, the first two of which are the same as in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF].

Assumption 1 (Reversibility). The signal X t is reversible with respect to the density π, i.e., π(x)P t (x, dx ′ ) = π(x ′ )P t (x ′ , dx).

See the discussion in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] on the possibility of relaxing the above assumption. For K ∈ Z + , define now the space of multi-indices M = Z K + to be

M = {m = (m 1 , . . . , m K ) : m j ∈ Z + , for j = 1, . . . , K},
whose origin is denoted 0 = (0, . . . , 0).

Assumption 2 (Conjugacy). For Θ ⊂ R l , l ∈ Z + , let h : X × M × Θ → R + be such that sup x∈X h(x, m, θ) < ∞ for all m ∈ M, θ ∈ Θ and h(x, 0, θ ′ ) = 1 for some θ ′ ∈ Θ. Then f x (•) is conjugate to distributions in the family F = {g(x, m, θ) = h(x, m, θ)π(x), m ∈ M, θ ∈ Θ},
i.e., there exist functions t : Y ×M → M and T :

Y ×Θ → Θ such that if X t ∼ g(x, m, θ) and Y t |X t = x ∼ f x , we have X t |Y y = y ∼ g(x, t(y, m), T (y, θ)).
Here g(x, m, θ) takes the role of "current" prior distribution, i.e., the prior on the signal state which was possibly derived from previous conditioning and propagations, and g(x, t(y, m), T (y, θ)) takes the role of the posterior, i.e., g(x, m, θ) conditional on the new data point y. The functions t(y, m), T (y, θ) provide the transformations that update the parameters based on y. In absence of data, the condition h(x, 0, θ ′ ) = 1 reduces g(x, m, θ) to π(x).

The third assumption weakens Assumption 3 in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] by assuming the dual process has finite activity on a discrete state space, and possibly has a deterministic companion.

Assumption 3 (Duality). Given a deterministic process Θ t ∈ Θ and a regular jump continuous-time Markov chain M t on Z K + with transition probabilities

p m,n (t; θ) := P(M t = n|M 0 = m, Θ 0 = θ), (4)
equation ( 1) holds with D t = (M t , Θ t ) and h as in Assumption 2.

The following result provides a full description of the propagation and update steps which allow to compute the filtering distribution.

Proposition 2.1. Let Assumptions 1-3 hold, and let m∈M w m g(x, m, θ) be a countable mixture with m∈M w m = 1. Then, for ψ t as in (3) we have

ψ t m∈M w m g(x, m, θ) = n∈M w ′ n (t)g(x, n, Θ t ), (5)
where

w ′ n (t) = m∈M w n p m,n (t; θ), (6)
and p m,n (t; θ) are as in (4). Furthermore, for ϕ y as in (2), we have

ϕ y m∈M w m g(x, m, θ) = m∈M ŵn,θ (y)g(x, t(y, m), T (y, θ)) (7)
where ŵm,θ (y) ∝ w m µ m,θ (y) and

µ m,θ (y) := X f x (y)g(x, m, θ)dx. (8)
Proof. First observe that ψ t (g(x, m, θ)) = n∈M p m,n (t; θ)g(x, n, Θ t ), which follows similarly to Proposition 2.2 in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF]. Then the claim follows by linearity using the fact that

ψ t i≥1 w i ξ i = i≥1 w i ψ t (ξ i ) so that ψ t m∈M w m g(x, m, θ) = m∈M w m ψ t (g(x, m, θ)) = m∈M w m n∈M p m,n (t; θ)g(x, n, Θ t ) = n∈M m∈M
w m p m,n (t; θ)g(x, n, Θ t )

Using now the fact that

ϕ y i≥1 w i ξ i = i≥1 w i µ ξ i (y) j w j µ ξ j (y) ϕ y (ξ i ), (9)
we also find that

ϕ y m∈M w m g(x, m, θ) = m∈M w m µ m,θ (y) n∈M w n µ n,θ (y)
g(x, t(y, m), T (y, θ))

where µ m,θ (y) = X f x (y)g(x, m, θ)dx = X f x (y)h(x, m, θ)π(x)dx.
The expression (5), together with [START_REF] Barbour | A transition function expansion for a diffusion model with selection[END_REF], provides a general recipe on how to compute the forward propagation of the current marginal distribution of the signal g(x, m, θ), based on the transition probabilities of the dual continuous-time Markov chain. Since the update operator (2) can be easily applied to the resulting distribution, Proposition 2.1 then shows that under these assumptions all filtering distributions are countable mixtures of elementary kernels indexed by the state space of the dual process, with mixture weights determined by the dual process transition probabilities p m,n (t; θ). When the latter happens to give positive mass only to points {n ∈ M : n ≤ m}, as is the case for a pure-death process, then the right-hand side of (5) reduces to a finite sum, and one can construct an exact filter with a computational cost that is polynomial in the number of observations, as shown in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF].

The above approach can be seen as an alternative to deriving the filtering distribution of the signal by leveraging on a spectral expansion of the transition function P t in (3), which typically requires ad hoc computations and does not lend itself easily to explicit update operations through [START_REF] Ascolani | Predictive inference with Fleming-Viotdriven dependent Dirichlet processes[END_REF]. Note also that expressions like (5) can be used, by taking appropriate limits of p m,n (t; θ) as t → 0, to identify the transition kernel of the signal P t itself, see, e.g., [START_REF] Barbour | A transition function expansion for a diffusion model with selection[END_REF][START_REF] Griffiths | Dual process in the two-parameter Poisson-Dirichlet Petrov diffusion[END_REF][START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF].

3 Recursive formulae for filtering and smoothing

In order to translate Proposition 2.1 into practical recursive formulae for filtering and smoothing, we are going to assume for simplicity of exposition that the time intervals between successive data collections t i -t i-1 equal ∆ for all i. For ease of the reader, we will therefore use the symbol P ∆ instead of P t i -t i-1 for the signal transition function over the interval ∆ = t i -t i-1 . We will also use the established notation whereby i|0 : i -1 indicates that the argument refers to time t i = i∆, and we are conditioning on the data collected at times from 0 to t i-1 = (i -1)∆.

Define the filtering density

ν i|0:i (x i ) := p(x i |y 0:i ) ∝ X i p(x 0:i , y 0:i )dx 0:i-1 , (10)
i.e., the law of the signal at time t i given data up to time t i , obtained by integrating out the past trajectory. Define also the predictive density

ν i+1|0:i (x i ) := p(x i+1 |y 0:i ) = X p(x i |y 0:i )P ∆ (x i+1 |x i )dx i , (11)
i.e, the marginal density of the signal at time t i+1 , given data up to time t i . This can be expressed recursively as a function of the previous filtering density p(x i |y 0:i ), as displayed. Finally, define the marginal smoothing density

ν i|0:n (x i ) := p(x i |y 0:n ) ∝ X n p(x 0:n , y 0:n )dx 0:i-1 dx i+1:n , (12)
where the signal is evaluated at time t i conditional on all available data. The first two distributions above are natural objects of inferential interest, whereas the latter is typically used to improve previous estimates once additional data become available. Finally, for Θ ∆ as in Assumption 3 and t(•, •), T (•, •) as in Assumption 2, define for i = 0, . . . , n the quantities

ϑ i|0:i := T (y i , ϑ i|0:i-1 ), ϑ i|0:i-1 := Θ ∆ (ϑ i-1|0:i-1 ), ϑ 0|0:-1 := θ 0 . (13)
Here, ϑ i|0:i-1 denotes the state of the deterministic component of the dual process at time i, after the propagation from time i -1 and before updating with the datum collected at time i, and ϑ i|0:i the state after such update.

The following Corollary of Proposition 2.1 extends a result of [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] (see also Theorem 1 in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] for an easier comparison in a similar notation).

Corollary 3.1. Let Assumptions 1-3 hold, and assume that

ν i-1|0:i-1 (x) = m∈M w (i-1) m g(x, m, ϑ i-1|0:i-1 ).
Then (11) can be written, through (3), as

ν i|0:i-1 (x) = ψ ∆ (ν i-1|0:i-1 (x)) = m∈M w (i-1) ′ m g(x, m, ϑ i|0:i-1 ), w (i-1) ′ m = n∈M w (i-1) n p n,m (∆; ϑ i-1|0:i-1 ), m ∈ M, (14)
with p n,m (∆; ϑ i|0:i ) as in (4). Furthermore, given the observation y i , (10) can be written, through

, as

ν i|0:i (x) = ϕ y i (ν i|0:i-1 (x)) = m∈M w (i) m g(x, m, ϑ i|0:i ), w (i) m ∝ µ n,ϑ i|0:i-1 (y i )w (i-1) ′ n , m = t(y i , n), n ∈ M, (15)
with µ m,θ as in [START_REF] Cox | A theory of the term structure of interest rates[END_REF].

Algorithm 1 outlines the pseudo-code for implementing the update and propagation steps of Corollary 3.1. How to use these results efficiently can depend on the model at hand. When the transition probabilities p m,n (t; θ) are available in closed form, their use could lead to the best performance, but can also at times face numerical instability issues (as is the case pointed out in Section 4 below).

When the transition probabilities p m,n (t; θ) are not available in closed form, one can approximate them by simulating N replicates of the dual component M t , and then regroup probability masses according to the arrival states as done in [START_REF] Comte | Multiplicative Kalman filtering[END_REF]. The dual is typically easier to simulate than the original process, given its discrete state space. For instance, pure-death or B&D processes are easily simulated using a Gillespie algorithm [START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF], whereby one alternates sampling waiting times and jump transitions for the embedded chain. Depending on the dual process, there might also be more efficient simulation strategies.

A different type of approximation of the propagation step [START_REF] Comte | Multiplicative Kalman filtering[END_REF] in Corollary 3.1 can be based on pruning the transition probabilities or the arrival weights (cf. ( 14)) under a given threshold, followed by a renormalisation of the weights. Both this approximation strategy and that outlined above assign positive weights only to a finite subset of M, hence they overcome the infinite dimensionality

Algorithm 1: Filtering Input: Y 0:n , t 0:n Result: ϑ i|0:i , M i|0:i and W i = {w (i) m , m ∈ M i|0:i } Initialise Set ϑ 0|0 = T (Y 0 , θ 0 ) with T as in Assumption 2 Set M 0|0 = {t(Y 0 , 0)} = {m * } and W 0 = {1} with t as in Assumption 2 Compute ϑ 1|0 from ϑ 0|0 as in (13) Set M * = B(M 0|0 ) and W * = {p m * ,n (∆, ϑ 0|0 ), n ∈ M * } with p m,n as in (4) for i from 1 to n do Update Set ϑ i|0:i = T (Y i , ϑ i|0:i-1 ) Set W i = w * m µ m,ϑ i|0:i-1 (Yi) n∈M * w * n µ n,ϑ i (Yi) , m ∈ M * with µ m,θ defined as in (8) Set M i|0:i = {t(Y i , m), m ∈ M * }
and update the labels in W i Copy ϑ i|0:i , M i|0:i and W i to be reported as the output Propagation Compute ϑ i+1|0:i from ϑ i|0:i of the dual process state space. In the next sections we will investigate such strategies for two hidden Markov models driven by Cox-Ingersoll-Ross and K-dimensional Wright-Fisher diffusions.

Set M * = B(M i|0:i ) and W * = m∈M i|0:i w (i) m p m,n (∆, ϑ i|0:i ), n ∈ M * end Note: M i|0:i = {m ∈ M : w (i) m > 0} ⊂ M
Next, with the purpose of describing the marginal smoothing densities [START_REF] Chaleyat-Maurel | Filtering the Wright-Fisher diffusion[END_REF], we need an additional assumption and some further notation.

Assumption 4 For h as in Assumption 3, there exist functions d :

M 2 → M and e : Θ 2 → Θ such that for all x ∈ X , m, m ′ ∈ M, θ, θ ′ ∈ Θ h(x, m, θ)h(x, m ′ , θ ′ ) = C m,m ′ ,θ,θ ′ h(x, d(m, m ′ ), e(θ, θ ′ )), ( 16
)
where

C m,m ′ ,θ,θ ′ is constant in x.
Denote by

← ϑ i , ← ϑ ′
i the quantities defined in (13) computed backwards. Equivalently, these are computed as in [START_REF] Chopin | An introduction to sequential Monte Carlo[END_REF] with data in reverse order, i.e. using y n:0 in place of y 0:n , namely

← ϑ i|i+1:T = Θ ∆ ( ← ϑ i+1|i+1:T ), ← ϑ i|i:T = T (y i , ← ϑ i|i+1:T ), ← ϑ T |T = T (y T , θ 0 )
The following result extends Proposition 3 and Theorem 4 of [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF]: Proposition 3.2. Let Assumptions 1-4 hold, and let ν 0 = π. Then, for 0 ≤ i ≤ n -1, we have

p(x i |y 0:n ) = m∈M, n∈M w (i) m,n g(x i , d(m, n), e( ← ϑ i|i+1:n , ϑ i|0:i )), with w (i) m,n ∝ ← - w (i+1) m w (i) n C m,n, ← ϑ i|i+1:n ,ϑ i|0:i , ← ω (i+1) m = n∈M ← ω (i+2) n µ n, ← ϑ i+1|i+2:n (y i+1 )p t(y i+1 ,n),m (∆; ← ϑ i+1|i+1:n ) w (i)
n as in [START_REF] Crawford | Transition probabilities for general birth?death processes with applications in ecology, genetics, and evolution[END_REF] and

C m,n, ← ϑ i|i+1:n ,ϑ i|0:i
as in (16).

Proof. Note that Bayes' Theorem and conditional independence allow to write [START_REF] Chaleyat-Maurel | Filtering the Wright-Fisher diffusion[END_REF] as

ν i|0:n (x i ) = p(x i |y 0:n ) ∝ p(y i+1:n |x i )ν i|0:i (x i )
where the right-hand side involves the filtering distribution, available from Corollary 3.1, and the so called cost-to-go function p(y i+1:n |x i ) (sometimes called information filter), which is the likelihood of future observations given the current signal state. Along the same lines as Proposition 3 in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] we find that p(y i+1:

n |x i ) = m∈M ← ω (i+1) m h(x i , m, ← ϑ i|i+1:n ) with ← ω (i+1) m
as in the statement. The main claim can now be proved along the same lines as Theorem 4 in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF].

The main difference between the above result and Theorem 4 in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] lies in the fact that the support of the weights [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] is denoted by ← M i|i+1:n ) can possibly be countably infinite and coincide with the whole of M. Indeed, which points of M have positive weight are determined by the transition probabilities of the dual process, which in the present framework is no longer assumed to make only downward moves in M. Section 6 will deal with this possibly infinite support for a concrete implementation of the inferential strategy.

{ ← ω (i+1) m , m ∈ M} (which is

Cox-Ingersoll-Ross hidden Markov models

The Cox-Ingersoll-Ross diffusion, also known as the square-root process, is widely used in financial mathematics for modelling short-term interest rates and stochastic volatility. See [START_REF] Chen | Pricing interest rate options in a two-factor Cox-Ingersoll-Ross model of the term structure[END_REF][START_REF] Cox | A theory of the term structure of interest rates[END_REF][START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF][START_REF] Göing-Jaeschke | A survey and some generalizations of Bessel processes[END_REF][START_REF] Heston | A closed-form solution for options with stochastic volatility, with applications to bond and currency options[END_REF]. It also belongs to the class of continuous-state branching processes with immigration, arising as the large-population scaling limit of certain branching Markov chains [START_REF] Kawazu | Branching processes with immigration and related limit theorems[END_REF] and as the time-evolving total mass of a Dawson-Watanabe branching measure-valued diffusion [17]. Let X t be a CIR diffusion on R + that solves the one-dimensional SDE

dX t = δσ 2 -2γX t dt + 2σ X t dB t , X 0 ≥ 0, (17)
where δ, γ, σ > 0, which is reversible with respect to the Gamma density π = Ga(δ/2, γ/σ 2 ). The following proposition identifies a B&D process as dual to the CIR diffusion.

Proposition 4.1. Let X t be as in (17), let M t be a B&D process on Z + which jumps from m to m + 1 at rate λ m = 2σ 2 (δ/2 + m)(θ -γ/σ 2 ) and to m -1 at rate µ m = 2σ 2 θm, and let

h(x, m, θ) = Γ(δ/2) Γ(δ/2 + m) γ σ 2 -δ/2 θ δ/2+m x m e -(θ-γ/σ 2 )x . ( 18
)
Then (1) holds with D t = M t .

Proof. The infinitesimal generator associated to (17) is

Af (x) = (δσ 2 -2γx)f ′ (x) + 2σ 2 xf ′′ (x),
for f : R + → R vanishing at infinity. Letting h(x, m) denote ( 18) constant with respect to θ, a direct computation yields

Ah(•, m)(x) = (δσ 2 -2γx) mx m-1 -x m (θ -γ/σ 2 ) Γ(δ/2) Γ(δ/2 + m) γ σ 2 -δ/2 θ δ/2+m e -(θ-γ/σ 2 )x + 2σ 2 x m(m -1)x m-2 + x m (θ -γ/σ 2 ) 2 -2mx m-1 (θ -γ/σ 2 ) × Γ(δ/2) Γ(δ/2 + m) γ σ 2 -δ/2 θ δ/2+m e -(θ-γ/σ 2 )x = δσ 2 mθ δ/2 + m -1 h(x, m -1) + 2γ(θ -γ/σ 2 ) δ/2 + m θ h(x, m + 1) -[2γm + δσ 2 (θ -γ/σ 2 )]h(x, m) + 2σ 2 m(m -1) θ δ/2 + m -1 h(x, m -1) + 2σ 2 (θ -γ/σ 2 ) 2 δ/2 + m θ h(x, m + 1) -4σ 2 m(θ -γ/σ 2 )h(x, m) = 2σ 2 θmh(x, m -1) + 2σ 2 (δ/2 + m)(θ -γ/σ 2 )h(x, m + 1) -[2γm + σ 2 (δ + 4m)(θ -γ/σ 2 )]h(x, m).
where it can be checked that

2σ 2 θm + 2σ 2 (δ/2 + m)(θ -γ/σ 2 ) = 2γm + σ 2 (δ + 4m)(θ -γ/σ 2 ).
Hence the r.h.s. equals

Bg(m) = λ m [g(m + 1) -g(m)] + µ m [g(m -1) -g(m)] with g(•) := h(x, •), λ m = 2σ 2 (δ/2 + m)(θ -γ/σ 2 )
, and µ m = 2σ 2 θm, which is the infinitesimal generator of a B&D process with rates λ m , µ m . The claim now follows from Proposition 1.2 in [START_REF] Jansen | On the notion(s) of duality for Markov processes[END_REF].

Assign now prior ν 0 = Ga(δ/2, γ/σ 2 ) to X 0 , and assume Poisson observations are collected at equally spaced intervals of length ∆, specifically Y |X t = x iid ∼ Po(τ x), for some τ > 0. By the wellknown conjugacy to Gamma priors, we have X t |Y = y ∼ Ga(δ/2 + y, γ/σ 2 + τ ). For simplicity, and without loss of generality, we can set τ = 1, which allows to interpret the update of the gamma rate parameter as the size of the conditioning data set. The filtering algorithm starts by first updating the prior ν 0 for the signal to

ν 0|0 := ϕ Y 0 (ν 0 ). If we observe Y 0 = (Y 0,1 , . . . , Y 0,k ) at time 0, then ν 0|0 is the law of X 0 | k j=1 y 0,j = m ∼ Ga(δ/2 + m, γ/σ 2 + k).
Then ν 0|0 is propagated forward for a ∆ time interval, yielding ν 1|0 := ψ ∆ (ν 0|0 ). In light of Proposition 4.1, an application of (5) to ν 0|0 yields the infinite mixture

ψ ∆ Ga(δ/2 + m, γ/σ 2 + k) = n≥0 p m,n (∆)Ga(δ/2 + n, γ/σ 2 + k), (19)
where p m,n (t) are the transition probabilities of M t in Proposition 4.1. Hence, the law of the signal is indexed by integers, i.e., points of the dual state space, where after the first update the distribution gives mass one to the sum of the observations k j=1 y 0,j = m, whereas after the propagation the mass is spread over the whole Z + by the effect of the dual process. We then observe Y 1 ∼ f X 1 , which is used to update ν 1|0 to ν 1|1 and has the effect of shifting the probability masses of the mixture weights. For example, the weight p m,n (∆) in ( 19) is assigned to n ∈ Z + , but after the update based on

Y 1 = (Y 1,1 , . . . , Y 1,k ′ ) it will be assigned to n + m ′ if k ′ j=1 y 1,j = m ′
, on top of being transformed according to [START_REF] Cappé | Inference in hidden Markov models[END_REF]. We then propagate forward again and proceed analogously.

When the current distribution of the signal, after the update, is given by a mixture of type m∈Z + w m Ga(δ/2 + m, γ/σ 2 + k), it is enough to rearrange the mixture weights after the propagation step as in [START_REF] Barbour | A transition function expansion for a diffusion model with selection[END_REF].

The main difference with qualitatively similar equations found in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF] is now given by the transition probabilities p m,n (t) in [START_REF] Etheridge | Some mathematical models from population genetics[END_REF], which are those of the B&D process in Proposition 4.1. Before tackling the problem of how to use the above expressions for inference, we try to provide further intuition of the extent and implications of such differences. To this end, consider the simplified parameterization α = δ/2, β = γ/σ 2 , σ 2 = 1/2, τ = 1, whereby one can check that the embedded chain of the B&D process of Proposition 4.1 has jump probabilities

p m,m+1 = k(α + m) k(α + m) + m(β + k) , p m,m-1 = 1 -p m,m+1 .
Here m, k are the same as in the left-hand side of [START_REF] Etheridge | Some mathematical models from population genetics[END_REF], so m/k is the sample mean. It is easily verified that p m,m+1 < p m,m-1 if m/k > α/β and viceversa. Therefore, the dual evolves on Z + so that it reverts m/k to the prior mean α/β. Indeed, the dual has ergodic distribution NB (α, β/(β + k)), whose mean is kα/β, i.e., such that m/k on average coincides with α/β.

Recall now that the dual process elicited in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] for the CIR model is D t = (M t , Θ t ), with M t a pure-death process with rates from m to m -1 equal to 2σ 2 θ and Θ t a deterministic process that solves dΘ t /dt = -2σ 2 Θ t (Θ t -γ/σ 2 ), Θ 0 = θ. This dual has a single ergodic state given by (0, β) (note that [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] use a slightly different parameterization, where the ergodic state (0, β) means that, in the limit for t → ∞, the gamma parameters are the prior parameters). In particular this entails the convergence of p m,n (t) in [START_REF] Etheridge | Some mathematical models from population genetics[END_REF] to 1 if n = 0 and 0 elsewhere as t → ∞. Whence the strong ergodic convergence ψ t (g(x, m, θ)) → π as t → ∞, whereby the effect of the observed data become progressively negligible as t increases. One could then argue that in the long run, the filtering strategy based on the pure-death dual process in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] completely forgets the collected data.

Hence one could expect filtering with long-spaced observations (relative to the forward process autocorrelation) to be similar to using independent priors at each data collection point. On the other hand, the B&D dual can be thought as not forgetting but rather spreading around the probability masses in such a way as to preserve convergence of the empirical mean to the prior mean. It is not obvious a priori which of these two scenarios could be more beneficial in terms of filtering, hence in Section 6.1 we provide numerical experiments for comparing the performance of strategies based on these different duals.

In view of such experiments, note that the transition probabilities of the above B&D dual are in principle available in closed form (cf. [START_REF] Bailey | The elements of stochastic processes with applications to the natural sciences[END_REF][START_REF] Crawford | Transition probabilities for general birth?death processes with applications in ecology, genetics, and evolution[END_REF]), but their computation is prone to numerical instability. Alternatively, we can approximate the transition probabilities p m,n (t) in ( 19) by drawing N sample paths of the dual started in m and use the empirical distribution of the arrival points. This can in principle be done through the Gillespie algorithm [START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF], which alternates sampling waiting times and jumps of the embedded chain. A faster strategy can be achieved by writing the B&D rates in Proposition 4.1 as λ m = λm + β and µ m = µm with

λ = 2σ 2 (θ -γ/σ 2 ), β = σ 2 δ(θ -γ/σ 2 ), µ = 2σ 2 θ,
where λ, µ represent the per capita birth and death rate and β is the immigration rate. Then write M t = A t + B t where A t is the population size of the descendant of autochthonous individuals (already in the population at t = 0), and B t the descendants of the immigrants. These rates define a linear B&D process, whereby [START_REF] Tavaré | The linear birth-and-death process: an inferential retrospective[END_REF] suggests simulating A t by drawing, given

A 0 = i, F ∼ Bin(i, g(t)), A t ∼ NBin(F, h(t)) + F, (20)
with h(t) = (λ -µ)/(λ exp{(λ -µ)t} -µ) and g(t) = h(t) exp{(λ -µ)t}, with the convention NBin(0, p) = δ 0 . Let now N s be the number of immigrants up to time s, which follows a simple Poisson process with rate β, so given N t the arrival times are uniformly distributed on [0, t]. Once in the population, the lineage of each immigrating individual follows again a B&D process and can be simulated using [START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection[END_REF] starting at i = 1. Summing the numerosity of each immigrant family at time t yields B t .

Wright-Fisher hidden Markov models

The K-dimensional WF diffusion is a classical model in population genetics which is also widely used in many areas of applied probability and statistics, as it can model temporally evolving probability frequencies. It takes values in the simplex

∆ K = x ∈ [0, 1] K : 1≤i≤K x i = 1
and, in the population genetics interpretation, it models the temporal evolution of K proportions of types in an underlying large population. Its infinitesimal generator on C 2 (∆ K ) is

A = 1 2 K i,j=1 x i (δ ij -x j ) ∂ 2 ∂x i ∂x j + 1 2 K i=1 (α i -θx i ) ∂ ∂x i (21)
for α = (α 1 , . . . , α K ) ∈ R K + , θ = K i=1 α i , and its reversible measure is the Dirichlet distribution whose density with respect to Lebesgue measure is

π α (x) = Γ(θ) K i=1 Γ(α i ) x α 1 -1 1 • • • x α K -1 K , x K = 1 - K-1 i=1 x i .
See for example [START_REF] Ethier | Markov processes: characterization and convergence[END_REF], Chapter 10. The transition density of this model is (cf., e.g., [17], eqn. (1.27))

p t (x, x ′ ) = ∞ m=0 d m (t) m∈Z K + :|m|=m MN(m; m, x)π α+m (x ′ ), (22)
where MN(m; m, x) = m m 1 ,...,m K K i=1 x m i i , and where d m (t) are the transition probabilities of the block counting process of Kingman's coalescent on Z + , which has an entrance boundary at ∞. Cf., e.g., [17], eqn. (1.12).

It is well known that a version of Kingman's typed coalescent with mutation is dual to the WF diffusion. This can be seen as a death process on Z K + which jumps from m to me i at rate

q m,m-e i = m i (θ + |m| -1)/2. ( 23
)
See, for example, [START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection[END_REF][START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection and the lambda coalescent limit[END_REF][START_REF] Griffiths | Diffusion processes and coalescent trees[END_REF]. See also [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF], Section 3.3. The above death process with transitions d m (t) is indeed the process that counts the surviving blocks of the typed version without keeping track of which types have been removed.

Recall now that a Moran model with N individuals of K types is a particle process with overlapping generations whereby at discrete times a uniformly chosen individual is removed and another, uniformly chosen from the remaining individuals, produces one offspring of its own type, leaving the total population size constant. See, e.g., [START_REF] Etheridge | Some mathematical models from population genetics[END_REF]. In the presence of mutation, upon reproduction, the offspring can mutate to type j at parent-independent rate α j . The generator of such process on the set B(Z K + ) of bounded functions on Z K + can be written in terms of the multiplicities of types n ∈ Z K + as

Bf (n) = 1 2 1≤i̸ =j≤K n i (α j + n j )f (n -e i + e j ) - 1 2 1≤i̸ =j≤K n i (α j + n j )f (n), ( 24
)
where an individual of type i is removed at rate n i , the number of individuals of type i, and is replaced by an individual of type j at rate α j + n j .

The following proposition extends a result in [START_REF] Carinci | Dualities in population genetics: A fresh look with new dualities[END_REF] (cf. Section 5).

Proposition 5.1. Let X t have generator [START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection and the lambda coalescent limit[END_REF], let N t ∈ Z K + be a Moran model which from n jumps to ne i + e j at rate n i (α j + n j )/2, and let

h(x, n) = Γ(θ + |n|) Γ(θ) K i=1 Γ (α i ) Γ (α i + n i ) x n i i , θ = K i=1 α i .
Then (1) holds with D t = N t and h as above.

Proof. From ( 21), since θ = K i=1 α i , we can write

2A = 1≤i≤K x i (1 -x i ) ∂ 2 ∂x 2 i - 1≤i̸ =j≤K x i x j ∂ 2 ∂x i ∂x j + 1≤i≤K (α i (1 -x i ) -x i 1≤j≤K,j̸ =i α j ) ∂ ∂x i = 1≤i̸ =j≤K x i x j ∂ 2 ∂x 2 i - 1≤i̸ =j≤K x i x j ∂ 2 ∂x i ∂x j + 1≤i≤K α i 1≤j≤K,j̸ =i x j ∂ ∂x i - 1≤i̸ =j≤K α j ∂ ∂x i .
Then one can check that

2Ah(x, n) = 1≤i̸ =j≤K n i (α i + n i -1) Γ(θ + |n|) Γ(θ) x n-e i +e j K h=1 Γ (α h ) Γ (α h + n h ) - 1≤i̸ =j≤K n i (α j + n j )x n Γ(θ + |n|) Γ(θ) K h=1 Γ (α h ) Γ (α h + n h ) = 1≤i̸ =j≤K n i (α j + n j )h(x, n -e i + e j ) - 1≤i̸ =j≤K n i (α j + n j )h(x, n). Hence we have (Ah(•, n))(x) = (Bh(x, •))(n)
where the right hand side is (24) applied to h(x, n) as a function of n. The claim now follows from Proposition 1.2 in [START_REF] Jansen | On the notion(s) of duality for Markov processes[END_REF].

Assign now prior ν 0 = π α to X 0 , and assume categorical observations so that P(Y = j|X t = x) = x j . By the well-known conjugacy to Dirichlet priors, we have X t |Y = y ∼ π α+δy , where α + δ y = (α 1 , . . . , α j + 1, . . . , α K ) if y = j. When multiple categorical observations with vector of multiplicities m ∈ Z K + are collected, we write π α+m . The filtering algorithm then proceeds by first updating ν 0 to ν 0|0 := ϕ Y 0 (ν 0 ) = π α+m , if Y 0 = (Y 0,1 , . . . , Y 0,n ) yields multiplicities m, then propagating ν 0|0 to ν 1|0 := ψ ∆ (ν 0|0 ). In light of the previous result, an application of ( 5) to a single distribution g(x, m, θ) = π α+m yields the mixture

ψ t (π α+m ) = n:|n|=|m| p m,n (t)π α+n , (25)
where p m,n (t) are the transition probabilities of N t in Proposition 5.1. We then observe Y 1 |X 1 , which is in turn used to update ν 1|0 to ν 1|1 , as so forth. We refer again the reader to [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF], Section 2.4.2, for details on qualitatively similar recursive formulae.

In [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF], the overall multiplicity |n| equals the original |m|, as an effect of the population size preservation provided by the Moran model. The space {n : |n| = |m|} is finite, which shows that Assumption 3 need not lead to filtering distributions being infinite mixtures. However, it is not obvious a priori how ( 25) compares in terms of practical implementation with the different representation obtained in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF], namely

ψ t (π α+m ) = n:|n|≤|m| pm,n (t)π α+n ( 26
)
where pm,n (t) are the transition probabilities of the death process on Z K + with rates [START_REF] Ferrante | On necessary conditions for the existence of finite-dimensional filters in discrete time[END_REF]. Similarly to what has already been discussed for the CIR case, the death process dual has a single ergodic state given by the origin (0, . . . , 0), which entails the convergence of pm,n (t) to 1 if n = (0, . . . , 0) and 0 elsewhere, implying the strong convergence ψ t (π α+m ) → π α in [START_REF] Genon-Catalot | A non-linear explicit filter[END_REF]. This is ultimately determined by the fact that Kingman's coalescent removes lineages by coalescence and mutation until absorption to the empty set.

At first glance, a similar convergence is seemingly precluded to [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF]. However, we note in the first sum of (24) that the new particle's type is either resampled from the survived particles or drawn from the baseline distribution, in which case the new particle is of type j with (parent-independent) probability α j / K i=1 α j . Hence a heuristic argument is that each particle will be resampled from the baseline distribution in finite time. Together with the fact that K j=1 π α+δ j α j / K i=1 α i = π α , which follows from Corollary 1.1 in [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF], and considering that the number of particles is finite, we can therefore expect that, as t → ∞, we have the convergence ψ t (π α+m ) → π α in [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF] in this case as well.

The transition probabilities p m,n (t) in [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF], induced by the Moran model, are not available in closed form. This poses a limit on the direct applicability of the presented algorithms for numerical experiments. The first alternative is then to approximate them by drawing N points from the discrete distribution on the dual space before the propagation, making use of the Gillespie algorithm to draw as many paths, and evaluating the empirical distribution of the arrival points. Further alternatives are suggested by the fact that an appropriately rescaled version of the Moran model converges in distribution to a WF diffusion (see, e.g., [START_REF] Etheridge | Some mathematical models from population genetics[END_REF], Lemma 2.39). This would mean spatially rescaling the Moran model in [START_REF] Ferrante | Finite dimensional filters for nonlinear stochastic difference equations with multiplicative noises[END_REF] and writing the resulting generator as (cf. ( 24))

C |n| f (x) = |n| 2 1≤i̸ =j≤K n i |n| α j + n j |n| f x - e i |n| + e j /|n| -f (x) ,
where p i = n i /|n|, which by using a Taylor expansion can be shown to yield [START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection and the lambda coalescent limit[END_REF], when |n| → ∞.

We could therefore use the WF diffusion to approximate the Moran dual transitions in [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF]. This strategy in principle needs to deal with the intractable terms d m (t) in the transition function expansion [START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF] of the diffusion, but in this regard, one can adopt the solution found by [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF].

It is also known that one could also construct a sequence of WF discrete Markov chains indexed by the population size which, upon appropriate rescaling, converge to the desired WF diffusion (see, e.g., [START_REF] Karlin | A second course in stochastic processes[END_REF], Sec. 15.2.F or [START_REF] Etheridge | Some mathematical models from population genetics[END_REF], Sec 4.1). Since two sequences that converge to the same limit can to some extent be considered close to each other, one could then consider a WF discrete chain indexed by |n| with a parameterization that would make it converge to [START_REF] Etheridge | A coalescent dual process in a Moran model with genic selection and the lambda coalescent limit[END_REF], and use it to approximate the Moran transition probabilities. This would permit a straightforward implementation, given WF discrete chains have multinomial transitions.

Numerical experiments

To illustrate how the above results can be used in practice and how they perform in comparison with other methods, we are going to consider particle approximations of the dual processes to evaluate the predictive distributions for the signal, denoted p(x k+1 |y 1:k ). We compare these distributions with the exact predictive distribution obtained through the results in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] and those obtained through bootstrap particle filtering which make use of [START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF]. Particle filtering can be considered the state of the art for this type of inferential problems, a general reference being [START_REF] Chopin | An introduction to sequential Monte Carlo[END_REF]. The experiments in this section provide an argument in favour of using the proposed approximations for filtering and smoothing in light of the fact that these tasks can be successfully tackled once one obtains an accurate inference on the signal prediction distribution, as shown in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF].

We first briefly describe the specific particle approximation on the dual space we are going to use.

To approximate a predictive distribution ν i|0:i-1 (x i ), the classical particle approximation used in bootstrap particle filtering can be described as follows:

• sample X m i-1 iid ∼ ν i-1|0:i-1 ;
• propagate the particles by sampling

X m i ∼ p t (X m i-1 , •) for m = 1, . . . , N , with p t as in (22); • estimate ν i|0:i-1 with νi|0:i-1 := 1 N N m=1 δ X m i .
The filtering distributions obtained from the dual processes considered in this work are mixtures of distribution of the form ν i-1|0:i-1

(x i-1 ) = m w m h(x i-1 , m)π(x i-1
). We can use the approximation νi-1|0:i-1

(x i-1 ) := 1 N N n=1 h(x i-1 , m (n) )π(x i-1 )
, where m (n) ∼ m w m δ m , which amounts to transposing the previous approximation to the dual space or to perform a particle approximation of the discrete mixing measure. The natural approximation of ν i|0:i-1 (x i ) is therefore as follows:

• sample m (n) iid ∼ m w m δ m ;
• propagate the particles by sampling n (n) ∼ p m (n) ,• (t) where p m (n) ,• (t) are the transition prob-abilities of the dual process;

• estimate ν i|0:i-1 (x i ) with νi|0:i-1 (x i ) := 1 N N n=1 h(x i , n (n) )π(x i ).
Here some important remarks are in order. The above dual particle approximation is a finite mixture approximation of a mixture which can be either finite or infinite. Hence the above strategy can be applied both to filtering given death-like duals but also given general duals on discrete state spaces. The quality of the dual particle approximation, in general, may differ from that obtained through the particle filtering approximation since the particles live on a discrete space in the first case and on a continuous space in the second. This is the object of the following sections, at least for these specific examples. Finally, the ease of implementation of the two approximations may be very different because simulating from the original Markov process may be much harder than simulating from the dual process. An example is the simulation of Kingman's typed coalescent, immediate as compared to the simulation from [START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF], which would be unfeasible without [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF].

Cox-Ingersoll-Ross numerical experiments

The CIR diffusion admits two different duals:

• the death-like dual given by D t = (M t , Θ t ), with M t a pure-death process on Z + with rates 2σ 2 θm from m to m -1 and Θ t a deterministic process that solves the ODE dΘ t /dt = -2σ 2 Θ t (Θ t -γ/σ 2 ), Θ 0 = θ. Cf. [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF], Section 3.1.

• the B&D dual M t on Z + with birth rates from m to m+1 given by λ m = 2σ 2 (δ/2+m)(θ-γ/σ 2 ) and death rates from m to m -1 given by µ m = 2σ 2 θm respectively. Cf. Proposition 4.1.

Note that the latter is time-homogeneous, the former is not. In general, temporal homogeneity is to be preferred since a direct simulation with a Gillespie algorithm in the inhomogeneous case would require a time-rescaling. However, for this specific case, there is a convenient closed-form expression for the transition density of the first dual, which can be used to simulate for arbitrary time transitions (see the third displayed equation at page 2011 in [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF]). The second dual, by virtue of the temporal homogeneity, can be simulated directly using a Gillespie algorithm. This may be slow if the event rate becomes large, but as suggested in Section 4 we can see it as a linear B&D process, and a convenient closed-form expression can be used to simulate arbitrary time transitions.

We compare these two particle approximations with an exact computation of the predictive distribution following [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF] and to a bootstrap particle filtering approach on the original state space of the signal, which is easy to implement for arbitrary time transitions thanks to the Gamma-Poisson expansion of the CIR transition density (see details in [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF], Section 5).

Figure 1 shows the comparison of the above-illustrated strategies, with prediction performed for a forecast time horizon of 0.05. The CIR parameters were specified to δ = 11, σ = 1, γ = 1.1. The starting distribution for the prediction is a filtering distribution for a dataset whose last Poisson observation equals 4 (so the starting distribution is a mixture of Gamma distributions roughly centred around 4). The density estimates for the bootstrap particle filter were obtained from a Gamma kernel density estimator with bandwidth estimated by cross-validation. This is expected to induce a negligible error because the target distribution is a finite mixture of Gamma distributions.

The figure suggests that the bootstrap particle filter is slow to converge to the exact predictive distribution. However, with 50 particles, both dual approximations are already almost indistinguishable from the exact predictive distribution. This shows that accurately approximating the mixing measure on the discrete dual space seems to require fewer particles than approximating the continuous distribution on the original continuous state space. Next, we turn to investigating the error on the filtering distributions, which combines successive particle approximations. Since the update operation can be performed exactly through [START_REF] Chen | Pricing interest rate options in a two-factor Cox-Ingersoll-Ross model of the term structure[END_REF], particle filtering using the dual process is conveniently implemented like a bootstrap particle approximation to a Baum-Welch filter with systematic resampling. We quantify the error on the filtering distributions by measuring the absolute error on the first moment and the standard deviation of the filtering distributions (with respect to the exact computation). We also include the error on the signal retrieval, measured as the absolute difference between the first moment of the filtering distributions and the value of the simulated "true" hidden signal. The mean filtering error is averaged over the second half of the sequence of observations to avoid possible transient effects at the beginning of the observation sequence and further averaged over 50 different simulated datasets. The parameter specification is again δ = 11, σ = 1, γ = 1.1, with a single Poisson observation at each of 200 observation times, with interval spacing equal to 0.1. Figure 2 shows that the puredeath particle approximation performs better than the B&D particle approximation, but the latter performs comparably to the bootstrap particle filter, possibly with a modest advantage.

Wright-Fisher numerical experiments

The WF diffusion admits two different duals:

• Kingman's typed coalescent with mutation dual, given by a pure-death process on Z K + with rates λ m,m-e i = m i (|α| + |m| -1)/2 from m to me i . Cf. [START_REF] Papaspiliopoulos | Optimal filtering and the dual process[END_REF], Section 3.3.

• a Moran dual process, given a homogeneous B&D process on Z K + with rates λ m,m-e i +e j = m i (α j + m j )/2 from m to me i + e j . Cf. Proposition 5.1.

Here both processes are temporally homogeneous and can thus be easily simulated using a Gillespie algorithm, with the only caveat that the simulation can be inefficient when the infinitesimal rates are large. Similar to the CIR case, there is a closed-form expression for the transition probabilities in the first case, which can be used for simulation purposes for arbitrary time transitions (see Theorem 3.1 in [START_REF] Papaspiliopoulos | Conjugacy properties of time-evolving Dirichlet and gamma random measures[END_REF]). Although this expression is easy to handle in the one-dimensional CIR case, this is more challenging in the multi-dimensional WF case, with significant numerical stability issues raised by the need to compute the sum of alternated series with terms that can both overflow and underflow. In [START_REF] Kam King | Exact inference for a class of hidden Markov models on general state spaces[END_REF], these hurdles were addressed using arbitrary precision computation libraries and careful re-use of previous computations applicable when data is equally spaced. The Gillespie simulation strategy presents no such restriction and may be significantly faster when the event rates remain low.

As mentioned in Section (5), no closed-form expression is available for the Moran dual and the Gillespie algorithm approach is the main option, likely resulting is a slow algorithm. Alternatively, it is possible to approximate the Moran dual process by a finite population Wright-Fisher process, with the quality of approximation increasing with the population size. The interest of this approximation is that the event rate is lower for the finite population Wright-Fisher process than for the Moran process. This is related to the fact that weak convergence of a sequence of WF chains to a WF diffusion occurs when time is rescaled by a factor of N (cf. [START_REF] Karlin | A second course in stochastic processes[END_REF], Sec. 15.2.F), whereas a Moran model whose individual updates occur at the times of a Poisson process with rate 1, needs a rescaling by a factor N 2 to obtain a similar convergence. In other words, in order to establish weak convergence to the diffusion, time must be measured in units of N generations in the WF chain and in units of N 2 generations in the Moran model. For this reason, the resulting Gillespie simulation will be faster using a WF chain approximation to the Moran model.

The above considerations also suggest another possibility. Since the Moran process converges weakly to a Wright-Fisher diffusion, the latter could also be used as a possible approximation instead of a WF chain. In this case, it is possible to sample directly from [START_REF] Favero | A dual process for the coupled Wright-Fisher diffusion[END_REF] for arbitrary time transitions using the algorithm in [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF]. Hence we would be using a WF diffusion to approximate the dual Moran transitions in [START_REF] Frydman | Asymptotic inference for the parameters of a discrete-time square-root process[END_REF].

A standard bootstrap particle filter performed directly on the Wright-Fisher diffusion state space also crucially relies on the algorithm of [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF] for the prediction step, without which approximate sampling recipes from the transition density would be needed.

In Figure 3, we compare prediction strategies for a K = 4 WF diffusion using (figure legend indicated in parenthesis):

• the closed-form transition of the pure death dual ("Exact" in Fig. 3);

• an approximation of the above using a Gillespie algorithm ("PD");

• the Moran dual using a Gillespie algorithm ("BD Gillespie Moran");

• the WF chain approximation of the Moran dual using a Gillespie algorithm ("BD Gillespie WF");

• the WF diffusion approximation of the Moran dual using [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF] ("BD diffusion WF");

• a bootstrap particle filtering approximation using [START_REF] Jenkins | Exact simulation of the Wright-Fisher diffusion[END_REF] ("Bootstrap PF").

Figure 3 shows that among these particle approximations of p(x k+1 |y 1:k ), the Wright-Fisher diffusion approximation of the Moran dual seems to converge slowest, followed by the bootstrap particle filter. Prediction was performed for a forecast time horizon equal to 0.1, with WF parameters α = (1.1, 1.1, 1.1, 1.1). The starting distribution for the prediction is a filtering distribution for a dataset whose last multinomial observation is equal to (4, 0, 9, 2) (so the starting distribution is a mixture of Dirichlet distributions roughly centred around (4/15, 0, 9/15, 2/15)). The density estimates for the bootstrap particle filter are obtained from a Dirichlet kernel density estimator with bandwidth estimated by cross-validation. This is expected to induce a negligible error because the target distribution is a finite mixture of Dirichlet distributions. The figure shows that the particle approximation of the pure death dual process using the closedform transition exhibits better performance. The bootstrap particle approximation has the fastest improvement relative to increasing the number of particles. Overall, the Moran dual performs better or comparably to bootstrap particle filtering.

Concluding remarks

We have provided conditions for filtering diffusion processes on R K which avoid computations on the uncountable state space of the forward process when a dual process is available and is a regular jump continuous-time Markov chain on a discrete state space. Motivated by certain diffusion models for which only duals with a countable state space are known (e.g. B&D-like duals for WF diffusions with selection), we have investigated the performance of filtering based on a B&D dual for the CIR diffusion and based on a Moran process dual for the WF diffusion. All approximation methods proposed appear to be valuable strategies, despite resting on different simulation schemes. The optimal strategy is bound to depend on the application at hand, together with several other details like the interval lengths between data collection times, and possibly be constrained by which of these tools are available. Overall, approximate filtering using B&D-like duals may perform better or comparably to bootstrap particle filtering, with the advantage of operating on a discrete state space. The computational effort for each of these strategies is also bound to depend on a series of factors the identification of which is beyond the scope of this contribution.

  is the support of the weights of ν i|0:i ; B(m) denotes the states reached by the dual process from m, and B(M) those reached from all m ∈ M.

Figure 1 :

 1 Figure 1: Comparison of the signal predictive distribution p(x k+1 |y 1:k ) obtained through the approximation approach to the death-process dual and the B&D dual, and through the bootstrap particle filter, with the exact predictive. The number of particles used for the approximations are 50, 100, 500, 1000 and are indicated in the panel labels.

Figure 2 :

 2 Figure 2: Mean filtering error as a function of the number of particles for the various particle approximation methods.

Figure 3 :

 3 Figure 3: Convergence of the WF predictive distribution (only the first dimension) with the number of particles for the various particle approximations.

Figure 4

 4 Figure4evaluates the filtering error for a WF process with K = 3 and parameters α = (1.1, 1.1, 1.1), given 20 categorical observations collected at each time, and 10 collection times with a time interval of 1. We consider increasing numbers of particles and use 100 replications to estimate the error. The figure shows that the particle approximation of the pure death dual process using the closedform transition exhibits better performance. The bootstrap particle approximation has the fastest improvement relative to increasing the number of particles. Overall, the Moran dual performs better or comparably to bootstrap particle filtering.

Figure 4 :

 4 Figure 4: Mean filtering error as a function of the number of particles for the various particle approximation methods.


Kevin Guilaume
email: kevin.guilaume@enedis.fr

Bertrand Lacombe
email: bertrand.lacombe@enedis.fr

Bastien Lacroix
email: bastien.lacroix@enedis.fr

Jean-Louis Dirion
email: dirion@mines-albi.fr

INVESTMENT EFFICIENCY ASSESSMENT ON THE ELECTRIC GRID BY GROUP OF CONTROL

The French DSO Enedis spends several hundred million euros every year to maintain and renew its existing network.

The investments are made according some policies. One of these policies, called Refurbishment, consists in renovating some segment lines of the overhead medium voltage power network.

The choices of the segment lines are made by some experts according mainly to the past outages.

To be sure that the choice has been appropriate, we would like to know if the investment is efficient or not. To develop an investment efficiency assessment, we have to compare the number of outages after treatment in a group of refurbished segment lines with a group of un-refurbished segment lines all other things being equal.

The method proposed in this paper consist of defining appropriate weights for the un-refurbished segment lines group to make it statistically comparable to the refurbished segment lines group.

.

When we apply the appropriate weights to the unrefurbished segment lines group, we obtain a control group that can be compared to the refurbished segment lines group. By counting the difference of the outages after treatment between the refurbished segment lines group and the control group, we can estimate the performance of the investment.

Using the statistical test "Exact Test of Fisher" [2], we can prove that the difference of outages is significant.

INTRODUCTION

With a network of 1.4 million km, Enedis, main DSO in France, spends several hundred million euros every year to maintain and renew its existing network in order to maintain power supply quality.

Business Units follow different policies depending on the type of asset considered. One of these policies, called Refurbishment, consists of renovating some segment lines of the overhead medium voltage power network.

A segment line is the smallest part of the grid with the same features (material, age, size of the cables…). Some experts carry out surveys to choose the segment lines on which we expect to maximize the improvement of electricity supply continuity by analysing, for example, past outages. When the choices are made, the intervention staff performs the renovation on the grid.

To be sure that the choice has been appropriate, we would like to know if the investment is efficient or not. Furthermore, Enedis uses the information about the efficiency of the investments made in the dialog with the public stakeholders who own the managed network.

Developing an investment efficiency assessment is difficult, as we have to compare the number of outages after treatment in a group of refurbished segment lines with a group of un-refurbished segment lines all other things being equal. Groups have to be comparable on a spatial and a temporal aspect.

First, let's look at the method used to define appropriate weights for the un-refurbished segment lines group. The method, called calibration estimator, is widely used in survey sampling theory to reweight the sample and make it representative of the population studied.

Then, we will see how to use a statistical test to compare the weighted un-refurbished segment lines group, called control group with the refurbished segment lines.

CALIBRATION ESTIMATOR Theory and vocabulary

The principle can be summarized in the next picture [Figure 1] where we reweight the un-refurbished segment lines group among two variables called control variables (VarA and VarB) in order to have the same margin (total weighted or frequency weighted of VarA and VarB) between the refurbished segment lines group, i.e. target group and the un-refurbished weighted segment lines group, i.e. control group.

Time and spatial aspects between the target and the control group

In order to count the outages among same size time periods, we would like to be sure to have comparable time intervals between the refurbished segment lines group and the weighted un-refurbished segment lines group before and after the refurbishment date.

As we don't have refurbishment date for the unrefurbished segment lines group, we will create a "fictive date" as it is shown in the below figure [Figure 2]. To deal with time aspects, we consider the number of outages before refurbishment or fictive date as a control variable: in [Figure 2], it means we make sure that on the red time intervals we have the same outage rate for both lines groups.

To deal with spatial aspects, we will add the department (the region) in the control variables.

Choices of control variables

Choosing appropriately control variables is essential to construct a control group representative of the target group.

In addition of the number of outages before refurbishment or fictive date and department, the technical experts of Enedis consider size, age and designation of the segment lines, structure and number of clients of the grid as control variables.

RESULT AND STATISTICAL TEST Contingency table

When the control group is well defined, we can compare the target group and the control group through a contingency table [Figure 3]. In this contingency table, in order to assess the efficiency of the refurbishment on 6,161 segment lines, we used the calibration estimator to affect weights to the 61,178 unrefurbished segment lines. Thus, we obtain a control group of 6,161 segment lines comparables to the 6,161 refurbished segment lines.

We can see that we have two time less outages on the refurbished segment lines than on the un-refurbished segment lines, all other things being equal considering the control variables.

In the next section we will check if the result is statistically efficient.

Exact Fisher test

To assess if the refurbishment is statistically efficient, we can compare the target group and the control group by doing statistical test like the Fisher test [2] widely used to test the link of two variables through a contingency table.

In the following screenshot [Figure 4], we can see the R output of the fisher test function: As the pvalue of the test is less than 0.05, we can reject the hypothesis that the number of incidents and the groups (control and target) are independent, so, the refurbishment is statistically efficient.

The output of the R fisher test function also gives a confidence interval of the odds ratio that is the inverse of the ratio (0.5 for 2 in our case).

We can conclude that the refurbishment is statistically efficient.

CONCLUSION

The algorithm proposed in this article is useful to prove the efficiency of our investments on the overhead medium voltage power network all other things being equal.

The efficiency is statistically proven. This is well appreciated by the public stakeholders in the region to whom Enedis reports regarding the investments made on the network.

We used it for one of ours policies but it can be used for others policies.

Our further actions will be to think about some declinations of this algorithm by cluster of segment lines to find the cluster which maximise the efficiency of our policy.

MISCELLANEOUS

Figure 1 :

 1 Figure 1: Calibration estimator principle

Figure 2 :

 2 Figure 2: Fictive date All the un-refurbished segment lines are grouped by age and given a fictive date which is the average of the refurbishment date for the same age refurbished lines.

Figure 3 :

 3 Figure 3: Example of result on a geographical subdivision of Enedis

Figure 4 :

 4 Figure 4 : Fisher's Exact test -R output

Acknowledgments

Thanks to IMT Mines Albi School and particularly Jessy IDEZ, Gregory FOLTYN, Florentin ROUBY for having implemented the algorithm in R.

Thanks to Isabel GARCIA BURREL for his trust and support.


Sophie Fortz
email: sophie.fortz@unamur.be

Xavier Devroey
email: xavier.devroey@unamur.be

Towards Feature-based ML-enabled Behaviour Location

Keywords: Variability Exploration, Feature, Machine Learning, Software Variability

Mapping behaviours to the features they relate to is a prerequisite for variability-intensive systems (VIS) reverse engineering. Manually providing this whole mapping is labour-intensive. In black-box scenarios, only execution traces are available (e.g., process mining). In our previous work, we successfully experimented with variantbased mapping using supervised machine learning (ML) to identify the variants responsible of the production of a given execution trace, and demonstrated that recurrent neural networks (RNNs) work well (≥ 80% accuracy) when trained on datasets in which we label execution traces with variants. However, this mapping (i) may not scale to large VIS because of combinatorial explosion and (ii) makes the internal ML representation hard to understand. In this short paper, we discuss the design of a novel approach: feature-based mapping learning.

INTRODUCTION

Variability-intensive systems (VIS) constitute a unique category within the realm of software systems, distinguished by the ability to adapt their characteristics and behaviour through the activation or deactivation of specific features. These systems encompass a diverse range of applications, including Software Product Lines (SPLs) [START_REF] Apel | Feature-Oriented Software Product Lines -Concepts and Implementation[END_REF][START_REF] Pohl | Software Product Line Engineering: Foundations, Principles and Techniques[END_REF], and business process families [START_REF] Marcello | Business process variability modeling: A survey[END_REF]. Yet, their behaviours are complex to analyse primarily due to combinatorial explosion and the absence of explicit variability-aware models [START_REF] Halin | Test them all, is it worth it? Assessing configuration sampling on the JHipster Web development stack[END_REF]. In this context, locating variations is essential for any reverse-engineering endeavour [START_REF] Klewerton | Reengineering legacy applications into software product lines: a systematic mapping[END_REF]. Black-box testing techniques can also benefit from this information to e.g., sample which variants should be tested first [START_REF] Halin | Test them all, is it worth it? Assessing configuration sampling on the JHipster Web development stack[END_REF].

We previously proposed VaryMinions [START_REF] Fortz | VaryMinions: leveraging RNNs to identify variants in event logs[END_REF] 1 to tackle the research question of identifying variant(s) that may have produced a given trace. VaryMinions relies on Recurrent Neural Networks (RNNs) mapping the occurrence of event messages to variants directly. Once the RNNs are trained, the learned mappings can be applied to new event logs. VaryMinions showed good performance, with a minimum of 80% accuracy in retrieving which variants could produce a specific event log.

Currently, to determine if a variant can generate a specific trace, we have to observe a substantial number of traces per variant to establish statistically significant patterns and draw confident conclusions. In practice, performing predictions across all potential variants requires enumerating and executing each variant, preferably multiple times. While this is doable for event logs associated with a limited number of variants [START_REF] Fortz | VaryMinions: leveraging RNNs to identify variants in event logs[END_REF], we may face combinatorial explosion while enumerating variants. Sampling is a common practice in such cases [START_REF] Halin | Test them all, is it worth it? Assessing configuration sampling on the JHipster Web development stack[END_REF], though it may miss important configurations impacting the ML model confidence [START_REF] Temple | Empirical assessment of generating adversarial configurations for software product lines[END_REF]. Mapping a trace to the entire configuration of a variant is quite demanding in terms of ML algorithms. We chose our Deep Learning (DL) models for their capacity to treat sequences. Yet, as in any DL model, interpretability is an issue, ultimately questioning the notion of learned features in the model's representation (latent space) [START_REF] Temple | Explicit or Implicit? On Feature Engineering for ML-Based Variability-Intensive Systems[END_REF]. Additionally, the cost of training and tuning these models is not negligible.

FROM VARIANTS TO FEATURES

To overcome these limitations, we shift our perspective from predicting which variants generate a specific execution trace to assessing the significance of (de)selecting particular features regarding that trace. We propose to independently predict the (de)selection of each feature, relying solely on the execution trace and the ground truth, represented as a vector of features, which can produce the trace. Each feature could be predicted using smaller and more interpretable ML models. But we should choose models that can output their confidence in their decision so that each feature can have three values: selected, deselected, or undetermined (i.e., irrelevant in the current context).

Consequently, we need to train a set of small ML models, one for each feature. These models associate the presence of specific messages with the (de)selection of a given feature. As illustrated in the left segment of Figure 1, we have chosen Support Vector Machines (SVMs) as our ML models, although other options are worth considering. For instance, the first SVM learns the (de)selection of the first feature, the second SVM handles the second one, and so on. Once these models are trained, they can be applied to new execution traces. However, rather than providing a binary decision (selected or deselected), we intend to leverage the confidence level associated with the decision, represented as a floating value between -1 and 1, with 0 being the most uncertain case (i.e., lying perfectly on the separation function). This more expressive approach justifies our choice of SVMs, as depicted in Figure 1. Based on this confidence score, we can establish thresholds for decisionmaking. A feature will be considered as selected if the confidence score falls between, for instance, 0.65 and 1, while it will be regarded as not selected if the confidence score ranges between -0.65 and -1. Any confidence score between -0.65 and 0.65 is deemed uncertain due to its proximity to the decision boundary, resulting in an "undetermined" decision. It's worth noting that the definition of these thresholds and their implications will require thorough evaluation.

The previous description, which includes the three possible outputs (selected, deselected, and "undetermined"), pertains to the prediction phase. Each model yields a confidence score for the (de)selection of its corresponding feature. Consequently, we obtain potential partial configurations (i.e., when at least one feature is predicted as "undetermined"). In some scenarios, it may be beneficial to get a complete configuration, achieved by replacing all the "undetermined" with random (de)selected values. However, independent learning of (de)selection for features can introduce inconsistencies with the underlying feature model of the system, which is a problem that also arises when attempting to set values for "undetermined" features. Consequently, we need to use a solver to check that all constraints are satisfied. If any constraints are violated, we must develop strategies to rectify the configuration, such as through repair methods [START_REF] Buchner | Automated inconsistency detection and solution proposals in cyber-physical system networks[END_REF][START_REF] Franz | Configfix: interactive configuration conflict resolution for the linux kernel[END_REF][START_REF] Jahanbanifar | Implicit and Explicit Semantics Integration in Proof-Based Developments of Discrete Systems: Communications of NII Shonan Meetings[END_REF][START_REF] Xiong | Range fixes: Interactive error resolution for software configuration[END_REF]. As illustrated in the right part of Figure 1, this repair process should consider the input trace to ensure that the output of the repair procedure remains capable of executing the trace. Nonetheless, the initial output partial configuration provides the "form" (i.e., the selection and deselection of features) of a configuration, or set of configurations, having a high likelihood of generating the input execution traces. In case of errors or unexpected behaviour, we can swiftly focus on the most relevant features and their interactions.

CONCLUSION

In this paper, we highlighted the role of mapping choices when using machine learning to relate configurations with their execution trace. While our initial direct mapping performed well (i.e., variantbased), it is challenging in terms of scalability, interpretability and computation. Learning a feature-to-trace mapping may alleviate these issues, but we must assess the proposed correction mechanism. We hypothesise that there is no perfect mapping and that we could eventually revisit the product line analysis cube [START_REF] Alexander Von Rhein | The PLA Model: On the Combination of Product-Line Analyses[END_REF] for ML tasks. Our long-term goal is to offer a detailed account of the trade-offs involved in the choices along the cube's edges and to help the community make informed mapping choices.

Figure 1 :

 1 Figure 1: Conceptual representation of the feature-based approach during the prediction phase.

VAMOS'24, February 05-07, 2024, Bern, Switzerland 2024. ACM ISBN 978-x-xxxx-xxxx-x/YY/MM. . . $15.00 https://doi.org/10.1145/nnnnnnn.nnnnnnn

the replication package is available at[START_REF] Fortz | VaryMinions[END_REF]

ACKNOWLEDGEMENTS

Sophie Fortz is supported by the Fonds de la recherche scientifique (FRS-FNRS) via a FRIA grant. Gilles Perrouin is an FRS-FNRS Research Associate. This research is partly supported by the EOS VeriLearn project under FRS-FNRS Grant No. O05518F-RG03, and the CyberExcellence project (No. 2110186) funded by the Public Service of Wallonia (SPW Recherche).