Gilles Celeux

Jean-Baptiste Durand
email: baptiste.durand@imag.fr

Selecting Hidden Markov Model State Number with Cross-Validated Likelihood

Keywords: Hidden Markov Models, Model Selection, Cross-Validation, Missing Values at Random, EM Algorithm

The problem of estimating the number of hidden states in a hidden Markov model is considered. Emphasis is placed on cross-validated likelihood criteria. Using crossvalidation to assess the number of hidden states allows to circumvent the well documented technical difficulties of the order identification problem in mixture models. Moreover, in a predictive perspective, it does not require that the sampling distribution belongs to one of the models in competition. However, computing cross-validated likelihood for hidden Markov models for which only one training sample is available, involves difficulties since the data are not independent. Two approaches are proposed to compute cross-validated likelihood for a hidden Markov model. The first one consists of using a deterministic halfsampling procedure, and the second one consists of an adaptation of the EM algorithm for hidden Markov models, to take into account randomly missing values induced by crossvalidation. Numerical experiments on both simulated and real data sets compare different versions of cross-validated likelihood criterion and penalised likelihood criteria, including BIC and a penalised marginal likelihood criterion. Those numerical experiments highlight a promising behaviour of the deterministic half-sampling criterion.

Introduction

Finite mixture distributions provide powerful models for statistical learning (see [START_REF] Mclachlan | Finite Mixture Models[END_REF]. In the context of Signal Processing, hidden Markov models (HMM), which are finite mixture models with Markov regime, have been widely used for modelling time series with homogeneous zones, where the observed data have a similar distribution, associated to the states of an unobserved finite state Markov process. For instance, HMM models have been used successfully in various applications such as Speech Recognition (see for instance [START_REF] Rabiner | A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition[END_REF] or DNA sequence analysis (see for instance [START_REF] Churchill | Stochastic Models for Heterogeneous DNA Sequences[END_REF].

For a review on hidden Markov models and their applications, see for example Ephraïm and Merhav, 2002). An important but difficult question to be solved when dealing with a finite mixture model is to choose a relevant number of components or a number of hidden states in the HMM context. Many criteria or procedures have been proposed to answer this open question (see McLachlan and Peel 2000, chapter 6 for a recent state of the art).

The main problem occurring with mixture order identification is that the validity conditions for the likelihood ratio tests and related penalised likelihood criteria do not hold. In particular, the vector parameter λ in R r for the true number of components is not identifiable in a space of higher dimension r ′ > r. For instance, it has been proved that, for general HMMs, the likelihood ratio statistic is stochastically unbounded even for bounded parameters [START_REF] Gassiat | The likelihood ratio test for the number of components in a mixture with Markov regime[END_REF]. However, for an independent mixture model, for which the observed data are independent, it has been shown that on a practical ground the BIC criterion of [START_REF] Schwarz | Estimating the dimension of a model[END_REF] has often a satisfactory behaviour (see [START_REF] Roeder | Practical Bayesian Density Estimation Using Mixtures of Normals[END_REF][START_REF] Roeder | Practical Bayesian Density Estimation Using Mixtures of Normals[END_REF][START_REF] Fraley | Model-Based Clustering, Discriminant Analysis, and Density Estimation[END_REF][START_REF] Fraley | Model-Based Clustering, Discriminant Analysis, and Density Estimation[END_REF]. Moreover, this criterion has been proved to lead to a consistent estimation of the number of components of a mixture model in the independent case when the likelihood of the model is bounded (see [START_REF] Kéribin | Consistent estimation of the order of mixture models[END_REF]. BIC, though, is not proved to be consistent for HMMs in the case of general observation distributions. Furthermore, it has been shown that BIC underpenalise the likelihood when determining the number of change-points in change-point processes, a problem closely related to HMMs [START_REF] Zhang | A modified Bayes information criterion with applications to the analysis of comparartive genomic hybridization data[END_REF]. Actually such a possible behavior as been noticed for estimating the number of components in a mixture model (see, for instance, [START_REF] Biernacki | Assessing a Mixture Model for Clustering with the Integrated Completed Likelihood[END_REF] In the independent case, many other criteria have been proposed and experimented in various situations (see Mclachlan and Peel, 2000 chapter 6). For instance the ICL criterion of [START_REF] Biernacki | Assessing a Mixture Model for Clustering with the Integrated Completed Likelihood[END_REF] has been shown to be relevant to choose a number of mixture components leading to a clustering structure with the greatest evidence (McLachlan and Peel 2000, chapter 6). Another promising criterion is the cross-validated likelihood criterion which has been proposed and experimented by [START_REF] Smyth | Model selection for probabilistic clustering using cross-validated likelihood[END_REF] in the case of independent mixtures. An advantage of this criterion is that it seems to avoid some of the theoretical difficulties occurring with penalised likelihood criteria and unrealistic assumptions regarding the distribution of the data. There have been a lot of papers from both theoretical and practical viewpoints on choosing the number of components in an independent mixture model. But, as far as we know, there have been a few practical studies on the choice of the number of hidden states for HMMs despite the fact that some penalised likelihood criteria have been analysed in the HMM context (see [START_REF] Boucheron | Inference in hidden Markov models, chapter order estimation[END_REF] for a review). Cross-validation can easily be applied to HMMs when the training data consist of several independent sequences (as in speech recognition for instance). However, applying cross-validation for HMMs estimated on a single sequence is problematic.

The aim of the present article is twofold. Firstly, several ways are proposed to solve the difficulty of implementing the cross-validated likelihood criterion in the HMM context, where the dependencies between the observations have to be taken into account in a proper way. Secondly, numerical experiments on both simulated and real data sets are proposed to compare various criteria to select the number of hidden states in a HMM model. Those criteria are various versions of the cross-validated likelihood criterion, the well known AIC [START_REF] Akaike | Information theory as an extension of the maximum likelihood theory[END_REF] and BIC criteria, the ICL criterion and the penalised marginal likelihood criterion proposed and studied in [START_REF] Gassiat | Likelihood ratio inequalities with application to various mixtures[END_REF] for HMMs.

The article is organised as follows. Section 2 is devoted to the presentation and computation of different versions of the cross-validated likelihood criterion. The computation of some of those criteria involves some difficulties when some of the data are removed at random from the learning set of dependent data. Those difficulties are overpassed in the present paper. Numerical experiments to compare several criteria, including our designed versions of cross-validated likelihood criteria, are presented in Section 3. Subsection 3.2 is devoted to the presentation of Monte Carlo numerical experiments on simulated data. Subsection 3.3 is devoted to the presentation of a real example concerning heart rate variability data for sleeping neonates. The article is concluded with a short discussion section, and the computational details of Section 2 are given in the Appendix.

2 Choosing the number of states of a HMM using crossvalidated likelihood

As mentioned in the introduction, assessing the number of components in mixture models encounters theoretical difficulties. A way to bypass those difficulties is to make use of resampling procedures. For instance, McLachlan and Peel (1997) proposed a parameterised bootstrap procedure to the assessment of the P -value of the likelihood ratio test in testing [START_REF] Smyth | Model selection for probabilistic clustering using cross-validated likelihood[END_REF] proposed to choose the number of components in a mixture model for independent data by maximising the cross-validated likelihood. The cross-validation approach is natural in model selection because it aims at estimating the predictive performance of a model. For instance, the cross-validated likelihood provides an out-of-sample estimate of the Kullback-Leibler divergence between the actual distribution of the data and a model distribution [START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF]. An advantage of the cross-validation approach in model selection is that it does not impose the unrealistic assumption that the actual distribution of the data belongs to one of the models in competition as, for instance, AIC and BIC do (see [START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF][START_REF] Spiegelhalter | Bayesian measures of model complexity and fit (with discussion)[END_REF]. In the mixture context, an other advantage of cross-validated likelihood is that it circumvents the above-mentioned technical difficulties encountered with penalised likelihood criteria. In this context, using the cross-validated likelihood and especially the half-sampling likelihood has been proved to be a valuable procedure to assess the number of components in finite mixture modelling [START_REF] Smyth | Model selection for probabilistic clustering using cross-validated likelihood[END_REF]. In the HMM context, the implementation of cross-validated procedures for determining an appropriate number of hidden states given the data is facing a specific difficulty: removing observations at random from the training data sequence would break the Markovian dependence between the states of the hidden Markov process. A simple method would consists of partitioning the original sequence into long contiguous blocks and treating them as independent sequences. However, such a strategy is much likely to produce partial sequences where some states of the Markov process are never reached, which would prevent it to find a sensible number of states. Now if, fortunately, a lot of sequences are available, a straightforward cross-validation approach can be implemented. It would consists of using the cross-validation procedure on the sequences rather the observations. This approach has been used in [START_REF] Robertson | Downscaling of daily rainfall occurence over Northeast Brazil using a hidden markov model[END_REF].

K = K 0 components versus K = K 1 components. More recently,
However, in this paper, we consider the more general situation where few sequences of data are available. And, in this section, we describe two novel procedures for computing the cross-validated likelihood in the context of hidden Markov processes.

In the sequel, we consider models for sequential data, i.e. hidden Markov chains. However, this study can be easily extended to non-sequential processes if the observed data are conditionally independent given the hidden process, for example to hidden Markov trees.A hidden Markov chain consists of an unobserved state process S = (S

Computing cross-validated likelihood criteria for HMMs

When concerned with computing cross-validated maximum likelihood, parameter λ has to be estimated from an incomplete observed sequence, since part of the training data set has been removed. Furthermore the model assessment requires the computation of the likelihood using the remaining data, which also forms an incomplete sequence. To deal with this difficulty, the first approach we propose is a particular half-sampling procedure.

It consists of removing from the original training sequence respectively the odd and even indices. The key point is that the resulting processes are still hidden Markov chains, which can be identified through the usual forward-backward recursion of [START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains[END_REF], i.e.

the standard implementation of the EM algorithm for HMMs. Since this implementation is subject to underflow when n is moderately large, we use the smoothing algorithm of [START_REF] Devijver | Baum's forward-backward Algorithm Revisited[END_REF], which is numerically stable. This approach provides an easy, rapid and efficient way to implement half-sampling. It preserves the Markovian structure of the model in a simple way, but imposes a particular form of half-sampling. The second approach we consider deals with the more general situation of v-fold cross-validation, where d = n/v among n data point are removed at random at each of the v steps of the cross-validation procedure. This approach has been proposed by [START_REF] Zhang | Model selection via multifold cross validation[END_REF] for the selection of linear models. It can be dealt with the EM algorithm of Dempster et al.

(1977) and it leads to derive a forward-backward recursion for HMMs that is dedicated to data removed at random.

Half-sampling from the odd and even subchains

Half-sampling is a two-fold cross-validation procedure for which the learning and test data sets, of equal size n/2, are usually chosen at random in the actual data set. However, the indices of those two sub-samples can be chosen in a deterministic way. If the two subsamples are designed by considering respectively the odd and even indices of the original data set, supposed to follow a HMM distribution, the resulting processes are still HMMs.

More precisely, the theorem below holds.

Theorem 1 Distribution of the odd (resp. even) subchain.

Let (Y 1 , . . . , Y 2n ) be a hidden stationary and ergodic Markov chain with K hidden states and with parameter λ = (π, p, θ 1 , . . . , θ K ). Then the odd (resp. even) subchain

(Y 1 , Y 3 , . . . , Y 2n-1 ) = (Y ′ 1 , . . . , Y ′ n ) (resp. (Y 2 , Y 4 , . . . , Y 2n ) = (Y ′′ 1 , . . . , Y ′′ n )
) is a hidden Markov chain with K hidden states and with parameter (π, p 2 , θ 1 , . . . , θ K ).

Remark: Note that this result would not hold if the Markov chain was periodic. However, since the Markov chain is supposed to be stationary and ergodic, it cannot be periodic.

Proof: For the odd sequence, the result follows from the fact that the hidden process (S 1 , S 3 , . . . , S 2n-1 ) is straightforwardly a Markov chain with K states, with initial distribution π and with transition probability matrix p 2 . Using the conditional independence of (Y 1 , Y 3 , . . . , Y 2n-1 ) given (S 1 , S 3 , . . . , S 2n-1 ) leads to the expected result. The argument is quite analog for the even sequence. The only difference is that the initial distribution of the even chain is πp. In the case where the original hidden chain is assumed stationary, both hidden subchains are Markov chains with the same distribution, since in this case

πp 2 = πp = π.
Consequently, the parameters of the odd and even HMM are the same. Their parameter can be estimated straightforwardly from the EM algorithm and the loglikelihood can be computed as a by-product of the EM algorithm. In the half-sampling procedure the role of the odd and even subsequences are permuted. When the vector parameter λ is estimated using the odd (resp. even) subsequence, its loglikelihood is calculated on the even (resp. odd) subsequence using the forward recursion of [START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains[END_REF]. Thus, the resulting OEHS (for Odd-Even Half-Sampling) criterion to assess the number of hidden states of the HMM model is the sum of those loglikelihoods. Since the dominant term of the complexity for each EM iteration is 2nK 2 , this is also the dominant term of the complexity for this half-sampling procedure.

Multifold cross-validation procedure

The implementation of the multifold cross-validation in the general case amounts to delete at random part of the observations and to estimate the vector parameter λ and to com- The details of the EM algorithm derivation are given in the Appendix. The main results are the following:

• The E step of the EM algorithm amounts to compute the quantities P (S t = j, S t+1 = k|Y Obs = y Obs ) and P (S t = j|Y Obs = y Obs ) for each hidden states j and k, at each time t;

• This can be done using the following forward-backward recursion. The forward recursion is based on the quantities αt (j) = P ({Y u = y u } u∈Obs,u≤t , S t = j) and the backward recursion is based on the quantities βt (j

) = P ({Y u = y u } u∈Obs,t<u≤n |S t = j).
As shown in the Appendix, those recursions amount to replace f θ k (y t ) with the value one when Y t is removed, in the standard forward-backward algorithm of Baum et al.

(1970). Moreover, our recursion has appealing interpretation: if Y t is removed and if Y t-1 and Y t+1 are observed, applying equation ( 10) followed by equation ( 9) in the Appendix shows that

αt+1 (k) = j p jk l p lj αt-1 (l)f θ k (y t+1 ) = j l p jk p lj αt-1 (l)f θ k (y t+1 ).
Thus, it can be noticed that the recursion involves the coefficients of matrix p 2 . This is due to the transition between S t-1 and S t+1 being ruled by p 2 in the Markov chain S. In the same manner, it is easily proved by induction that if n ′ successive y t are removed between two observed y t , the forward recursion associated with those two observed y t involves the transition matrix p n ′ +1 . This is illustrated in Figure 1.

• The above implementation of the forward-backward recursion is subject to underflow when n is moderately large, as discussed in [START_REF] Devijver | Baum's forward-backward Algorithm Revisited[END_REF]. This statement generally applies to HMMs, not just to partially observed HMMs. For this reason, we derive in the Appendix an implementation that is a smoothing algorithm inspired by Devijver's algorithm.

Both implementations can be used to compute the loglikelihood of a given parameter for a partially observed HMM, using equation [START_REF] Devijver | Baum's forward-backward Algorithm Revisited[END_REF] of the Appendix for the first implementation, and ( 18) for the second one.

t = 2 t = 3 S 3 Y 1 t = n time t = 4 S 4 Y 4 S 1 Y n S 2 Y 2 Y 3 A 3 S n t = t 0 = 1
Figure 1: Deletion of observations with random indices within a hidden Markov chain. If the random variables are deleted from t = t 0 to t = t 0 + n ′ , the forward algorithm computes αt+n ′ +1 from αt by computing p n ′ +1 . This is illustrated here for n ′ = 2 and t 0 = 1.

• The M step, given in the Appendix, is slightly different from the standard M step of [START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains[END_REF], and has also a direct interpretation in the case of observation distributions in the exponential family.

Remark: As shown in the Appendix, the adaptation of the general multi fold crossvalidation scheme to the HMM context appears to be possible without too much difficulty.

However, great variability can be expected from such a resampling procedure since useful information to estimate transition probabilities is lost. Moreover, it is quite doubtful that a single run of a multi fold cross-validated procedure would lead to a stable and reliable assessment of the number of hidden states of a HMM. In [START_REF] Durand | Modèles à structure cachée : inférence, sélection de modèles et applications[END_REF], numerical experiments highlight the poor behaviour of such single runs of multi fold cross-validation and show that, typically, such procedures have a tendency to overestimate the number of states in a HMM. To reduce this great variability, it is highly recommended to make use of Monte Carlo repetitions of the general multi fold cross-validation scheme. Denoting M the number of Monte Carlo repetitions of the multi fold cross-validation procedure, this leads to compute the cross-validated loglikelihood of a model as the mean of the M cross-validated loglikelihood values derived from the M cross-validation procedures. The EM algorithm detailed in the Appendix can easily be extended to the context of several incomplete sequences. Thus, the cross-validation methodology that consists of deleting observations at random could easily be extended to the case where several sequences are available, though other cross-validation scheme would also be possible.

Numerical Experiments

The aim of this section is to give elements on the practical ability of cross-validated likelihood criteria to estimate in a relevant way the number of hidden states in a HMM.

In that purpose, the performances of penalised likelihood criteria AIC, BIC, ICL and a Penalised Marginal Likelihood (PML) criterion are compared from numerical experiments with the performances of two cross-validation procedures presented in the previous section to assess the number of hidden states. Before presenting those numerical results, the compared criteria and the experiment conditions are described in the next subsection. In Section 3.2, numerical experiments on simulated data are reported. In Section 3.3, an application of the HMM model to analyse heart rate variability of neonates during sleep is presented.

Criteria in competition and experiment conditions

The considered cross-validated likelihood criteria are the Odd/Even half-sampling procedure (OEHS) and the Monte Carlo half-sampling procedure (MCHS) described in Section where

D(Y ) = 2E[ln p(Y ) -ln p(Y |m, λm )]
is the expected deviance on a test sample Y . Assuming that the data arose from a distribution belonging to the collection of models in competition and using asymptotic arguments, Akaike (1974) proposed to estimate this difference with 2ν m where ν m is the number of free parameters of the model m. This leads to the so-called AIC criterion

AIC(m) = 2 ln p(y|m, λm ) -2ν m . (1)
An other point of view consists of basing the model selection on the integrated likelihood of the observed data in a Bayesian perspective (see [START_REF] Kass | Bayes factors[END_REF]. This integrated likelihood is

p(y|m) = p(y|m, λ m )π(λ m )dλ m , (2)
π(λ m ) being a prior distribution for parameter λ m . A classical asymptotic approximation of the logarithm of the integrated likelihood is the BIC criterion of [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. It is

BIC(m) = ln p(y|m, λm ) - ν m 2 ln(n). (3)
This approximation needs regularity conditions on the likelihoods of the model collection M and is accurate when the prior distribution π(λ m ) is centered around the maximum likelihood estimate λm (see Raftery 1995). Notice that it has been argued that this formulation may only be appropriate in circumstances where it was really believed that one and only one of the competing models is in fact true (Bernardo and Smith 1994, chapter 6).

One of the interest of the BIC criterion as compared to the AIC criterion is that under regularity conditions, BIC is expected to be consistent while AIC is not. (For a general theoretical comparison between AIC and BIC, an interesting reference is [START_REF] Yang | Can the strengths of AIC and BIC be shared? A confict between model identification and regression estimation[END_REF]. The consistency of BIC estimator for HMM number of states is still far from being established [START_REF] Boucheron | Inference in hidden Markov models, chapter order estimation[END_REF]. However, under regularity conditions, [START_REF] Gassiat | Likelihood ratio inequalities with application to various mixtures[END_REF] has proved the consistency of the following PML criterion to estimate the number of hidden states of a HMM

PML(m) = n i=1 ln p(y i |m, λm ) - ν m 2 ln(n), (4)
where p(y i |m, λm ) denotes the marginal distribution of y i for model m. PML as the same form as BIC where the likelihood has been replaced with the pseudo likelihood of λ m .

When using a HMM model, it can be advantageous to choose the number of hidden states in order to get the mixture model giving rise to partitioning data with the greatest evidence. With that purpose in mind, [START_REF] Biernacki | Assessing a Mixture Model for Clustering with the Integrated Completed Likelihood[END_REF] proposed to approximate the integrated likelihood of the complete data (y, s) (or integrated completed likelihood), with a BIC-like approximation leading to the criterion

ICL(m) = ln p(y, ŝ | m, λm ) - ν m 2 ln n, (5)
where the sequence of missing states has been replaced by its most probable value ŝ, for parameter estimate λ, derived through the Viterbi algorithm (see for instance Ephraïm and Mehrav, 2002).

Experiment conditions

The simulation and estimation protocols are described here- completed. This phase is expected to reduce the dependence of the estimates with respect to the starting position. For each of the ten runs, the parameter value with highest likelihood is used as initial value for the EM algorithm, which is run for fifty iterations. The parameter value with highest likelihood, among the ten different runs, is chosen as initial position for a further long run of the EM algorithm, which is stopped either when 1,000 iterations are completed or when the relative loglikelihood increase is below 10 -6 . The MCHS criterion is computed using ten Monte Carlo repetitions.

Finally, it is important to note that the cross-validation procedures OEHS and MCHS can be initiated from the MLE derived from the whole observed data sequence. This strategy is expected to reduce the number of iterations required for convergence of the EM algorithm. However, this can lead to suboptimal estimates (compared to random initial values). Moreover, in a cross-validation procedure, the estimates to be assessed are not supposed to be dependent of the samples used to assess its performance. This is why we use the full sequence MLE as an eleventh initial parameter for the algorithm.

This mixed strategy empirically has led to the highest likelihood of the final parameter, compared to those based of either ten random initial values, or only the full sequence MLE.

The identification of partially observed HMMs has been performed with new Matlab routines 1 .

Simulated data

Thirty sequences of length 350 have been simulated, using a two-dimensional stationary Gaussian HMM model with three states, defined by the parameter λ 0 composed with p =   0.7 0.15 0.15 0.1 0.7 0.2 0.1 0.1 0.8

  µ 1 = 0.0 0.0 µ 2 = 2.0 1.9 µ 3 = 2.0 -1.9
and the variance matrices Σ 1 = Σ 2 = Σ 3 = Id, where Id refers to the identity matrix.

Figure 2 (a) displays the thirty simulated data sets on the same picture. When estimating the Gaussian HMM model, the variance matrices are assumed to be of the form γId, where the unknown scalar γ is independent from the Markov chain state. In this experiment, the number of hidden states K varies between K = 1 and K = 7. The frequencies of choosing a K value with each criterion are reported in Table 1.

Criterion

Number It appears that except PML and AIC, all the criteria select the generating model.

Criterion PML has a slight tendency to underestimate K while AIC has a marked tendency to select more complex models.

Further experiments, performed using a generating HMM model with five hidden states less separated than those defined by λ 0 , are now presented. Thirty sequences of length

      , µ 1 = 0 0 µ 2 = 2.7 2.7 µ 3 = 2.7 -2.7 µ 4 = 2.0 1.5 µ 5 = 2.0 -1.5 ,
and the variance matrices Σ 1 = . . . = Σ 5 = Id. Moreover, for the same model, thirty sequences of length 1,400 have been simulated. Those data sets are displayed in Figure 2 (b) and (c) respectively.

In the estimation procedure, the number of maximal iterations has been increased from 1,000 to 1,500 for the sequences of length 350, and to 2,000 for the sequences of length With a larger sequence of length 1400, all the criteria except ICL and PML, select most often the same number of hidden states as the generating model, as it appears in Table 3.

The most satisfactory results are obtained with criterion MCHS and to a smaller extend with criterion BIC. Criteria OEHS and AIC have a tendency to overestimate the number of hidden states, contrary to PML, which seems to highly underestimate the number of hidden states, as ICL does to a smaller extend. Actually, as shown in Figure 2(b) and (c), the five mixture components are poorly separated. From this point of view, the behaviour of ICL, which favours no overlapping solutions by its very nature, is sensible. Concerning the behaviour of PML, since the mixture components are overlapping, it seems difficult to diagnosis the presence of five hidden states without taking into account the Markovian dependency between data. However, the PML criterion does not.

To measure the asymptotic ability of the competing criteria to select the right number of hidden states, 30 sequences of length 15,000 have been simulated from the same model. Those results tend to highlight a correct behaviour of the three criteria for large samples, although BIC and ICL have a more marked maximum on the number of states of the generating model. It is worth noticing that for relatively small sample sizes, the chosen number of states with most criteria is smaller than the true number of states of the generating model, as it appears from Table 2. This behaviour is no surprising: for small sample sizes, a more parsimonious model can be expected to produce more stable prevision.

Heart rate variability in sleeping neonates

We consider the problem of heart rate variability analysis in sleeping neonates addressed in Consequently, we use a transformed signal where the number of heartbeats per second is computed. This transformation induces a loss of information, but leads to an easy segmentation of the signal into biologically interpretable periods. Since in our sample this quantity varies between 1 and 4, we model this discrete signal by a HMM with multinomial emission distributions, and with stationary state distribution. The parameter is estimated using the same procedure as in 3.2, except that three random values are drawn as initial parameters for the EM algorithm. Figure 3 presents the values of each selection criterion for one sequence of length 5,356.

Both BIC and ICL select a three-state model. The MCHS selects a two-state model.

Strictly speaking, the OEHS selects a five-state model, but the value of the criterion for this model is not much different from that for a two-state model, which would actually be selected by a modeller.

The MLE associated with a three-state model is the following: π = 0.73 0.15 0.12 ; p =   0.99 0.01 0.00 0.07 0.27 0.66 0.00 0.80 0.20

  ; θ1 = θ2 = θ3 =  
0.00 0.03 0.88 0.09 0.01 0.00 0.99 0.00 0.00 0.00 0.00 1.00

  ,
which means that P (Y t = 1|S t = 1) = 0.00, P (Y t = 2|S t = 1) = 0.03, P (Y t = 3|S t = 1) = 0.88,

P (Y t = 4|S t = 1) = 0.

Discussion

It has been proposed to assess the number of hidden states in a HMM with the crossvalidated likelihood. Owing to the dependence of the data sequence occurring with HMM models, the computation of cross-validated likelihood induced some difficulties. Those difficulties have been circumvented using a deterministic half-sampling scheme or solved using a new version of the EM algorithm to estimate the parameter of a HMM with data missing at random. This leads to define two cross-validation criteria of different nature to assess the number of hidden states in a HMM.

Both criteria, the so-called OEHS and MCHS have been compared with penalised likelihood criteria AIC, BIC, PML and ICL from numerical experiments. The conclusions of those numerical experiments are the followings.

• Criteria AIC, BIC and ICL seem to have a behaviour in the HMM context analogous to their behaviour in the independent mixture context (see McLachlan and Peel, 2000, chapter 6). Criterion AIC has a tendency to underpenalise the complexity of a model, ICL favours models that give rise to partitioning the data with the greatest evidence from the hidden states, and BIC seems to perform well if a HMM gives a reasonable representation of the observed process. However, for the heart rate variability data sets, BIC seems to overestimate the complexity of the models, as it could happen for real data sets. It provides a hidden representation of the data in three sets highly less useful than the two-state representation provided by the cross-validated criteria as it appears from Figure 4.

• In the HMM context, only the PML criterion has been proved to be consistent under regularities conditions [START_REF] Gassiat | Likelihood ratio inequalities with application to various mixtures[END_REF]. However, it seems that PML converges very slowly to the optimal solution. Moreover, in practical situations, it seems to have a high tendency to overpenalise the complexity of a HMM model when the sequence length is not very large. Actually, it can be remarked that the maximal pseudo likelihood is always smaller than the maximal likelihood for any HMM. It means that a good fit is less rewarded with the pseudo likelihood than with the likelihood.

• The cross-validation criterion MCHS seems to have a satisfactory behaviour. From our experiments, MCHS may have a slight tendency to select more complex models than BIC does on simulated datasets, but a more parsimonious one on the considered real dataset. Its main disadvantage is that it is much more time consuming than the other criteria. Only a two fold cross-validation criterion has been experimented; however, a Monte Carlo ten (for instance) fold cross-validation criterion would expectedly lead to analogous and more stable performance without producing a deterioration of the CPU time.

• Criterion OEHS can be regarded as a good surrogate to the MCHS criterion. It is highly less time consuming and seems to have analogous performance.

It can be argued that periodic subsampling schemes (as OEHS) will not work when the hidden process is also periodic, which is a very particular case. For such chains, some diagonal coefficients of the transition matrix will be null. Thirty sequences were simulated. The results are given in Table 5. A less efficient initialisation procedure was used, which explains the differences with Table 1. These results

show that the performance of OEHS remains reasonable when the sojourn time into each state is very short.

Q(λ, λ (m-1) ) = E λ (m-1) [ln L Y O bs ,Y M is ,S (λ)|Y Obs = y Obs ].
It follows from equation ( 6) that

Q(λ, λ (m-1) ) = t j E λ (m-
E λ (m-1) [ln f θ j (Y t )I {St=j} |Y Obs = y Obs ] = E λ (m-1) [ln f θ j (Y t )|S t = j]P λ (m-1) (S t = j|Y Obs = y Obs )
since, given S t = j, Y t does not depend on any other random variable, and particularly on the Y Obs process. Thus, equation ( 7) can be written

Q(λ, λ (m-1) ) = t∈Obs j ln f θ j (y t )P λ (m-1) (S t = j|Y Obs = y Obs ) + t∈Mis j E λ (m-1) [ln f θ j (Y t )|S t = j]P λ (m -

E step

The E step of the EM algorithm consists of computing Q(λ, λ (m-1) ). From [START_REF] Churchill | Stochastic Models for Heterogeneous DNA Sequences[END_REF], this amounts to compute the quantities P λ (m-1) (S t = j, S t+1 = k|Y Obs = y Obs ) and P λ (m-1) (S t = j|Y Obs = y Obs ) for each hidden states j and k, at each time t. This can be done using the following forward-backward recursion. In the sequel, P λ (m-1) is denoted P for the sake of simplicity.

The forward recursion is based on the quantities αt (j) = P ({Y u = y u } u∈Obs,u≤t , S t = j). It is initiated at t = 1 with α1

(j) = P (Y 1 = y 1 , S 1 = j) = π j f θ j (y 1 ) if 1 ∈ Obs; P (S 1 = j) = π j if 1 ∈ M is.
Then the (α t (j)) j are computed inductively for increasing values of t.

P ({Y u = y u } u∈Obs ) = j αn (j) (11)
The backward recursion is based on the quantities βt

If t + 1 ∈ M is, βt (j) is computed by integrating P ({Y u = y u } u∈Obs,t<u≤n |S t = j) over y t+1 . From equation [START_REF] Durand | Modèles à structure cachée : inférence, sélection de modèles et applications[END_REF], this leads to βt (j) =

y t+1 k p jk βt+1 (k)f θ k (y t+1 ) dy t+1 = k p jk βt+1 (k) y t+1 f θ k (y t+1 ) dy t+1 = k p jk βt+1 (k), (13)
since βt+1 (k) does not depend on y t+1 .

The above implementation of the forward-backward recursion is subject to underflow when n is moderately large, as discussed in [START_REF] Devijver | Baum's forward-backward Algorithm Revisited[END_REF].

For this reason, we recommend the following implementation, which is a smoothing algorithm inspired by Devijver's algorithm. The corresponding recursions are based on the quantities α t (j) = P (S t = j|{Y u = y u } u∈Obs,u≤t ) and

β t (j) = P ({Y u = y u } u∈Obs,t<u≤n | S t = j) P ({Y u = y u } u∈Obs,t<u≤n | {Y u = y u } u∈Obs,u≤t ) . Thus α t (j) = αt (j) P ({Y u = y u } u∈Obs,u≤t ) (14)
The forward-backward recursion of the smoothing algorithm is the following. For each state k and each time t > 1, we define α ′ t (j) as

α ′ t (j) = αt (j) P ({Y u = y u } u∈Obs,u≤t-1 ) . ( 15
)
As a consequence from equation ( 14),

α t (j) = α ′ t (j)P ({Y u = y u } u∈Obs,u≤t-1 ) P ({Y u = y u } u∈Obs,u≤t ) . ( 16
)
Then from equations ( 9) and ( 10), we obtain the recursion

α ′ t+1 (k) = j p jk α t (j)f θ k (y t+1 ) if t + 1 ∈ Obs j p jk α t (j) if t + 1 ∈ M is.
From equation [START_REF] Gassiat | Likelihood ratio inequalities with application to various mixtures[END_REF],

and from equation ( 16), α t+1 (k) is computed as

α t+1 (k) = α ′ t+1 (k) l α ′ t+1 (l)
.

As a consequence from equation ( 17) applied recursively, the likelihood is given by

P ({Y u = y u } u∈Obs ) = t j α ′ t (j) (18)
The quantity β t (j) in the backward recursion is related to βt (j) as follows:

β t (j) = βt (j) P ({Y u = y u } u∈Obs,t<u≤n |{Y u = y u } u∈Obs,u≤t ) .

For each state k and each time t < n, we define β ′ t (j) as β ′ t (j) = β t (j)P ({Y u = y u } u∈Obs,u≤t+1 ) P ({Y u = y u } u∈Obs,u≤t ) .

Consequently, from equations ( 19) and ( 20), β ′ t (j) = βt (j)P ({Y u = y u } u∈Obs,u≤t+1 )

P ({Y u = y u } u∈Obs ) . (21)
Then from equations ( 12) and ( 13), it leads to the recursion

β ′ t (j) = k p jk β t+1 (k)f θ k (y t+1 ) if t + 1 ∈ Obs k p jk β t+1 (k) if t + 1 ∈ M is.
Using equations ( 17) and ( 20), β t (j) is computed as

β t (j) = β ′ t (j) l α ′ t+1 (l)
.

The quantities P (S t = j|Y Obs = y Obs ), denoted by ξ t (j), and P (S t = j, S t+1 = k|Y Obs = y Obs ), denoted by γ t (j, k), are derived as in [START_REF] Devijver | Baum's forward-backward Algorithm Revisited[END_REF]. From the classical equations ξ t (j) = αt (j) βt (j) P (Y Obs = y Obs ) and γ t (j, k) = αt (j)p jk f θ k (y t+1 ) βt+1 (k) P (Y Obs = y Obs ) ,

we derive the following expression for the smoothed probabilities:

ξ t (j) = α t (j)β t (j) and γ t (j, k) = α t (j)p jk f θ k (y t+1 )β t+1 (k)

l α ′ t+1 (l)
.

M step

The maximisation of the quantity Q(λ, λ (m-1) ) with respect to λ leads to the standard update formula for p pjk = t P λ (m-1) (S t = j, S t+1 = k|Y Obs = y Obs ) t P λ (m-1) (S t = j|Y Obs = y Obs ) .

The estimation formula for θ j depends on the family {f θ } θ∈Θ . The formulas of Redner

  pute the likelihood from two different incomplete sequences. The whole training sequence y of length n is decomposed at random into an observed subsequence y Obs and a removed subsequence y Mis where Obs represents the subset of {1, . . . , n} corresponding to the observed value indices, and where M is is the set {1, . . . , n}\Obs of the removed value indices. To estimate the parameter λ in such context, where two types of missing values exist, namely the hidden states and the removed observations, we resort to the EM algorithm.

2. 1 .

 1 2 with half training and half test points. Penalised likelihood criteria A classical approach to the model assessing problem consists of penalising the fit of a model by a measure of its complexity. A convenient measure of fit is the deviance of a model m ∈ M, which is d(y) = 2[ln p(y) -ln p(y|m, λm )] where p(y) denotes the true distribution of the observed data y = (y 1 , . . . , y n ), p(y|m, λ m ) denotes the distribution under the model m parameterised with λ m , and λm is the maximum likelihood estimate (MLE) of λ m . A common way of choosing a penalisation term is to evaluate how large would be the deviance difference on average over learning and test sets of same size. That is, the penalisation would be an estimation of nD(Y ) -E(d(Y ))

  after. A procedure now described, inspired from Biernacki et al. (2003), is expected to greatly reduce the initial position dependence of the EM algorithm. This procedure makes use of ten different initial parameter values chosen at random. For each initial position, fifty iterations of the stochastic EM à la Gibbs algorithm of Robert et al. (1993) are first

Clairambault

  et al.(1992). The signal has been modelled using a Gaussian HMM in[START_REF] Celeux | Estimation de chaînes de Markov cachées : méthodes et problèmes[END_REF]. Their aim was to identify several periods in the neonate sleep, which can be characterised by states of the heart rate variability and can be interpreted as quiet sleep, intermediate sleep and agitated (or so-called active) sleep. The problem of selecting the number of hidden states has not been considered yet, and the authors chose an a priori two-state HMM corresponding to calm and active sleep modelling. Though the hypothesis of a Gaussian emission distribution family leads to operational restoration of the periods, this model does not appear to fit the data accurately.

  Howerer, it is interesting to assess the behaviour of OEHS for diagonal coefficients close to zero. For this reason, further experiments have been performed in[START_REF] Durand | Modèles à structure cachée : inférence, sélection de modèles et applications[END_REF], using a three-state model with following transition matrices: matrix p = p 1 of Section 3.2, and

For any parameter λ and

  any occurrence of the hidden sequence s, the completed data loglikelihood of λ is defined byln L y O bs ,y M is ,s (λ) = jk I {st=j,s t+1 =k} + K j=1 ln π j I {st=j}(6)where I {.} denotes the indicator function. Let λ (m-1) be the parameter estimate at iteration m -1 of the EM algorithm. The function Q of λ , to be maximised in the M step of the EM algorithm, is the conditional expectation of the completed data loglikelihood for the current parameter value λ (m-1) ,

1 )

 1 (S t = j|Y Obs = y Obs ) + t j k ln p jk P λ (m-1) (S t = j, S t+1 = k|Y Obs = y Obs ) + j ln π j P λ (m-1) (S t = j|Y Obs = y Obs ).

  (j) = P ({Y u = y u } u∈Obs,t<u≤n |S t = j). It is initiated at t = n with βn (j) = 1. Then the ( βt (j)) j are computed inductively for decreasing values of t. If t + 1 ∈ Obs, βt (j) is obtained from the ( βt+1 (k)) k as in Baum et al. (1970), i.e. βt (j) = k p jk βt+1 (k)f θ k (y t+1 ).

P

  ({Y u = y u } u∈Obs,u≤t+1 ) P ({Y u = y u } u∈Obs,u≤t ) ,

Figure 2 :

 2 Figure 2: (a) Simulated data from 30 HMMs of length 350 with K = 3 hidden states; (b) Simulated data from 30 HMMs of length 350 with K = 5 hidden states; (c) Simulated data from 30 HMMs of length 1,400 with K = 5 hidden states.

29

 29

Figure 3 :

 3 Figure 3: Model selection for the heart rate variability sample. The values of the criteria OEHS, MCHS, BIC, PML and ICL are represented as a function of the number of hidden states K.

Figure 4 :

 4 Figure 4: Hidden state restoration for the heart rate variability dataset with the Viterbi algorithm (applied to a 200-point-wide window). Top: two-state HMM model. Bottom: three-state HMM model.

  , . . . , s n ) of S, the Y t are conditionally independent and Y t follows a distribution with density f θs t belonging to a parametric family {f θ } θ∈Θ .

	1 , . . . , S n ) with finite
	values {1, . . . , K} and an observed process Y = (Y 1 , . . . , Y n ) such that
	1. S is a homogeneous Markov chain, which is supposed stationary and ergodic, with
	stationary state distribution π = (π i ) 1≤i≤K and with transition probability matrix
	p = (p ij ) 1≤i,j≤K ;
	2. Given a realisation s = (s 1 The set of parameters of the hidden Markov chain model, denoted λ in the following, is
	usually estimated with the maximum likelihood method by the EM algorithm (Baum et
	al., 1970).

Table 1 :

 1 Frequencies of choosing a number of hidden states for the three-state simulated HMM with the criteria MCHS, OEHS, AIC, BIC, ICL and PML. The frequencies are computed over thirty simulated sequences of length 350.

			of selected hidden states
		1 2	3	4	5 6 7
	MCHS	--	96.67 %	3.33 %	---
	OEHS	--	100.00 % -	---
	AIC	--	50.00 %	50.00 % ---
	BIC	--	100.00 % -	---
	ICL	--	100.00 % -	---
	PML	-23.33 % 76.67 %	-	---

Table 2 :

 2 Frequencies in percentage of choosing a number of hidden states for the fivestate simulated HMM with the criteria MCHS, OEHS, AIC, BIC, ICL and PML. The frequencies are computed over thirty simulated sequences of length 350.

	Criterion		Number of selected hidden states
		1 2	3	4	5	6	7
	OEHS	--	50.00 % 26.67 % 16.67 % 6.67 % -
	MCHS	--	33.33 % 40.00 % 23.33 % 3.33 % -
	AIC	--	6.67 %	50.00 % 40.00 % 3.33 % -
	BIC	--	90.00 % 10.00 % -	-	-
	ICL	-6.67 %	76.67 % 16.67 % -	-	-
	PML	-30.00 % 70.00 % -	-	-	-
	From Table 2, it appears that no criterion selects the generating model in the majority
	of cases, and most criteria favours the three-state solution, which seems reasonable as
	shown in Figure 2(b). From this set of numerical experiments, only cross-validation criteria
	and AIC indicate the five hidden state solution as a plausible solution. Moreover, in

1,400. A posteriori, this was not necessary since for most sequences, convergence was achieved before 1,000 iterations. the present case, AIC seems to have the most satisfactory behaviour, but it is probably favoured here by its tendency to underpenalise the complexity of a model.

Table 3 :

 3 Frequencies in percentage of choosing a number of hidden states for the fivestate simulated HMM with the criteria MCHS, OEHS, AIC, BIC, ICL and PML. The frequencies are computed over thirty simulated sequences of length 1,400.

Table 4 :

 4 Frequencies in percentage of choosing a number of hidden states for the five-state simulated HMM with the criteria OEHS, AIC, BIC, ICL and PML. The frequencies are computed over 30 simulated sequences of length 15,000.The results are represented in Table4. The procedure MCHS, which is much more CPU time costly, had not been included in those large sample size numerical experiments.It appears from table 4 that OEHS, BIC and ICL selected five-state models in the every case. AIC selected a six-state model for two sequences, and PML selected a six-state model for one sequence (with values of the criteria very close to those of a five-state and of a three-state model) and a five-state model for the other sequences.

	Criterion	Number of selected hidden states
		1 2 3 4 5	6	7
	OEHS	----100.00 % -	-
	AIC	----93.33 %	6.67 % -
	BIC	----100.00 % -	-
	ICL	----100.00 % -	-
	PML	----96.67 %	3.33 % -

  09, etc. Note that the actual frequency of {Y t = 1} is approximately 3.7 × 10 -4 . which is vaguely defined. The HMM model aims at associating the hidden states with stages of sleep. With a three-state model, states 2 and 3 appear to be unstable, as shown by the hidden state restoration on Figure4and by the diagonal coefficients of the estimated transition probability matrix. Thus, a sojourn into state 2 or state 3 cannot be interpreted as a well-defined stage of sleep, not even as intermediate sleep. In contrast, the states of the two-state model are stable and are more likely to be interpreted as stages of sleep. Therefore, a two-state model seems to be preferable to a three-state model. In the former model, states 1 and 2 would correspond to quiet sleep (characterised by slow heart rate) and active sleep (high heart rate), respectively. The three-state model leads to the splitting of state 2 of the two-state model into two new states that are poorly separated. Oscillations in the state process are the consequence of selecting an additional unnecessary state, rather than actual oscillations in the stage of sleep.

	The MLE associated with a two-state model is:
	π = 0.65 0.35 ; p =	0.99 0.01 0.02 0.98	; θ1 = θ2 =	0.00 0.03 0.89 0.08 0.00 0.00 0.62 0.38	,

which essentially consists of merging states 2 and 3 of the three-state model. This is also highlighted by the hidden state restoration performed by the Viterbi algorithm providing the globally most likely sequence of hidden states, as shown in Figure

4

. The PML selects a one-state model, and is a decreasing function of the number of hidden states. Such a one-state model is clearly inappropriate to model periods in neonatal sleep. From a biological viewpoint, two or three stages of sleep are defined in neonates (unlike in adults, who have five stages of sleep): quiet sleep, active sleep, and sometimes intermediate sleep,

Table 5 :

 5 Frequencies of choosing a correct (=), too low (-) or too high (+) number of states for three-state simulated HMMs with the criteria OEHS, MCHS and BIC, using three transition matrices and thirty sequences of length 350.

	Appendix: EM algorithm for HMMs with deleted observations

  1) [ln f θ j (Y t )I {St=j} |Y Obs = y Obs ] ln p jk P λ (m-1) (S t = j, S t+1 = k|Y Obs = y Obs ) + j ln π j P λ (m-1) (S 1 = j|Y Obs = y Obs ). (7) This expression can be simplified by noting that if t ∈ Obs, i.e. if Y t is considered as observed, E λ (m-1) [ln f θ j (Y t )I {St=j} |Y Obs = y Obs ] = ln f θ j (y t )P λ (m-1) (S t = j|Y Obs = y Obs ), and, if t ∈ M is, i.e. if Y t has been removed from the training sample,

	+
	t	j	k

  If t + 1 ∈ Obs, αt+1 (k) is obtained from the (α t (j)) j in a standard fashion as in[START_REF] Baum | A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains[END_REF], i.e. ∈ M is, αt+1 (k) is computed by integrating P ({Y u = y u } u∈Obs,u≤t+1 , S t+1 = k)

			αt+1 (k) =	p jk αt (j)f θ k (y t+1 ).	(9)
				j
	If t + 1 over y t+1 . From equation (9), this leads to
	αt+1 (k) =	y t+1 j	p jk αt (j)f θ k (y t+1 ) dy t+1		(10)
	=	j	p jk αt (j)	y t+1	f θ k (y t+1 ) dy t+1 =	j	p jk αt (j),
	since αt (j) does not depend on y t+1 and since f θ k is a probability density function. The
	likelihood is given by

  can be easily extended to the dependent case. The natural parameterΦ j = E θ [T (Y )] is estimated byFor instance, for the Gaussian family with mean µ and variance Σ, this leads to the (j)(y t -μj )(y t -μj ) t + [

								(m-1)
								j
			ξ t (j)					.
			t
	updating formulas
	μj =	t∈Obs	ξ t (j)y t + [ t∈Mis ξ t (j)	ξ t (j)]µ	(m-1) j
	Σj =		t∈Mis	t ξ t (j)][Σ	(m-1) j	+ (µ	(m-1) j	-μj )(µ	(m-1) j	-μj ) t ]
	and Walker (1984) for independent mixtures in the exponential family of distributions
		f θ (y) =	b(y) a(θ)	e θ t T (y)

t∈Obs ξ t (j)T (y t ) + [ t∈Mis ξ t (j)]Φ t∈Obs ξ t t ξ t (j)

.

Acknowledgements

The author are thankful to J. Clairambault and P. Gonçalvès for helpful discussions concerning the application.


Dr Jean-Alexandre Long
email: jalong@chu-grenoble.fr

Jean-Luc Descotes

Christian Letoublon

Development of the Miniaturised Endoscope Holder LER (Light Endoscope Robot) for Laparoscopic Surgery

Introduction: Our presentation focuses on our 1 st experiments with an innovatively-designed robot endoscope holder for laparoscopic surgery, which has a limited size and cost.

Materials and Methods:

A compact robot (LER) placed on the patient's skin that can be used in the lateral and dorsal supine position was tested on cadavers and laboratory pigs in order to allow successive modifications. The current control system is based on voice recognition. The amplitude of vision is 360 degrees with an angle of 160 degrees. Twenty three procedures were performed.

Results:

The tests made it possible to bring the prototype forward on a variety of aspects, including reliability, steadiness, ergonomics and dimension. The ease of installation, which takes only five minutes and the easy handling of the robot made it possible for 21 out of 23 procedures to be performed without the need of an assistant. Voice recognition was chosen as the control system.

Conclusion:

The LER robot is a robotised camera holder guided by the surgeon's voice that can eliminate the need for an assistant to hold the camera during laparoscopic surgery. The ease of installation and manufacture should make it an effective and inexpensive system in the lateral and dorsal supine positions. Ongoing randomized, prospective clinical trials will soon validate this robot prior to marketing.

Introduction

In laparoscopy, static camera-holding arms do exist, but do not allow the necessary changes in field of vision, particularly on zoom movements [START_REF] Moran | Robotic surgery: urologic implications[END_REF] . Light endoscope-holder robots were introduced to do away with the need for an assistant. AESOP®, EndoAssist ® and LapMan® are three of the systems currently on the market. They have been proven effective in ensuring a stable image and reducing operation times [START_REF] Kavoussi | Comparison of robotic versus human laparoscopic camera control[END_REF][START_REF] Mettler | One year of experience working with the aid of a robotic assistant (the voice-controlled optic holder AESOP) in gynaecological endoscopic surgery[END_REF] but have a limited distribution range, due to their large dimension and high price [START_REF] Aiono | Controlled trial of the introduction of a robotic camera assistant (EndoAssist) for laparoscopic cholecystectomy[END_REF][START_REF] Ballester | Comparison of task performance of robotic camera holders EndoAssist and AESOP[END_REF] . The aim of TIMC-GMCAO Laboratory was to create a miniaturised light endoscope-holder robot, offered at a lower cost. We are presenting the results of the trials carried out with the LER prior to clinical assessment.

Materials and Methods

The LER (light endoscope robot) developed by TIMC-GMCAO Laboratory consists of a compact camera-holder robot (figures 1 and 2) resting directly on the patient's abdomen, and an electronic box containing the electricity supply and robot controllers [START_REF] Berkelman | Development of a compact-cable driven laparoscopic endoscope manipulator[END_REF][START_REF] Berkelman | A compact, compliant laparoscopic endoscope manipulator[END_REF][START_REF] Berkelman | Design, Control and Testing of a Novel Compact Laparoscopic Endoscope Manipulator[END_REF] . The robot moves an optical trocar (tested diameter : 10 mm) making it possible guide and move the endoscope, the control of the zoom being included. Three miniaturised back-driveable motors are integrated. One motor is used to control the endoscope's insertion depth (motor 1) (figure 3), the second enables the core needle to rotate on an axis (motor 2) (figure 4), and the third enables endoscope pan-tilt (motor 3) (figure 5). It is possible to move the robot manually when the motors are off.

Gears and spoked wheels are used to control the rotations and a cable wound around a pulley device and held by a spring, makes it possible to control the endoscope's insertion depth (figure 2). The LER has three ranges of movement: 360° rotation around the vertical axis, an 80° tilt from the vertical position; and a horizontal incline of the endoscope 20 cm along its axis. A patent was filed in May 2002. A pedal, a microphone-equipped portable PC and vocal recognition software (Microsoft SpeechStudio®) are connected to the control box.

Vocal recognition can be replaced by remote control. The robot is connected to the control box by two cables.

The compact architecture was designed so as to not hamper the surgeon's movement and enable rapid set-up, while keeping costs low [START_REF] Berkelman | Development and testing of a compact endoscope manipulator for minimally invasive surgery[END_REF] . A fully sterilisable robot place directly on the patient's abdomen was ultimately chosen to save time in putting together the robot. The return to the highest position is enabled by a spring.

Motor 1

Results

In two procedures, it was necessary to stop the robot and resume using camera control because of two robot breakdowns: the base ring came out of position; and once one motor overheated. The accidents were corrected on the following version of the LER.

-Robot steadiness:

The robot had to be fixed to the operating table using an articulated arm (figure 6).

The tests on the pelvi-trainer had given hope that it would be possible to lay the robot on the patient, but after several hours of surgery, the device left marks on the animals' skin. Clamping the robot made it possible to correct that problem, while also improving the stability of the image. The articulated arm does present a drawback in that it can extend over or near the operating area. As a result, it is necessary to reflect on how to position the LER before making an incision, so that the set-up does not hinder the movement of the instruments.

-Command interface:

The interface customarily used was the voice recognition system, as it proved reliable in recognising orders, even with surgeons who had no experience using the robot. The sound environment, which was particularly noisy at times, did not create any problems, as described for other systems [START_REF] Ballester | Comparison of task performance of robotic camera holders EndoAssist and AESOP[END_REF] . Training was preferable beforehand so as to ensure that the software would recognise a specific voice.

Manual commands using a mini-keypad clipped onto an instrument can also be used with high precision and good ergonomics. However, the need to move the keypad onto another instrument when changing tools proved a hindrance.

-Set-up and removal time:

The quick and easy set-up for the robot varied from 5 to 40 minutes, as different systems were tested, and depending on experience. During the last procedures, the set-up times were the shortest. The positioning of the articulated arm used to clam the robot had to be given consideration before the procedure, so that it would not affect the positioning of the trocars. Removal was very quick and a test conversion to laparotomy on a laboratory pig, following healing of the splenic pedicule during a splenectomy, showed that the system can be easily removed in 30 seconds.

-

Magnitude of intra-abdominal field of vision:

The 360° movement made it possible to explore the entirety of the abdominal cavity.

The robot's speed of movement was satisfactory (75° per second, 80 mm/s), making it possible to focus on a point of interest quickly in continuous mode. The discontinuous movement function made it possible to change the angle of vision by just a few degrees, thereby enabling small, finer readjustments during already-meticulous movements. We were hampered by the inadequate insertion depth of the endoscope during a uretro-vesical anastomosis, performed as part of radical prostatectomy. The lack of zoom to the very end of the pelvis was due to the length of the compact spring, which reduced the opportunity to move the endoscope in depth (figure 7). This problem has been resolved.

-

Comfort of use

The surgeons using the robot were satisfied with the ergonomics and comfort of the LER.

During the experiments, the standard headset was replaced by a bluetooth headset to reduce the amount of cables.

-

Robot dimensions and trocar positioning

The advantage with the LER lies in its compactness. The assistants were not hampered by the robot when setting up or moving around during the procedure. The clamp did not hinder the surgeon's movements either. As regards the position of the trocars, the 11-cm-diameter base did not require any change in trocar position during the procedures that we performed using laparoscopy. In contrast, in lateral supine position, during the first procedure, we were hampered by the clamp arm, which had been placed on a trocar insertion point.

-Position of patient.

The initial robot had been created for use in dorsal supine position (figure 8). The need to clamp the robot extended the potential patient positions to include the lateral supine position (figure 9). The stability of the whole system was satisfying. The diameter of the robot's base made it uneasy to introduce the trocars during a lomboscopy, which does not appear to be an indication. With the transperitoneal approach, we did not encounter any problems. The trials on animals and cadavers were compelling in this position.

The set up in dorsal supine position did not raise any particular problems.

-System Reliability Reliability was the main problem during the trials, but has since been resolved. 2 motors had to be changed as they had deteriorated following an electrical trial.

Discussion

Through use of the LER, the robot proved reliable on cadavers and laboratory pigs. The reliability problems that appeared were solved during the final tests. The patient security issues appear to have been overcome and it appears possible to set the robot directly on the patient. The surgical procedures were never forced to deviate from conventional laparoscopy.

The most noteworthy benefits were the stability of the useful image during anastomosis and the possibility of operating alone. However, this required robot set-up times ranging from 5

to 40 minutes, depending on the team's level of experience.

The recognised benefits of endoscope-holder robots are the ability to do without an aid, the image stability, the lesser degree of fatigue for the surgeon and the lower frequency of camera soiling [START_REF] Long | Use of robotics in laparoscopic urological surgery: state of the art[END_REF] . However, their large dimensions and, above all, their price, have slowed down their spread. The aim is to produce an affordable system that can be cost-effective, as it does away with the need for an operating assistant and enables quick movements. For this reason, it is necessary that the system be set up quickly, as it appears difficult that the robot should be able to provide any improvement in procedure time compared to an assistant-facilitated procedure, even if the literature does report lower times with EndoAssist 4 and AESOP [START_REF] Long | Use of robotics in laparoscopic urological surgery: state of the art[END_REF] .

It is difficult to estimate in advance the price of the LER, as the price of a prototype does not include the development and marketing costs.

The contribution of the LER as compared to existing systems lies above all in its compact nature, which makes it possible not to clutter the operating table and room area, while not hampering the assistant's movements. Trocar position was not hampered by the robot's base diameter during laparoscopy.

Feasibility has been affirmed by the pre-clinical tests and many changes have taken place. Its ability to withstand time and repeated sterilisation procedures is yet unknown.

Conclusion

The LER is a light endoscope-holder robot with miniature design, set on the patient's abdomen, in order to make the system less cumbersome and costly than the robots currently on the market.

Pre-clinical trials have shown that LER use is feasible in laparoscopic digestive and urological surgery on cadavers and animals.

Clinical trials will tell whether the system is suited to man.

Figure 1 :

 1 Figure 1: LER, or light endoscope robot

Figure 2 :

 2 Figure 2: LER robot scheme

Figure 3 :

 3 Figure 3: Action of Motor 1 enabling zoom movements as the cable winds around the pulley.The return to the highest position is enabled by a spring.

Figure 4 :Figure 5 :

 45 Figure 4: Encoscope rotation movements around trocar axis (motor 2)

Figure 6 :

 6 Figure 6: Robot firmly set

Figure 7 :Figure 8 :

 78 Figure 7: Zoom (spring and cable)


E Marenić
email: marenic@insa-toulouse.fr

D Brancherie
email: delphine.brancherie@utc.fr

M Bonnet
email: mbonnet@ensta.fr

Asymptotic analysis based modeling of small inhomogeneity perturbation in solids: two computational scenarios

Keywords: defect, asymptotic analysis, elastic moment tensor

Asymptotic analysis based modeling of small inhomogeneity perturbation in solids: two computational scenarios

Eduard Marenić, D. Brancherie, M

Motivation, introduction and problem definition

The role played by defects in the initiation and development of rupture is crucial and has to be taken into account in order to realistically describe the behavior till complete failure. The difficulties in that context revolve around (i) the fact that the defect length scale is much smaller than the structure length scale, and (ii) the random nature of their position and size. Even in a purely deterministic approach, taking those defects into consideration by standard models imposes to resort to geometrical discretisations at the defect scale, leading to very costly computations and hindering parametric studies in terms of defect location and characteristics.

Our current goal is to design an efficient two-scale numerical strategy which can accurately predict the perturbation in terms of stress caused by an inhomogeneity in elastic (background) material. To make it computationally efficient, the analysis uses only a mesh for the defect-free structure, i.e. the mesh size does not depend on the (small) defect scale.

We consider a linearly elastic body occupying a smooth bounded domain Ω ⊂ R d (with the spatial dimensionality d = 2 or 3), whose boundary Γ is partitioned as Γ = Γ D ∪Γ N support a prescribed traction t and a prescribed displacement ū, while a body force density f is applied in Ω. On the basis of this fixed geometrical and loading configuration, we consider two situations, namely (i) a reference solid characterized by a given elasticity tensor C, which defines the background solution u, and (ii) a perturbed solid constituted of the same background material except for a small inhomogeneity whose material is characterized by C , which defines a perturbed solution u a . The aim of this work is to formulate a computational approach allowing to treat case (ii) as a perturbation of the background solution (i), in particular avoiding any meshing at the small inhomogeneity scale. This will be achieved by applying known results on the asymptotic expansion of the displacement perturbation with respect to the small characteristic size a of the inhomogeneity to case (ii).

Small inhomogeneity asymptotics

The difference of the mentioned displacement fields gives a displacement perturbation v a := u a -u. An asymptotic analysis of v a with respect to the characteristic defect size a provides a way to evaluate the influence of the location, size, shape and material characteristics of defects on the solution u a . We focus in this work on the outer expansion given by [START_REF] Ammari | Polarization and moment tensors. With applications to inverse problems and effective medium theory[END_REF]Thm. 11.4]

v a (x) = -∇ (1) G(z, x) : A(B, C, ∆C) : ∇u(z)a d + o(a d ), x = z, (1)
where the leading contribution is given as multiplicative decompsition in terms of three key 'ingridients': elastostatic Green's tensor G, elastic moment tensor (EMT) associated with the inhomogeneity A, and background strain (see [START_REF] Marenić | Application of asymptotic analysis to the two-scale modeling of small defects in mechanical structures[END_REF]).

G(ξ, x) represents the displacement field being a response at ξ ∈ Ω of the background body subjected to a unit point force applied at x with homogeneous boundary conditions. Using known symmetry relationship and decomposition G = G ∞ + G c to full-space and complementary parts, respectively, we propose two distinct computational scenarios: (a) fixed inclusion location z and varying evaluation point x, or (b) fixed evaluation point x and varying inclusion location z.

3 Computational scenarios to compute inclusion perturbation It turns out [START_REF] Marenić | Application of asymptotic analysis to the two-scale modeling of small defects in mechanical structures[END_REF] that for scenario (a) one needs to solve the d(d + 1)/2 problems (one for u and others for ∇G c ) to compute the perturbation anywhere in solid x ∈ Ω for virutally any eliptic inclusion. That is, it's shape, size and material can vary with virtually no additional cost.

In the inverse approach, scenario (b), the evaluation point x will be a point from critical zone. The perturbation in x will be computed for virutally any eliptic inclusion positioned anywhere in solid, z ∈ Ω. (Note that x = z for both scenarios).

Proposed computational scenarios can be used as numerical strategy for predicting the (stress) perturbation caused by small isolated inhomogeneities in elastic solids. Meshing the isolated inhomogeneity (with the attendant mesh preparation and computational costs) is avoided, the whole analysis relying on a mesh that is suitable for the defect-free configuration. Proposed approach will be illustrated on few numerical examples.