Patrick Pamphile
email: patrick.pamphile@univeriste-paris-saclay.fr

Gilles Celeux
email: gilles.celeux@inria.fr

Bayesian Estimations for Weibull Competing Risk Model with Masked Causes and Heavily Censored Data

Keywords: Competing risks, Weibull distribution, Masked causes, Censored data, EM algorithms, Bayesian sampling, BRM algorithm

In a reliability or maintenance analysis of a complex system, it is important to be able to identify the main causes of failure. Therefore a Weibull competing risks model is generally used. However, in this framework estimating the model parameters is a difficult ill-posed problem. Indeed, the cause of the system failure may not be identified and the experiment may be censored by the duration of the study. In addition, the other causes are naturally censored by the first one. In this paper, we propose a new method for estimating the parameters of the Weibull competing risks model, with masked causes and heavily censored data. We use a restoration of missing data through a Bayesian sampling of parameters with a weakly informative prior distribution. The mean of the posterior distribution can thus estimated by importance sampling. The proposed method is not an iterative method and therefore can be parallelized.

Experiments based on simulated data and a reliability data set show that the prediction performance of the proposed method is superior to the EM algorithm and its variants, as well as the Gibbs sampler, for low to very heavy censoring rates and even in the case of weakly separated components.

Introduction

Competing risk models are widely used to model lifetimes in medical or industrial applications. Indeed, they correspond to the lifetime of a patient or a system with multiple causes of failure. In order to have an effective diagnosis, it is essential to identify the distributions related to main causes of failures. Unfortunately, when a failure occurs, the cause is often unidentified. In fact, exact diagnosis of the failure cause could be impossible or too resource-consuming to be performed on every failure. In such cases, the modeler is faced with a competing risk model with hidden causes of failure (see among others [START_REF] Crowder | Classical competing risks[END_REF]).

Statistical inference for competing risks models has been widely discussed by many authors (see, for example [2, 3] or [START_REF] Sarhan | Statistical analysis of competing risks models[END_REF]). In these references, different competing failure time distributions such as exponential, gamma and Weibull are considered. When the causes of the failures are hidden, one is faced with an ill-posed problem and the standard maximum likelihood method is no longer satisfactory. The EM algorithm or its stochastic variant S-EM are then used [5]. But, In the case of real data, it is usual to have heavily censored (e.g. 70%

of censored data). The EM and S-EM algorithms are then no longer effective in fitting a competing risk model with masked causes to very heavily censored data. The Bayesian method is then an alternative solution cf. [START_REF] Berger | Bayesian analysis for the poly-Weibull distribution[END_REF][START_REF] Ranjan | Classical and Bayesian estimation for the parameters of a competing risk model based on minimum of exponential and gamma failures[END_REF]8,9] or [START_REF] Azizi | Statistical inference for competing risks model under progressive interval censored weibull data[END_REF].

In addition, the case of competing Weibull distributions is widely used. Indeed, Weibull distributions exhibit a decreasing, constant or increasing hazard function, which makes them suitable for modeling complex failure data (see for instance [START_REF] Murthy | Weibull Models[END_REF] and reference therein). Unfortunately, the estimator of the shape parameters has no closed-form expression : it is necessary to use iterative numerical approximations sensitive to the starting point and the censoring rate. For Bayesian methods, there is no conjugated prior distribution for Weibull distribution when the two parameters are unknown. The estimation of the posterior distribution or its moments requires, here again, approximations by iterative numerical methods. The Gibbs algorithm then converges very slowly cf. [START_REF] Casella | Explaining the Gibbs sampler[END_REF].

In parametric estimation with missing data, Bacha and al [START_REF] Bacha | Estimation de durées de vie fortement censurées[END_REF] have proposed the Bayesian restoration maximization (BRM) method. It is a Bayesian sampling -importance -resampling method [START_REF] Gelfand | Sampling-based approaches to calculating marginal densities[END_REF]: the parameters are simulated from the prior distribution and the missing data are then restored. The expectation of the posterior distribution is then estimated by an importance sampling technic [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. One of the traps of importance sampling is to use a proposal distribution too far from target distribution. Therefore, Rubin suggests a re-sampling cf. [START_REF] Rubin | A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: The SIR algorithm[END_REF]. In the case of a Weibull model competing with very heavily censored data, that is not sufficient.

In this work, we first of all propose to improve the BRM method by amending the proposal distribution to make it more effective in terms of variance. In The article is organized as follows. The Weibull competing risks model with masked failure causes is introduced in Section 2. Section 3 deals with the handling of missing data. The use of the EM and S-EM algorithms is reviewed.

Section 4 is devoted to Bayesian estimation. The Gibbs sampler and our further improvements to the BRM method are presented. In Section 5, the different procedures are compared through simulations and a reliability data set. A discussion section and a conclusion end this paper.

The Weibull competing risks model

Competing risks (CR) models encompass any failure process in which there is more than one distinct causes of failure. In a reliability framework the objective is then to estimate the failure rate of the main causes of the system failure, when the exact cause of system failure is unknown.

In this article, we assume that there are two dominant failure mechanisms in competition. This assumption is often realistic in practice and allows us to simplify the presentation. The following assumptions are made throughout this paper. Let us consider n identical systems put on a life test and let T 1 , . . . , T n be the lifetime of these systems. Assume that

T i = min k=1,2 T ik (1)
where (T ik , k = 1, 2) are independent latent failure times corresponding to two competing causes of failure; T ik has the Weibull probability density function (pdf)

f (t | β k , η k ) = β k η k t η k β k -1 e -t η k β k , (2)
where β k > 0 and η k > 0 are respectively the shape and scale parameters of the k-th Weibull distribution. Its corresponding reliability function R W and failure rate λ W are for t ≥ 0 The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM-algorithm[END_REF] or data augmentation methods [START_REF] Tanner | From EM to data augmentation: the emergence of MCMC Bayesian computation in the 1980s[END_REF] incorporate the missing data into the likelihood to improve the bias of the maximum likelihood estimator.

R(t | β k , η k ) = e -t η k β k (3) λ(t | β k , η k ) = f (t | β k , η k ) R(t | β k , η k ) = β k η k t η k β k -1 . (4
For masked Weibull competing risks, it is possible to consider the full missing data structure with the label and failure times of each cause, z = (o i : e i = 1; t ik , k = 1, 2, i = 1, . . . , n) ; or the simple missing data structure with only the labels, z = (o i : e i = 1, i = 1, . . . , n). In the following, the the methods of estimation that we propose make use of the full or the simple missing structure.

This important point will be specified explicitly for all the presented algorithms.

Identifiability

The question of parameter identifiability is a prerequisite for estimation. A parametric model is identifiable if a model cannot be obtained by two different

parameterizations: if f (t | θ) = f (t | θ ′ ) then θ = θ ′
with the exception of one permutation. There is a lot of works on identifiability conditions for the competing risks model. For the model under study, [START_REF] Liang | Identifiability of masking probabilities in competing risks models with emphasis on Weibull models[END_REF] have proved that a Weibull masked CR is identifiable if and only if the shape parameters (β k ) are not equal. In the following we will assume that β 1 < β 2 .

Graphical analysis

The presence of two main competing causes of failure can be related to system configuration (e.g two critical components connected in series) or using sim-ple graphical analyses of failure data. In particular, Weibull quantiles-quantiles plot (QQ-plot) is a quick and simply confirmatory graphical method cf.

[3] : for standard single Weibull distribution, the Weibull QQ-plot is lined up; whereas for Weibull CR, the Weibull QQ-plot is convex. In addition, a crude point estimate of parameters can be obtained from the Weibull QQ-plot: the smallest β is 100 the predominant mode for small mission times. On the other hand, the largest β corresponds to the predominant mode for high mission times (see [3]). Therefore, with β 1 < β 2 , the tangent on the left of the QQ-plot allows to estimate (β 1 ; η 1 ) and the tangent on the right (β 2 ; η 2 ), see Figure 1. These estimates are not precise, nor robust but nonetheless can be used as starting values for 105 the various algorithms proposed hereafter.

Figure 1: Weibull QQ-plot for competing risks data: the QQ-plot is convex and for small t, R W CR ≃ R(t | β 1 ; η 1 ) and for large t, R W CR ≃ R(t | β 2 ; η 2 ). Thus, the left tangent (right resp.) allows to estimate (β 1 ; η 1 ) ((β 2 ; η 2 ) resp.) cf.

[3]

Missing data and maximum likelihood methods

Let us denote θ = (θ 1 , θ 2 ), with θ 1 = (β 1 , η 1 ), θ 2 = (β 2 , η 2 ), the parameters of the two competing Weibull distributions and y = (t, e), the observed data.

When causes of failure are masked, the likelihood function of the competing risks model is

L(θ | y) = i: ei=1 f (t i | θ 1 )R(t i | θ 2 ) + f (t i | θ 2 )R(t i | θ 1 ) • i: ei=0 R(t i | θ 1 )R(t i | θ 2 ) = i: ei=1 λ(t i | θ 1 ) + λ(t i | θ 2 ) • n i=1 R(t i | θ 1 )R(t i | θ 2 ). (5)
A direct maximization of likelihood is possible in the case of two competing causes, with a quasi Newton algorithm cf. [START_REF] Dennis | Quasi-Newton methods, motivation and theory[END_REF]. It will be proposed in the numerical experiment section and we will see that the estimator is strongly 110 biased for heavily censored data and weakly separated components which is a pathological case for maximum likelihood estimation.

The EM algorithm

The Weibull competing risks model at hand is a typical hidden structure data model for which the EM algorithm is appropriate. The EM algorithm consists of maximizing iteratively the expectation of the completed likelihood conditionally on the observed data, cf. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM-algorithm[END_REF]:

E-step: let θ (b)
being the current value of the CR model parameter, compute:

Q(θ, θ (b) ) = E f (z | y,θ (b) ) [ln L(θ | y, z)] ; (6)
M-step: then

θ (b+1) = arg max θ Q(θ, θ (b) ).
For example, considering the simple missing data structure, we have z = o and

ln L(θ | y, o) = i:ei=1 ln λ(t i | θ oi ) + n i=1 ln R(t i | θ 1 ) + ln R(t i | θ 2 ).
Thus,

Q(θ, θ (b) ) = 2 k=1 Q(θ k | θ (b) )
where

Q(θ k | θ (b) ) = i: ei=1 p(o i = k | y; θ (b) k ) ln λ(t i | θ (b) k ) + n i=1 ln R(t i | θ (b) k ) (7)
where θ

(b)
k is the current value of the Weibull parameters of the kth cause at iteration (b) and

p(o i = k | y; θ (b) k ) = λ(t i | θ (b) k ) 2 ℓ=1 λ(t i | θ (b) ℓ ) , (8)
is the posterior probability of the cause of failure. It leads to the maximization

of Q(θ k | θ (b) ) for each cause k = 1, 2 cf.
Eq. ( 7):

1

β (b+1) k + i: ei=1 p(o i = k | y; θ (b) k ) ln t i i: ei=1 p(o i = k | y; θ (b) k ) - n i=1 t β (b+1) k i ln t i n i=1 t β (b+1) k i = 0 (9)
and

η (b+1) j =       n i=1 t β (b+1) k i n i=1 p(o i = k | y; θ (b) k )       1 β (b+1) k . ( 10
)
Remarks:

The advantages and drawbacks of EM algorithm are well-documented (see

120
for instance [START_REF] Mclachlan | The EM algorithm and extensions[END_REF]). Its main advantage is that the likelihood is increasing at each iteration. The E and M steps lead generally to closed form and simple formulas easy to program: the complexity of the M step is analogous to the complexity of the ML estimation of the parameters of a single Weibull distribution. Its main drawbacks are that it can converge painfully and its solutions can be highly 125 dependent on its initial position, especially in a multivariate context.

The Stochastic EM algorithm

A stochastic version of EM (S-EM) leads to simpler equations in the M-step and it is by far less sensitive than EM to its starting value cf. [START_REF] Celeux | Stochastic versions of the EMalgorithm: an experimental study in the mixture case[END_REF]. The principle of the S-EM algorithm is to complete, at each iteration, the missing data by drawing them according to their conditional distribution knowing the observed data and the current value of the parameters.

Remarks:

The stochastic EM algorithm generates an ergodic Markov chain. Thus a sequence of parameter estimates via S-EM is expected to visit the whole parameter space with long sojourns in the neighborhood of the local maxima of the likelihood function. Therefore, S-EM algorithm offers a relevant starting values for EM algorithm, see [START_REF] Biernacki | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF].

Missing data restorations

Two versions for restoring missing data are possible according that the simple or the full is considered. In the full version, the origins of the observed failure times and all the censored failure times are simulated at each iteration (full restoration). In the simple version, only the origins of the observed failure times are simulated (simple restoration).

Full restoration

In the full restoration version, the missing data restoration is as follows: E-step: for i = 1, . . . , n compute the posterior probability of the cause of failure

p(o i = k | y, θ (b) ), k = 1, 2 cf. Eq (8).

S-step:

Simulating the causes of failure,. the i-th failure is assigned to one of the causes; (o

(b) i , i = 1, . . . , n) are simulated according to p(o i = k | y, θ (b) ), k = 1, 2 : o (b) i ∼ p(o | y, θ (b) ). ( 11
)
Simulating the failure times for the censored data,. t i with e i = 0, the failure times t ik of the two causes cf. Eq. ( 1), k = 1, 2, are simulated according to the Weilbull distribution with parameter θ

(b)
k , with the constraint of being greater than t i :

t (b) ik ∼ f (t | θ (b) k ) R(t i | θ (b) k ) . ( 12
)
Simulating the censored failure times for the observed data,.

t i with e i = 1, it is assigned to the cause o (b) i , (i.e t ik = t i if o (b) i = k); the failure time of the other competing cause, (t ik , k ̸ = o (b) i ), is simulated according to the Weibull distribution with parameter θ (b) k , k ̸ = o (b)
i , with the constraint of being greater than t i :

t (b) ik        = t i if o (b) i = k; ∼ f (t | θ (b) k ) R(t i | θ (b) k ) if o (b) i ̸ = k. (13)
M-step: it consists of computing the ML estimate of each component, (θ

(b+1) k , k = 1, 2), from the completed sub-samples ( t (b)
ik , i = 1, . . . , n; k = 1, 2), by a standard procedure, see for instance [START_REF] Murthy | Weibull Models[END_REF].

Simple restoration

In the simple restoration version, he missing data restoration is as follows:

E-step: for i = 1, . . . , n compute the distribution p(o i = k | y, θ (b) ), k = 1, 2 cf. Eq (8);
S-step: Simulating the causes of failure, the i-th failure is assigned to one of the causes; (o 11); M-step: it consists of computing the ML estimates of each component, (θ

(b) i , i = 1, . . . , n) are simulated according to p(o i = k | y, θ (b) ), k = 1, 2 , cf. Eq. (
(b+1) k , k = 1, 2)
, from the censored sub-samples (t i such that e i = 1 and o (b) i = k ; t i such that e i = 0) by a standard procedure, see for instance [START_REF] Murthy | Weibull Models[END_REF].

Remark:

When restoring labels in S-step (full or simple), it can happen that no failure is attributed to a cause. Then, the S-step is rerun until each cause has a minimum number of failures.

Bayesian inference

The Bayesian estimation is an alternative to the ML estimation which is relevant when the information provided by the observed data is poor. It is worth recalling that we have four parameters to estimate θ = (β 1 , η 1 , β 2 , η 2 ) and the censoring rate is about 70%. The Bayesian estimation consists of incorporating prior information to reduce the uncertainty of the parameters. With the aid of a prior probability distribution π(θ), statistical inference is made from the posterior distribution of the parameters knowing the data. This posterior distribution, π(θ | y), is obtained using Bayes formula,

π(θ | y) = π(θ)L(θ | y) π(θ)L(θ | y)dθ ∝ π(θ)L(θ | y). ( 14
)
This distribution combines the prior distribution for the parameters π(θ) with the likelihood which contains information on the parameter provided by the data. Thus, Bayesian inference is expected to be useful to produce "regularized" maximum likelihood estimation. It is used in the present article from this point of view.

Prior information may be derived from data of previous studies or expert advices. It quantifies uncertainty about the parameters and it is formalized in the prior distribution π(θ). Therefore, we consider Bayesian inference in order to derive a pointwise estimate of θ which overcomes the difficulties of the EM algorithm for maximizing the likelihood in the challenging context of heavily censored masked Weibull competing risks model. A pointwise estimate of the parameter can be derived from this posterior distribution. Such Bayesian estimators are proper to produce regularized estimation of the model parameters protected from spurious and insensible local maxima of the EM algorithm.

Prior distribution

The choice of priors distribution is the key problem in the Bayesian inference.

The prior distribution expresses both the engineer's expertise on parameter values and statistical uncertainty. The same joint prior has been chosen for the parameters of the two competing Weibull distributions.

Shape parameter

Thanks to its shape parameter, β, Weibull distributions have the ability to assume the characteristics of many different types of lifetimes. The range of β in various applications are addressed in [START_REF] Barringer | Weibull databases[END_REF] and [START_REF] Jiang | A study of weibull shape parameter: Properties and significance[END_REF] for instance. The same beta distribution, symmetric on range interval [b min ; b max ], has been chosen as prior for shape parameters (β k ):

π(β)= Γ(p + q) Γ(p)Γ(q) (β -b min ) p-1 (b max -β) q-1 (b max -b min ) p+q+1 • 1 [0;1] β -b min b max -b min , (15)
where

Γ(t) = +∞ 0 t t e -t dt. Thus β -b min b max -b min ∼ B(p, q). (16)
The hyperparameters (b min , b max ) reflect engineer expertise on the parameters (β k ). They shall be fixed in accordance to the field failure data. The hyperparameters b min and b max were chosen according to [START_REF] Barringer | Weibull databases[END_REF]: b min = 0.5 and b max = 10. Hyperparameters (p, q) reflect uncertainty about (β k ). A vague informative prior was chosen by setting p = q = 1.1.

Scale parameter

For scale parameters (η k ), the prior density of η k given β k is a Generalized Inverse Gamma distribution (GIG):

π(η k | β k ) = β k a β k b k k Γ(b k ) η -(β k b k +1) k e - a k η k β k , (17)
where (a k , b k ) are hyperparameters. Thus η k β k has an inverse gamma distribution.

Gibbs sampling

A standard solution to get an approximation of the full posterior distribution through Monte Carlo Markov Chains (MCMC) algorithms is the Gibbs sampling, see for instance [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF] or [START_REF] Casella | Explaining the Gibbs sampler[END_REF]. The Gibbs algorithm consists of simulating iteratively all unknown quantities (parameters and missing data) according to their full conditional distributions cf. [START_REF] Roberts | Simple conditions for the convergence of the Gibbs sampler and Metropolis-Hastings algorithms[END_REF]:

For b = 0, . . . B -1: o (b+1) i ∼ p(o | y, θ (b) ), for i = 1, . . . , n; β (b+1) k ∼ π(β | η (b) k , {t i : o (b+1) i
= k}, e), for k = 1, 2;

θ (b+1) k ∼ π(θ | β (b+1) k , {t i : o (b+1) i = k}, e), for k = 1, 2.
Under mild conditions, and after a large number of iterations, ( θ (1) , . . . , θ (B) ) is a Markov chain of which distribution converges to π(θ | y), the posterior distribution. Therefore, after a burning period of M iterations the Gibbs estimator is computed from the B -M last terms:

θ(GIBBS) = arg max b∈{M +1,...,B} log π( θ (b) | y). (18)
Remarks:. Note that most of the conditional distributions are quite simple, except for β k : however, our choice of the prior ensures that the distribution is log-concave and thus the adaptatve rejection method of [START_REF] Gilks | Adaptive rejection sampling for Gibbs sampling[END_REF] can be used to simulate the full conditional distribution of β k , k = 1, 2. The Gibbs sampler is described more precisely in the appendix A.

The Gibbs sampler produces a Markov chain ( θ (1) , . . . , θ (B) ), whose stationary law is the posterior distribution. However, the convergence is not assured and it may exhibit poor mixing properties in the case of heavily censored data cf. [START_REF] Casella | Explaining the Gibbs sampler[END_REF].

BRM method and improvements

Actually, in the present context, computing a good approximation of the posterior distribution appears to be quite difficult, but estimating the posterior mean or mode is an easier task. We then proposed a non-iterative method that is effective in the case of very heavily censored data, see [START_REF] Bacha | Estimation de durées de vie fortement censurées[END_REF] for the Weibull CR model and [START_REF] Ducros | Bayesian estimation of Weibull mixture in heavily censored data setting[END_REF] for a mixture of Weibull distributions. We sketch the likelihood maximization algorithm with missing data completed by a Bayesian restoration (BR-LM algorithm):

Algorithm BR-LM
For a large number of runs 1 ≤ b ≤ B

Step 1 Bayesian sampling: the parameter θ Step 3 Maximization: the completed log-likelihood is maximised :

θ (b) = arg max θ ln L(θ | y, z (b) );
Step 4 Estimation: the mean of the posterior distribution is then obtained using a importance sampling technique cf. [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]: for a proposal distribution ρ

µ = E π(θ|y) [θ] = E ρ(θ) π(θ|y) ρ(θ) θ .
Here ρ, an estimate of the distribution of ( θ (b) , b = 1, . . . , B) is used as proposal distribution:

ρ(θ) = 1 B • h θ B b=1 K θ -θ (b) h θ , ( 19
)
where K is a Gaussian kernel. The window size, h θ , minimizing the mean integrated square error is given in [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF]. Then assign a weight to each θ (b) :

w b = π(θ (b) )L(θ (b) | y) ρ(θ (b)
) .

Defining the normalized weights

w b = w b / k w k .
then importance sampling estimate of µ is

µ = 1 B b w b • θ (b)
.

Remarks:

One of the traps of importance sampling is to use a proposal distribution too far from target distribution (i.e. θπ(θ|y) ̸ = µ ρ(θ) with ρ(θ) > 0 ). Many weights w b are then almost zero, with some are very high. Rubin suggests a re-sampling cf. [START_REF] Rubin | A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: The SIR algorithm[END_REF]. In the case of a Weibull model competing with very heavily censored data, that is not sufficient.

Here we propose a first improvement of the BR-LM algorithm. In order to correct the weight problem, we propose to use the EM algorithm on every θ (b) .

The following EM step is added to the BR-LM algorithm

Step 3' EM: The EM algorithm is initialized with θ (b) .

The Step 3', allows to avoid having many almost zero weights and thus to improve the variability of the estimator.

Remarks:

To be as pragmatic as possible, we have chosen a very uninformative prior distribution and consequently the estimation of the expectation of likelihood completed with the Bayesian sampler can be very bad. The use of the EM algorithm then improves the estimation of the expectation of completed likelihood.

In step 2 of the BR-LM algorithm, the completed log-likelihood is maximized.

One can also directly maximize the log posterior distribution completed with missing data :

Algorithm BR-PM

Step 1 Bayes: the parameter θ Step 3 Maximization: the completed log-posterior distribution is maximised :

θ (b) = arg max θ E f (z | y,θ (m) ) [ln π(θ | y, z)] (20)
Details of the maximization can be found in the Appendix.

Step 3' EM-MAP: The EM algorithm is easily modified to produce the maximum a posteriori (MAP) cf. [START_REF] Mclachlan | The EM algorithm and extensions[END_REF], see the appendix B. Let us denote θ

(b)
M AP the result of M iterations of the algorithm EM-MAP with θ (b) as initialization.

Step 4 Estimation: the mean of the posterior distribution is then estimated using a importance sampling technique with the empirical distribution of ( θ

(b) M AP , b = 1, • • • , B) as proposal distribution.
Remarks: The direct maximization of the expectation of the completed logposterior distribution requires the use of a numerical method. However, we note that this maximization appears as a regularization of the maximization of the completed log-likelihood of BR-LM, with the prior distribution as a penalization function. The EM-MAP algorithm, like the EM algorithm for BR-LM, is used to improve the law proposition used in preferential sampling. As will be seen in the next section, this improves the estimation.

Numerical experiments

Numerical experiments are carried out first on simulated data and then on a set of reliability data.

Simulated data sets

We simulated the operation of a serial system of two independent components. The component lifetimes are distributed according to Weibull's law: T 1 ∼ W(β 1 = 1.5 , η 1 = 2 500) green: The sample size was n = 200 and the censorship rates 70% and 20% have been used. It can be seen in Figure 2 that the number of observed lifetimes is not evenly distributed between the two components for the simulated data set. For heavily censored data (e.g. 70% censorship), there is very few data for component 1 and it is extremely difficult to estimate its distribution. As a matter of fact, there is few published works with such high censoring rates, although they are common in industrial contexts.

t 1 < t 2 < censoring T 2 ∼ W(β 2 = 4 , η 2 = 1 000) red : t 2 < t 1 < censoring white: t 1 , t 2 > censoring
The performance of the estimators is evaluated in terms of bias and root mean square error (RMSE) from 500 replications of the dataset :

Bias rel ( θ) = Bias( θ) θ RM SE rel ( θ) = RM SE( θ) θ .
From these 500 replications, we compare the performance of the estimators provided by (i) the EM algorithm initiated with the S-EM algorithm with B = 100 iterations, (ii) the Gibbs sampling estimates with 5 000 iterations and 4 700 discarded "burn-in" iterations, (iii) the algorithm BR-LM and (iv) the algorithm BR-PM with B = 50 000 runs.

For comparison, a high (70%) and low (20%) level of censorship were tested.

The results are shown in Table 1, for a 70% censorship level and in Table refTab20, for a 20% censorship level. The shaded cells correspond to the optimal estimate either in terms of relative bias or RMSE. It can be noted the Bayesian restoration algorithms outperform EM with S-EM algorithms and the Gibbs sampling, even with moderately censored data. Indeed, even when there is no censorship on the right in the data collection process, there is always missing data: on the one hand the cause of the failure is masked, on the other hand the other competing cause is censored by the cause involved in the failure.

The algorithms BR-LM and BR-PM provide estimations with similar biases; however, more often the algorithm BR-PM provides estimates with a lower RSME.

As mentioned previously (see Figure 2), the high rate of censorship makes (β 1 , η 1 ) estimation much more difficult than that of (β 2 , η 2 ). It is important to note this because W(β 1 , η 1 ) is the distribution of the most reliable component.

Real data

Data from aircraft windshields maintenance are analyzed (cf. [START_REF] Murthy | Weibull Models[END_REF]). The windshield of an aircraft is a complex equipment, essentially composed of several layers of materials, including a very resistant outer skin with a heating coating just below it. All the layers are laminated at high temperature and under high pressure. Windshield failures are usually due to damage to the outer layer or failure of the heating system. During a maintenance operation, 153 aircraft windshields were inspected: 88 were classified as faulty and the remaining 66 as still operational. Their operating times (in thousands of hours) since they were commissioned were recorded. (cf. Table 3). The Weibull QQ-plot suggests two main competing causes of failures (see Figure 3). The poor performance of the Gibbs algorithm illustrates the fact that it converges very slowly (cf. [START_REF] Casella | Explaining the Gibbs sampler[END_REF]). By comparing the intervals we have noticed that the algorithm BR-PM provides more accurate intervals for the scale parameter of the first mode.

Discussion

Algorithmic issues

For iterative algorithms, especially for S-EM and Gibbs sampler, the labels of components may be randomly switched. This label switching problem that complicates the estimation must then be addressed. The solution chosen is the one proposed by [START_REF] Celeux | Bayesian inference for mixture: The label switching problem[END_REF]: algorithm is post-processing by a k-means clustering.

If no failure or less than five are assigned to a cause, then the algorithm is shut down. This has an impact on the RMSE.

Full versus simple restoration

Table 7 provides comparisons between Full and Simple restoration for the BR-PM algorithm (cf. section 3.3). It can be seen that the Simple restoration version gives poorer results in terms of bias. is Bayesian sampling. The completed likelihood maximization step then acts as an importance re-sampling. The EM step then acts as a regularization step.

Component 1 Component 2

Step 1 : Replications of S-EM estimators.

Step 2 : EM algorithm with starting value from step 1 Step 3' : EM algorithm with starting value from step 2

k ) = (p -1) ln(β k -b min ) + (q -1) ln(b max -β k ) + ln(β k ) + β k b k ln(a k ) -β k b k ln(η (b) k ) - a k η (b) k β k + #{o (b) i = k} ln(β k ) + β k i: oi=k ln t i η (b) k - n i t i η (b) k β k = (p -1) ln(β k -b min ) + (q -1) ln(b max -β k ) + (#{o (b) i = k} + 1) ln(β k ) + β k b k ln a k η (b) k + i: oi=k ln t i η (b) k - a k η (b) k β k - n i t i η (b) k β k (22) ∂ ln f ∂β k (β k | t, o (b) , η (b) k ) = p -1 β k -b min - q -1 b max -β k + 1 + #{o (b) i = k} β k + b k ln a k η k + i:oi=k ln t i η k -ln a k η k a k η k β k - n i ln t i η k t i η k β k (23) ∂ 2 ln f ∂β 2 k (β k | t, o (b) , η (b) k ) = - p -1 (β k -b min ) 2 - q -1 (b max -β k ) 2 - #{o (b) i = k} + 1 β 2 k -ln a k η (b) k 2 a k η (b) k β k - i ln t i η (b) k 2 t i η (b) k β k (b)
Hence, the adaptative rejection method proposed by [START_REF] Gilks | Adaptive rejection sampling for Gibbs sampling[END_REF] to sample univariate log-concave probability density function can be used here. For a simple restoration we get :

∂E f (o | y,θ (m) ) [ln L(θ | y, o)] ∂η k (η k ) + ∂ ln π ∂η k (η k | β k ) = i:ei=1 p(o i = k | y; θ (m) ) ∂ ln λ ∂η k (t i | θ k ) + n i=1 ∂ ln R ∂η k (t i | θ k ) + ∂ ln π ∂η k (η k | β k ) =- β k η k i:ei=1 p(o i = k | y; θ (m) ) + β k η β k +1 k n i=1 t β k i - b k • β k + 1 η k + β k η β k +1 k a β k k = β k η k n i=1 t β k i + a β k k η β k k - i:ei=1 p(o i = k | y; θ (m) ) -b k - 1 β k (24)
Thus

∂E f (o | y,θ (m) ) [ln L(θ | y, o)] ∂η k (η k ) + ∂ ln π ∂η k (η k | β k ) = 0 ⇔ n i=1 t β k i + a β k k = η β k k i:ei=1 p(o i = k | y; θ (m) ) + b k + 1 β k .
Therefore [ln(t i ) -ln(η k )] t

η k =       n i=1 t β k i + a β k k i:ei=1 p(o i = k | y; θ (m) ) + b k + 1 β k       1/β k . ( 25
β k i + p -1 β k -b min - q -1 b max -β k + 1 β k + b k ln(a k ) -b k ln(η k ) + 1 η β k k [ln(η k ) -ln(a k )] a β k k = i:ei=1 p(o i = k | y; θ (m) ) + 1 β k + p -1 β k -b min - q -1 b max -β k + ln(η k ) 1 η β k k n i=1 t β k i + a β k k -b k - i:ei=1 p(o i = k | y; θ (m) ) - 1 η β k
β k + p -1 β k -b min - q -1 b max -β k - i:ei=1 p(o i = k | y; θ (m) ) + b k + 1 β k       n i=1 ln(t i ) t β k i + ln(a k ) a β k n i=1 t β k i + a β k k       + 1 β 2 k ln       n i=1 t β k i + a β k

  e i = 0 if t i is right censored. ; e = (e 1 , . . . , e n ) y i observed data: y i = (t i , e i ); y = (y 1 , . . . , y n ) o i cause of the i-th failure; o = (o 1 , . . . , o n ) z i missing data : simple missing ; the cause of failure z i = o i full missing ; the cause of failure and the lifetime of each cause z i = (o i ; (t ik )) ; z = (z 1 , . . . , z n ) z i restored missing data f, R and λ probability density function, reliability and failure rate of a Weibull distribution θ k = (β k , η k ) shape and scale Weibull parameters of the k-th cause of failure θ set of the CR model parameters = (θ 1 , θ 2 ) = (β 1 , η 1 , β 2 , η 2 ) CR model likelihood with the observed data L(θ|y, z) CR model likelihood with the complete data π(θ) prior distribution π(θ|y) posterior distribution knowing the observed data π(θ|y, z) posterior distribution knowing the observed data y completed by the missing data z E f [X] expectation for X ∼ f .

  addition a new proposal distribution is obtained from the maximization of the mean of the posterior distribution completed by the missing data. The efficiency of the proposed methods was evaluated by a large number of simulations for different levels of censoring rate. Simulations have shown that the new proposed methods are effective both in terms of relative bias and relative root mean square error.

  : the full missing data z(b) are restored according to their conditional distribution knowing the observed data and the simulated θ

  : the full missing data z(b) are restored according to their conditional distribution knowing the observed data and the simulated θ

W 1 ∼

 1 W(β 1 = 1.5 ; η 1 = 2 500) and W 2 ∼ W(β 2 = 5 ; η 2 = 1 000). The two Weibull distributions are weakly separated: this is a pathological case for maximum likelihood estimation. a) Components pdf b) synthetic data : (t 1 , t 2 )

Figure 2 :

 2 Figure 2: Distributions and synthetic data with n = 200, and 70% right censoring. The twoWeibull distributions are weakly separated: this is a pathological case for maximum likelihood estimation. Moreover with a very high censorship rate, there is very few data for component 1 and it is extremely difficult to estimate its distribution.

Tables 4 -

 4 6 provide 95% bootstrap or credible intervals for the EM with S-EM algorithms, the Gibbs sampler and the BR-PM algorithm we propose. Simulations have shown us that EM with S-EM provide estimates with high variance.

Figure 3 :

 3 Figure 3: Weibull QQ-plot for aircraft windshield cf. [11]; there are two main competing causes of failures.

Figure 4 :Step 1 :

 41 Figure 4: EM algorithm with S-EM estimate as starting values.

Figure 5 :

 5 Figure 5: BR-LM with EM algorithm

  3. Scaleπ(η k | y, o (b) , β (b) , η

.

  Thus, π(η k | y, o (b) , β (b) , η (b) (-k) ) is the generalized inverse gamma GIG #{o (b) i = k} + b k ; Maximization step of EM for the posterior distribution In Algorithm 4, the maximization step is performed by equating to zero the partial derivatives of conditionnel expectation of the completed log posterior with respect to each parameter cf. Eq. (20): E f (z | y,θ (m) ) [ln L(θ | y, z)] + ln π(θ).

  i )t β k i + ln(a k )a β k k + i:ei=1 p(o i = k | y; θ (m) ) ln(t i ) + b k ln(a k )(26)Upon substituting Eq. (25) in Eq. (26), it follows∂E f (o | y,θ (m) ) [ln L(θ | y, o)] ∂β k (β k ) + ∂ ln π ∂β k (β k ) + ∂ ln π ∂β k (η k | β k ) = 0 ⇔ i:ei=1 p(o i = k | y; θ (m) ) + 1

  p(o i = k | y; θ (m) ) + b k + i = k | y; θ (m) ) ln(t i ) + b k ln(a k ) = 0(27)References

  , . . . , o n ), with o i ∈ {1, 2} and the failure times of the k-th cause t ik , k ̸ = o i cf. Eq. (1).

	2.1. The missing data structure
	Masked competing risks model is a missing data structure model. Conse-
	quently, statistical inference is to be made on incomplete data and maximum
	likelihood could resort to EM-like algorithms.
	The observed data are as follows: (t, e) = ((t 1 , e 1 ), . . . , (t n , e n ))) where t i is
	a failure time due to cause o i ∈ {1, 2} if e i = 1 or a right censored time if e i = 0.
	The cause of failure o i ∈ {1, 2} is not observed.
	The missing data are of different natures. They include the unknown labels
	o = (o 1
	)

Table 1 :

 1 70% censoring rate: the performance of the estimators is evaluated in terms of bias and root mean square error (RMSE) from 500 replications of the data set with the parameters below.

	True	β 1	β 2	η 1	η 2	MTTF1	MTTF2
	parameters	1.5	5	2 500	1 000 2 256.86	906.40
	ML using BFGS quasi-newton method with WQQPlot starting value
	Average	2.74	3.72	1 857.41 966.00	1 698.53	872.35
	Bias rel	82.38% -25.63% -25.70% -3.40% -24.74%	-4.99%
	RMSE rel 145.58% 23.11% 49.82% 42.49%	54.10%	11.01%
			EM with S-EM starting value
	Average	2.25	6.17	1 822.76 1 189.16 1 649.58	1 119.99
	Bias rel	50.02% 23.39% -27.09% 18.92% -26.91%	23.56%
	RMSE rel	70.06% 43.03% 55.89% 112.82% 59.26%	131.06%
				Gibbs sampling
	Average	2.02	5.45	1561.07 895.39	1383.42	826.16
	Bias rel	34.99% 36.21% -37.56% -10.46% -38.70%	-8.85%
	RMSE rel	4.82%	3.69%	7.61%	0.67%	7.47%	0.72%
		BR-LM estimation with full restoration
	Average	1.86	5.99	2 024.01 1 119.38 1 829.01	1 035.53
	Bias rel	23.96% 19.85% -19.04% 11.94% -18.96%	14.25%
	RMSE rel	30.35% 28.85% 33.98% 35.57%	35.33%	34.12%
		BR-PM estimation with full restoration
	Average	1.88	4.42	2 220.02 1 073.46 1 864.74	978.45
	Bias rel	25.33% -19.26% -19.20% 7.35%	-17.37%	7.97%
	RMSE rel	25.19% 19.97% 9.64%	25.54%	9.44%	25.68%

Table 2 :

 2 20% censoring rate: the performance of the estimators is evaluated in terms of bias and root mean square error (RMSE) from 500 replications of the data set with the parameters below.

	True	β 1	β 2	η 1	η 2	MTTF1	MTTF2
	parameters	1.5	5	2 500	1 000 2 256.86	906.40
	ML using BFGS quasi-newton method with WQQPlot starting value
	Average	2.49	4.47	2 003.10 942.37 1 804.42	859.21
	Bias rel	66.06% -10.65% -19.88% -5.76% -20.05%	-5.22%
	RMSE rel 75.10% 17.83% 72.15% 3.92%	74.28%	4.30%
			EM with S-EM starting value
	Average	1.70	5.56	2 315.01 1 010.64 2 103.20	933.06
	Bias rel	13.36% 11.16% -7.40%	1.06%	-6.81%	2.94%
	RMSE rel 22.90% 15.43% 42.61% 2.98%	45.17%	3.32%
				Gibbs sampling
	Average	1.67	5.38	2 277.25 1 041.92 2 102.62	868.46
	Bias rel	11.30% 7.64%	-8.91%	4.19%	-6.83%	-4.19%
	RMSE rel 32.38% 12.78% 24.67% 3.58%	24.38%	4.64%
		BR-LM estimation with full restoration
	Average	1,58	5.14	2391.82 1 018.00 2 111.16	915.23
	Bias rel	5,28%	2.80%	-4.33%	1.20%	-6.45%	0.97%
	RMSE rel 15.47% 12.57% 13.89% 2.29%	16.85%	2.60%
		BR-PM estimation with full restoration
	Average	1.62	4.82	2 371.24 1 008.48 2 094.71	917.40
	Bias rel	8.08% -3.60% -6.35%	0.85%	-7.18%	1.21%
	RMSE rel 13.28% 14.30% 11.74% 2.64%	12.20%	3.46%

Table 3 :

 3 Aircraft windshields: lifetimes (×1 000 h). "-" indicated that the windshied had not failed, the lifetime is then censored. The censorship rate is 57.5%

	0.04	-0.046 -0.14	-0.15 -0.248 -0.28	0.301
	0.309 -0.313 -0.389 -0.487 0.557 -0.622	-0.9
	0.943 -0.952 -0.996 -1.003 -1.01	1.07	-1.085
	-1.092 1.124 -1.152 -1.183 -1.244 1.248 -1.249
	-1.262 1.281	1.281	1.303	-1.36	1.432 -1.436
	1.48	-1.492 1.505	1.506	1.568	-1.58	1.615
	1.619	1.652	1.652 -1.719 1.757 -1.794 1.795
	1.866	1.876	1.899	1.911	1.912	1.914 -1.915
	-1.92 -1.963 -1.978 1.981	2.01	2.038 -2.053
	-2.065 2.085	2.089	2.097 -2.117 2.135 -2.137
	-2.141 2.154 -2.163 -2.183	2.19	2.194	2.223
	2.224	2.229	-2.24	2.3	2.324 -2.341 2.349
	2.385 -2.435 -2.464 2.481 -2.543 -2.56 -2.592
	-2.6	2.61	2.625	2.632	2.646	2.661	-2.67
	2.688 -2.717 -2.819 -2.82	2.823 -2.878	2.89
	2.902	2.934	-2.95	2.962	2.964	3.0	-3.003
	-3.102 3.103	3.114	3.117	3.166 -3.304 3.344
	3.376	3.385	3.443	3.467	3.478 -3.483	-3.5
	3.578	3.595 -3.622 -3.665 -3.695 3.699	3.779
	3.924 -4.015 4.035	4.121	4.167	4.24	4.255
	4.278	4.305	4.376	4.449	4.485	4.57	4.602
	-4.628 4.663	4.694 -4.806 -4.881 -5.14

Table 4 :

 4 Aircraft windshields: n = 153 and the censorship rate is 57.5%. Parameters 95% bootstrap intervals with the Algorithm S-EM and EM

	θ	β 1	η 1	β 2	η 2
	25%	0.636 393 645.2 2.834 3 526.2
	97.5% 0.643 583 570.7 2.838 3 527.8
	Table 5: Aircraft windshields: n = 153 and the censorship rate is 57.5%. Parameters 95%
	CI's with the Gibbs sampling
	θ	β 1	η 1	β 2	η 2
	25%	0.554	511 .0	2.660	3 429.4
	97.5% 0.825 395 756.0 3.361 5 3 805.4

Table 6 :

 6 Aircraft windshields: n = 153 and the censorship rate is 57.5%. Parameters 95%

	credible intervals with the Algorithm BR-PM
	θ	β 1	η 1	β 2	η 2
	25%	0.642 5 386 701.2 2.838 3 527.7
	97.5%	0.670	394 590.4 2.851 3 534.1

Table 7 :

 7 Comparison Full versus Simple restoration for BR-PM algorithm.

	Synthetic data

)

  ∂E f (o | y,θ (m) ) [ln L(θ | y, o)]

				∂β k		(β k ) +	∂ ln π ∂β k	(β k ) +	∂ ln π ∂β k	(η k | β k )
	=	i:ei=1	p(o i = k | y; θ (m) )	∂ ln λ ∂β k	(t i | θ k )
		+	n i=1	∂ ln R ∂β k	(t i | θ k ) +	∂ ln π(β k ) ∂β k	+	∂ ln π ∂β k	(η k | β k )
	=	i:ei=1	p(o i = k | y; θ (m) )	1 β k	+ ln(t i ) -ln(η k )
		-	1 η β k k	n i=1

Bayesian restoration

The EM algorithm is rather stable but tends to be biased in the case of heavily censored data. As can be seen in Figure 4, the initialization with the S-EM algorithm does not correct the bias. With the help of Figure 5


Chloé Mercier

Frédéric Alexandre

Thierry Viéville

Reinforcement Symbolic Learning

Keywords: Reinforcement Symbolic Learning, Ontology Edit Distances, Models for Learning Sciences

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Understanding how humans solve problems and learn is a key issue in education, and one of the transversal competencies sometimes referred to as "21st-century skills". At the cognitive level, we need to consider both exploration and exploitation strategies directed either by a stimulus-driven behavior or toward different concurrent goals in a goal-driven behavior [START_REF] Alexandre | A global framework for a systemic view of brain modeling[END_REF]. At a computational level, the typical paradigm is reinforcement learning (RL) with a final reward (success of the task) and some intermediate rewards (discoveries of affordances, partial result regarding the goal), more precisely intrinsically motivated reinforcement learning [START_REF] Singh | Intrinsically Motivated Reinforcement Learning: An Evolutionary Perspective[END_REF]. The latter family of models is extensively used in cognitive neurosciences to model high-level executive control functions (e.g., [START_REF] Domenech | Executive control and decision-making in the prefrontal cortex[END_REF][START_REF] Rmus | The role of executive function in shaping reinforcement learning[END_REF]).

In a recent development, we have introduced the operationalization of a creative problem-solving task, via the construction of an ontology 1 [START_REF] Mercier | Formalizing Problem-Solving in Computational Thinking : an Ontology approach[END_REF]. The state is defined by the configuration of the objects manipulated during the task and some other observables regarding the learner 2 . The key point is that, in our case, the internal state of the subject and the external state of the task material constitute a complex structure, modeled here by a set of "statements"; that is to say, in the ontological vocabulary, entities typed by classes and linked together Supported by Inria, AEx AIDE https://team.inria.fr/mnemosyne/en/aide. 1 Concretely the ontology role is to help specifying the representation of the world as conceived by the learner, i.e., what is observed and what is to be inferred during the learning activity, to solve the task. This also helps to verify the specification coherence of the model, which allows us to infer some assumptions about non observable elements of the learning process, as detailed in the supplementary material accessible here https://gitlab.inria.fr/line/aide-group/creacog. 2 See here for a video illustration.

through relationships labeled by properties. Unlike usual mechanisms based on Markov chains, the state space to consider in our RL setup is thus not reduced to an unordered finite enumeration. We would like to study here to what extent reinforcement learning could be designed on such state spaces.

In "symbolic reinforcement learning" (e.g. [START_REF] Garnelo | Towards Deep Symbolic Reinforcement Learning[END_REF]), deep neural networks transform raw perceptual data into a symbolic representation which is then fed to a symbolic module that might perform, for instance, action selection. Other approaches such as [START_REF] Ma | Interpretable Reinforcement Learning With Neural Symbolic Logic[END_REF] propose architectures where a numerical reinforcement algorithm communicates with a reasoner. Here we would like to explore another track and make the reinforcement algorithm work directly on the symbolic data space itself in a more integrative way, which we propose to call "reinforcement symbolic learning".

Symbolic state space specification

General framework: An agent interacts with its environment. At a given discrete time, it perceives a part of the environment, i.e., a stimulus, including a reward. It infers elements (e.g., causes) from this input cue, including the computation of the next action. In our case, we consider a potentially hypermnesic agent for which any previous input, including rewards, might be part of the internal state. This choice is directly linked to the notion of episodic memory, an episode here being represented by an event sequence (a list of times of occurrence of atomic events). This frees us from the "Markovian" constraint taking into account not just one step in the past, at the cost of a multiplicative increase of the state space. This is going to be manageable thanks to the hierarchical structure of our state space, which encompasses a lot of information without the need for an exhaustive enumeration. Furthermore, the state space no more reduces to a flat enumeration of state values, while the number of possible ontology construction of size S on a vocabulary if size V is an order of magnitude higher (i.e. roughly in O(V3 S)).

Input structure: The input is a hierarchical data structure time sequence, with information regarding the task and the learner. Syntactically this corresponds to tuples of named values {• • • , name : value, • • • }; each value has a "type" (determined by a predefined schema 3 making explicit the value set), and may be either another tuple or a literal or the meta-value undefined (i.e., expected by the schema, but not present in the given context). Literal values are taken among a finite enumeration of qualitative values (e.g., a color set) or quantitative values (i.e., finite precision bounded values).

Comparing two inputs: An input s is thus a tree data structure, equipped 4 with a a partial semi-order compatible with an extended semi-distance. This means that two values may be equal, indistinguishable (i.e., too close to be ordered, thus equal or not), comparable or incomparable (i.e., too different to be compared). This mechanism not only allows to define a distance between two inputs (as the minimal cost of editing sequences transforming one input into another), but also to make explicit which node has been added, deleted, or changed. Thus, it offers a "geodesic", i.e., a path of transformation from one structure to another, allowing us to interpolate intermediate input structure between both of them.

Inferring other elements from input: Each data structure is translated in terms of RDF statements as follows. Each tuple is a "subject" and each named value corresponds to a relationship labeled by a "property", the value being the "object" targeted by the relationship. This transformation allows us to generate an ontology 5 , offering the possibility to perform inferences, thus implementing the learner behavior at a pure symbolic level. Conversely, each RDF ontology graph may be mapped back onto the data structure, with tuples having the value undefined if the corresponding statement is absent from the ontology after reasoning, and defined as the property object•s otherwise -hence the need for a pre-defined schema.

3 Reinforcement learning on symbolic state space.

Let us consider a concrete example on the most common algorithm setup, i.e., Qlearning with -greedy exploration, namely an action-state value function Q : S × A → R. At each step this function is updated, using the weighted average of the old value and the new information, while the action is chosen to either maximize the reward value, given the state, or with a small probability randomly explore new actions. This algorithm is known as "model-free", however, we are using it in a non-conventional way with inferences we generate from prior information on symbolic states, bringing it somewhat closer to model-based algorithms.

The key point is that, given the largeness of the state space, each state value is very likely different from another, so that one state value is very likely visited once, making it impossible to use the usual update rule on tabulated values Q[s t , a t ]. However, given a new state value s t and reward r t+1 , we can easily update all preceding Q-values in a neighborhood.

Considering an exponential weighting of radius ρ, for a learning factor α, a discount factor γ, this writes, for all Q-values:

Q[s, a t ]+= α e -d(s,st)/ρ (r t+1 + γ max a Q(s t+1 , a) -Q[s, a t ]) where d(s, s t ) stands for the predefined edit distance between both states. During an epoch of T steps, it means that we have to compute O(T 2 ) edit distances, each of these operations having a user-defined positive cost, related to the literal extended semi-distances. The key point is that we consider restricted edit distances preserving the tree filiation, computable in polynomial time [START_REF] Ouangraoua | A Constrained Edit Distance Algorithm Between Semiordered Trees[END_REF], which would not have been the case otherwise, or if considering the tree as a general graph or ontology portion. 5 In the sense of RDF and OWL2. and Q-value updates, but this complexity depends only on the trajectory length, rather than on the state space itself.

The computation of the maximum reward prediction max a Q(s t+1 , a), requires to both (i) interpolate the Q-value given available tabulated values and (ii) enumerate action candidates, including unprecedented actions. For (ii) consider a set of potential actions a k including previous actions, predefined prototypical actions, and putative actions generated by an external process. For (i) we can use an exponential interpolation Q(s, a k ) = st,at e -(d(s,st)+d(a k ,at))/ρ Q[s t , a t ] st,at e -(d(s,st)+d(a k +,at))/ρ in coherence with the previous design choice.

Conclusion

We have here all ingredients to apply well-established Q-learning mechanisms not only on an enumeration of indexed states, but on a rich semantic structure, which was the goal of this preliminary work. Meanwhile, in a companion study, we explore ways to map such a semantic structure onto a neuronal vector space [START_REF] Mercier | Ontology as neuronalspace manifold: towards symbolic and numerical artificial embedding[END_REF]. We are making a prospective presentation here to share the scientific idea, it goes without saying that this is only an open issue to be developed.

In the sense of https://json-schema.org.

We consider edit operations given an input (l+) adding, (l-) deleting or (l#) changing a value in a list, (t+) defining, (t-) undefining or (t#) changing a value in a tuple,