ITERATIVE CALIBRATION METHOD FOR A DIRECT NEURAL INTERFACE USING A MARKOV MIXTURE OF EXPERTS WITH MULTIVARIATE REGRESSION

Abstract

This invention relates to a method of calibrating a direct neural interface with continuous coding. The observation variable is modelled by an HMM model and the control variable is estimated by means of a Markov mixture of experts, each expert being associated with a state of the model.

During each calibration phase, the predictive model of each of the experts is trained on a sub-sequence of observation instants corresponding to the state with which it is associated, using an REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares) regression model.

A second predictive model giving the probability of occupancy of each state of the HMM model is also trained during each calibration phase using an REW-NPLS regression method. This second predictive model is used to calculate Markov mixture coefficients during a later operational prediction phase.

Claims

1. A method of calibrating a direct neural interface that will receive a plurality of electrophysiological signals acquired by means of a plurality of sensors, during a plurality of observation windows associated with observation times, and to provide command signals describing a trajectory to be made, said electrophysiological signals being preprocessed to obtain an observation tensor (X.sub.t) at each observation time (t), changes in the observation tensor being modelled by a hidden Markov Model (HMM), said interface using a mixture of a plurality K of experts, each expert (E.sub.k) being associated with a hidden state (k) of the HMM model and being defined by a multi-linear predictive model, predictions of each expert being combined by means of mixture coefficients (.sub.k.sup.t, k=1, . . . , K) to provide an estimate (Y.sub.t) of a control tensor representing said command signals, wherein: said calibration is made during a plurality of calibration phases, at predetermined instants of said trajectory, each calibration phase (u) corresponding to a plurality (L) of successive observation times, and making use of a tensor X.sub.u representing an observation tensor at said plurality of successive times, a tensor Y.sub.u representing a set command tensors at these same times and a matrix Z.sub.u giving states of the HMM model at these same times, and in that each calibration phase comprises: a step (a) in which tensors X.sub.u Y.sub.u, are extracted using the matrix Z.sub.u, observation tensors X.sub.u.sup.k and command tensors Y.sub.u.sup.k, k=1, . . . , K, relative to the different states (k) of the HMM model; a step (b) in which for each expert E.sub.k, tensors and are calculated corresponding to tensors X.sub.u.sup.k and Y.sub.u.sup.k respectively, after being centred and normalised, then a covariance tensor of tensor and a cross-covariance tensors of tensors and being modified by adding covariance and cross-covariance tensors and respectively derived from a previous calibration phase, weighted by a forget factor, .sup.k; a step (c) in which the multi-linear predictive expert models are updated E.sub.k, k=1, . . . , K using a multivariate regression of partial least squares with exponential recursive weighting (REW-NPLS) regression, starting from covariance and cross-covariance tensors modified in the previous step; a step (d) in which the tensor corresponding to the tensor X.sub.u, after being centred and normalised, the covariance tensor of tensor and the cross-covariance tensor of .sub.u and Z.sub.u are calculated, the covariance tensor of tensor and cross-covariance tensor of .sub.u and Z.sub.u being modified by adding to them the covariance tensor of .sub.u1 and the cross-covariance tensor of .sub.u1 et Z.sub.u1, derived from the previous calibration phase, weighted by a forget factor ; and a step (e) in which a second multi-linear predictive model is updated giving a state vector ({circumflex over (z)}.sub.t) of the HMM model as a function of the centred and normalised input tensor, .sub.u, the components of the state vector providing possibilities that the interface is in each of the different states k=1, . . . , K at an observation time, said update being made by an REW-NPLS regression starting from covariance and cross-covariance tensors modified in the previous step, the mixture coefficients (.sub.k.sup.t, k=1, . . . , K) of different experts then being obtained using the second multi-linear predictive model.

2. The method of calibrating a direct neural interface according to claim 1, wherein in each calibration phase, a number of transitions is determined for each pair of states observed during this calibration phase, and the transition matrix of the HMM model is updated using these numbers and the transition matrix of the HMM model calculated during the previous calibration phase, weighted by a forget factor.

3. The method of calibrating a direct neural interface according to claim 1, wherein in step (b), a modified covariance tensor cov(.sub.u.sup.k, .sub.u.sup.k) of the observation tensor for state k, centred and normalised, .sub.u.sup.k, is calculated by cov(.sub.u.sup.k,.sub.u.sup.k)=.sup.k cov(.sub.u1.sup.k,.sub.u1.sup.k)+.sub.u.sup.kx.sub.1.sub.u.sup.k in which x.sub.1 is a tensor product for a first mode, and in that the modified cross-covariance tensor cov(.sub.u.sup.k, .sub.u.sup.k) in which .sub.u.sup.k is a command tensor related to the centred and normalised state k is calculated by cov(.sub.u.sup.k, .sub.u.sup.k)=.sup.k cov(.sub.u1.sup.k, .sub.u1.sup.k)+.sub.u.sup.kx.sub.1.sub.u.sup.k.

4. The method of calibrating a direct neural interface according to claim 3, wherein in step (c), for each expert E.sub.k, k=1, . . . , K, an optimal rank, f.sub.opt.sup.k, in a latent variables space for the REW-NPLS regression is determined as that which minimises an error err.sub.u.sup.k,f.sup.k=.sup.k err.sub.u1.sup.k,f.sup.k+Y.sub.u.sup.k.sub.u1.sup.k,f.sup.k in which err.sub.u.sup.k,f.sup.k is a prediction error of the command tensor related to state k, by the expert E.sub.k during the calibration phase u for a rank f.sub.k=1, . . . , F in which F is a predetermined maximum rank, and .sub.u1.sup.k,f.sup.k is a prediction of Y.sub.u.sup.k by the expert E.sub.k, obtained by the multi-linear predictive model of the same expert during the previous calibration phase, for the same rank.

5. The method of calibrating a direct neural interface according to claim 1, wherein in step (d), a modified covariance tensor cov(.sub.u, .sub.u) of the observation tensor, centred and normalised, .sub.u, is calculated by cov(.sub.u, .sub.u)=.Math.cov(.sub.u, .sub.u)+.sub.ux.sub.1.sub.u in which x.sub.1 is a tensor product for a first mode, and in that the cross-covariance tensor cov(.sub.u, Z.sub.u) is modified by cov(.sub.u, Z.sub.u)=.Math.cov(.sub.u, Z.sub.u)+.sub.ux.sub.1Z.sub.u.

6. The method of calibrating a direct neural interface according to claim 5, wherein in step (e), an optimal rank, f.sub.opt, in a latent variables space for the REW-NPLS regression is determined as that which minimises an error err.sub.u.sup.f=.Math.err.sub.u1.sup.f+Y.sub.u.sub.u1.sup.f in which err.sub.u.sup.f is a prediction error of the command tensor during the calibration phase u for a rank f=1, . . . , F in which F is a predetermined maximum rank, and .sub.u1.sup.f is a prediction of Y.sub.u by a mixture of experts weighted by mixture coefficients obtained by the second multi-linear predictive model during the previous calibration phase, for the same rank.

7. The method of calibrating a direct neural interface according to claim 1, wherein the electrophysiological signals are electrocorticographic signals.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Other characteristics and advantages of the invention will become clear after reading a preferred embodiment of the invention, described with reference to the appended figures among which:

[0026] FIG. 1, described above, schematically represents a direct neural interface using a Markov mixture of experts, known in prior art;

[0027] FIG. 2 schematically represents a direct neural interface using a Markov mixture of experts according to one embodiment of the invention;

[0028] FIG. 3 is a flowchart representing a method of iteratively calibrating a direct neural interface using a Markov mixture of experts according to one embodiment of the invention.

DETAILED PRESENTATION OF PARTICULAR EMBODIMENTS

[0029] In the following, we will consider a direct neural interface with continuous decoding using a Markov mixture of experts as described in the introductory part.

[0030] The electrophysiological signals output from the different electrodes are sampled and assembled by data blocks, each block corresponding to a sliding observation window with width T. For example, the electrophysiological signals are sampled at a frequency of about 590 kHz, and the size of data blocks is 590 samples (T=1 s), and the offset T from one observation window to the next is 59 samples (T=100 ms). Each observation window is defined by an observation time (epoch) at which the window in question starts.

[0031] The electrophysiological signals are then preprocessed. In particular, this preprocessing may comprise elimination of the average taken on all electrodes, and a time-frequency analysis is then made on each observation window.

[0032] The time-frequency analysis can use a decomposition into wavelets, particularly Morlet wavelets.

[0033] A frequency smoothing or decimation can also be made on these results of the time-frequency analysis.

[0034] Thus, each observation window, or observation instant t, is associated with an order 3 tensor of observation data, resulting in the generation of an order 4 input tensor: the first mode corresponds to successive observation windows, the second mode corresponds to space, in other words the sensors, the third mode corresponds to time within an observation window, in other words the positions of the wavelets, and the fourth mode corresponds to the frequency, in other words to the number of frequency bands used for the decomposition into wavelets on an observation window.

[0035] More generally, the input tensor (or the observation tensor) will be order n+1, the first mode always being the mode relative to observation times (epochs). The input tensor (or observation tensor) is denoted X and its dimension is NI.sub.1 . . . I.sub.n.

[0036] Similarly, the trajectory of the imagined movement, observed or performed, is described by an output tensor (or command tensor), order m+1, denoted Y, with dimension NJ.sub.1 . . . J.sub.m, of which the first mode corresponds to successive times at which the commands will be applied (in general, this first mode also corresponds to the observation windows), the other modes corresponding to commands on the different effectors or different degrees of freedom of a multi-axis robot.

[0037] More precisely, the output tensor supplies N consecutive command data blocks, each of the blocks being used to generate command signals for the different effectors or degrees of freedom. Thus, it will be understood that the dimension of each data block can depend on the envisaged usage case and particularly on the number of degrees of freedom of the effector.

[0038] In the following, X.sub.t will be used to denote the observation tensor at time t. Consequently, this tensor is order n and its dimension is I.sub.1 . . . I.sub.n. It takes its values in a space X.sup.I.sup.1.sup. . . . I.sup.n in which is the set of reals. Similarly, Y.sub.t will be used to denote the command tensor at time t. Consequently, this output tensor is order m and dimension J.sub.1 . . . J.sub.m. It takes its values in a space Y.sup.J.sup.1.sup. . . . J.sup.m.

[0039] It will be assumed that the space X is formed from the union of a plurality K of regions, not necessarily discontiguous. In other words,

[00003]$\underset{}{X}=\underset{k=1}{\overset{K}{.Math.}}.Math.X_k$

in which X.sub.k, k=1, . . . , K are the regions in question.

[0040] The direct neural interface is based on a mixture of experts (ME), each expert E.sub.k operating on the elementary region X.sub.k of the input characteristics space and being capable of predicting the command tensor Y.sub.t at time t starting from the observation tensor, X.sub.t, when this tensor belongs to X.sub.k. In other words, each region X.sub.k has a corresponding command tensor prediction model Y.sub.t starting from the observation tensor X.sub.t. An expert E.sub.k can thus be considered as a multi-linear application of X.sub.k in Y.

[0041] It will also be assumed that each expert E.sub.k is associated with a hidden state k of a first order Markov Model (HMM). The different experts are combined by means of combination coefficients dependent on the hidden state at the time of the prediction. Definitively, starting from an input (or observation) tensor X.sub.t, the neural interface estimates the output (or command) tensor Y.sub.t, using:

[00004]$\begin{array}{cc}.Math.{\underset{\_}{Y}}_t=\underset{k=1}{\overset{K}{.Math.}}.Math.{}_k^t\left({\underset{\_}{}}_k.Math.{\underset{\_}{X}}_t+{\underset{\_}{}}_k\right)& \left(2\right)\end{array}$

in which .sub.k, .sub.k are a prediction coefficients tensor and an expert prediction bias tensor respectively, E.sub.k and .sub.k.sup.t is the weighting coefficient (also called the gating coefficient) for this expert and at the time t considered. The weighting coefficient .sub.k.sup.t is simply the conditional probability that the HMM model is in the state k knowing the previous input tensors X.sub.1:t=X.sub.1, X.sub.2, . . . , X.sub.t.

[0042] The set of coefficients and prediction biases, collectively called prediction parameters of the different experts, is designated by .sub.e={(.sub.k,.sub.k)|k=1, . . . , K}. The set of parameters for the combination of the different experts (including the parameters of the subjacent HMM model) is designated by .sub.g={A,{d.sub.k|k=1, . . . , K}, } in which A is the transition matrix between states, with size KK, in other words the matrix for which the elements a.sup.ij (independent of time by assumption of the HMM model) are defined by a.sup.ij=p(z.sub.t=j|z.sub.t1=i) in which z, represents the state of the model at time t and z.sub.t1 represents the state of the model at the preceding time; {d.sub.k|k=1, . . . , K} the parameters used to determine the conditional probability of observing the input tensor X.sub.t knowing the state z.sub.t=k, and is a vector with size K giving probabilities of occupancy of the different states at the initial time, in other words .sub.i=p(z.sub.t=i).sub.t=0.

[0043] FIG. 2 schematically represents a direct neural interface using a Markov mixture of experts according to one embodiment of the invention.

[0044] The electrophysiological signals output from the different electrodes 205 are processed in the processing module 210. This processing module performs sampling, optionally eliminates the average taken on all electrodes, and a time-frequency analysis is then made on each observation window. The time-frequency analysis can possibly be followed by frequency smoothing or decimation, as indicated above. Consequently, the processing module 210 outputs an input tensor X.sub.t at each observation window, characterised by an observation time t.

[0045] The direct neural interface also makes use of a hidden state machine based on an HMM model, 240, that can be in K possible states, and a plurality of experts, 220, each expert being associated with a hidden state.

[0046] Each expert E.sub.k makes an estimate Y.sub.t.sup.k of the output vector Y.sub.t and supplies it to a combination module, 230, that mixes the estimates according to expression (2). The combination coefficients (gating coefficients) are determined in the estimating module, 250, by .sub.k.sup.t=p(z.sub.t=k|X.sub.1:t), k=1, . . . , K, in which z.sub.t represents the state of the model at time t.

[0047] The estimate .sub.t gives the different kinematic movement parameters such as the position and/or the speed and/or the acceleration along the required trajectory.

[0048] The neural interface is calibrated firstly during an initialisation phase and then, during subsequent calibration phases, at given time intervals (on-line calibration).

[0049] These subsequent calibration phases are necessary due to the lack of stationarity of signals generated in the human brain. They can be done at regular intervals along the trajectory or successive trajectories to be executed.

[0050] A calibration phase u uses a plurality of input data blocks corresponding to a plurality L of successive times output from the observation times sequence. X.sub.u will denote the tensor with dimension LI.sub.1 . . . I.sub.n containing this plurality of input blocks. It also makes use of a plurality of output blocks giving kinematic set parameters, taken at the same times. Y.sub.u will denote the set tensor with dimension LJ.sub.1 . . . J.sub.m containing this plurality of input blocks.

[0051] The switches 261, 262 are in the closed position during a calibration phase. The calibration module 270 receives the input tensor X.sub.u, the set tensor Y.sub.u, and the matrix Z.sub.u with size KL, of which each of the columns represents the state of the machine at the time considered. More precisely, if at time the machine is in state z=k, all elements in the .sup.th column of the matrix Z.sub.u are null except for the element on row k that is equal to 1.

[0052] The calibration module uses X.sub.u, Y.sub.u, Z.sub.u to calculate the updated parameter sets .sub.e and .sub.g.

[0053] At the end of this calibration phase, the switches 261, 262 are open and the updated parameter sets .sub.e and .sub.g are output to a plurality of experts, 220, and to the machine 240, respectively.

[0054] Innovatively, the predictive model of each expert E.sub.k is trained using an REW-NPLS (Recursive Exponentially Weighted N-way Partial Least Squares) regression with a forget factor .sup.k (0<.sup.k<1) on calibration data blocks, X.sub.u, Y.sub.u, corresponding to the hidden state z.sub.u=k, denoted X.sub.u.sup.k and Y.sub.u.sup.k. More precisely, for each expert E.sub.k and for each calibration phase u, the data blocks X.sub.u.sup.k and Y.sub.u.sup.k are determined so that it can be trained in a supervised manner.

[0055] This is done as described in application FR-A-3061318, namely:

[0056] In a first step, the tensors X.sub.u.sup.k and Y.sub.u.sup.k are normalised starting from a number of observations said to be effective N.sub.u.sup.k=.sup.kN.sub.u1.sup.k+N.sub.u.sup.k in which N.sub.u.sup.k is the number of data blocks corresponding to the hidden state z.sub.u=k during the calibration phase u. The weighted sum s(X.sub.u.sup.k) and the weighted quadratic sum sq(X.sub.u.sup.k) are calculated:

[00005]$\begin{array}{cc}s\left({\underset{\_}{X}}_u^k\right)={}^k.Math.\underset{=1}{\overset{N_{u-1}^k}{.Math.}}.Math.x_{u-1;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k+\underset{=1}{\overset{N_u^k}{.Math.}}.Math.x_{u;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k& \left(3.Math.\text{-}.Math.1\right)\\ \mathrm{sq}\left({\underset{\_}{X}}_u^{{}_k}\right)={}^k.Math.\underset{=1}{\overset{N_{u-1}}{.Math.}}.Math.{\left(x_{u-1;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k\right)}^2+\underset{=1}{\overset{N_u^k}{.Math.}}.Math.{\left(x_{u;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k\right)}^2& \left(3.Math.\text{-}.Math.2\right)\end{array}$

in which and represent the elements of tensors X.sub.u1.sup.k and X.sub.u.sup.k. respectively. The average value

[00006]$\left({\underset{\_}{X}}_u^k\right)=\frac{s\left({\underset{\_}{X}}_u^k\right)}{N_u^k}$

and the standard deviation

[00007]$\left({\underset{}{X}}_u^k\right)=\sqrt{\frac{s.Math.q\left({\underset{}{X}}_u^k\right)-{\left(\left({\underset{}{X}}_u^k\right)\right)}^2}{N_u^k}}$

are then deduced. The elements of the centred and normalised input tensor, , are:

[00008]$\begin{array}{cc}{\tilde{x}}_{u;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k=\frac{x_{u;,i_1,.Math.\mathrm{.Math.}.Math..Math.i_n}^k-\left({\underset{\_}{X}}_u^k\right)}{\left({\underset{\_}{X}}_u^k\right)}& \left(4\right)\end{array}$

[0057] Similarly, the weighted sum s(Y.sub.u.sup.k) and the weighted quadratic sum sq(Y.sub.u.sup.k) are calculated:

[00009]$\begin{array}{cc}s\left({\underset{}{Y}}_u^k\right)={}^k.Math.\underset{=1}{\overset{N_{u-1}^k}{.Math.}}.Math.y_{u-1;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k+\underset{=1}{\overset{N_u^k}{.Math.}}.Math.y_{u;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k& \left(5.Math.\text{-}.Math.1\right)\\ \mathrm{sq}\left({\underset{}{Y}}_u^k\right)={}^k.Math.\underset{=1}{\overset{N_{u-1}^k}{.Math.}}.Math.{\left(y_{u-1;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k\right)}^2+\underset{=1}{\overset{N_u^k}{.Math.}}.Math.{\left(y_{u;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k\right)}^2& \left(5.Math.\text{-}.Math.2\right)\end{array}$

The average value

[00010]$\left({\underset{\_}{Y}}_u^k\right)=\frac{s\left({\underset{}{Y}}_u^k\right)}{N_u^k}$

and the standard deviation

[00011]$\left({\underset{}{Y}}_u^k\right)=\sqrt{\frac{\mathrm{sq}\left({\underset{}{Y}}_u^k\right)-{\left(\left({\underset{}{Y}}_u^k\right)\right)}^2}{N_u^k}}$

are then deduced, followed by the elements of the centred and normalised output tensor, :

[00012]$\begin{array}{cc}{\tilde{y}}_{u;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k=\frac{y_{u;,j_1,.Math.\mathrm{.Math.}.Math..Math.j_m}^k-\left({\underset{}{Y}}_u^k\right)}{\left({\underset{}{Y}}_u^k\right)}& \left(6\right)\end{array}$

[0058] The covariance tensor and the cross-covariance tensor are defined during the calibration phase u, starting from the centred and normalised tensors and :

[00013]$\begin{array}{cc}\mathrm{cov}\left(,\right)=.Math.{}_1.Math.& \left(7.Math.\text{-}.Math.1\right)\\ \mathrm{cov}\left(,\right)=.Math.{}_1.Math.& \left(7.Math.\text{-}.Math.2\right)\end{array}$

in which x.sub.1 designates the tensor product according to the first mode (temporal mode with index ) as described in the above-mentioned application FR-A-3061318.

[0059] These covariance tensors are modified taking account of the covariance tensors of the previous calibration phase, by weighting them with the forget factor .sup.k:

cov(,)=.sup.k cov(,)+x.sub.1(8-1)

cov(,)=.sup.k cov(,)+x.sub.1(8-2)

[0060] This calculation takes account of the data from a previous calibration phase, with a forget factor, to update the covariance tensors. The forget factors relative to the different experts can be chosen to be identical .sup.k=, k=1, . . . , K, in which is then the common forget factor. Alternatively, they can be chosen to be distinct so as to offer more prediction flexibility to the different experts.

[0061] Starting from the covariance tensors cov(, ) and cov(, ) working iteratively on rank f.sup.k=1, . . . , F in the latent variables space (in which F is a maximum predetermined number of dimensions in the latent variables space), we obtain a set of projection vectors {w.sub.u,1.sup.k,f.sup.k, . . . , w.sub.u,n.sup.k,f.sup.k}.sub.f.sub.k.sub.=1.sup.F with dimensions I.sub.1, . . . , I.sub.n respectively, a set of prediction coefficient tensors, {.sub.u.sup.k,f.sup.k}.sub.f.sub.k.sub.=1.sup.F, and a prediction bias sensor, {.sub.u.sup.k,f.sup.k}.sub.f.sub.k.sub.=1.sup.F, using the REW-NPLS method described in the above-mentioned application FR-A-3061318.

[0062] The optimum rank, f.sub.opt.sup.k, for each expert E.sub.k, during a calibration phase u is determined by calculating the prediction errors Y.sub.u.sup.k starting from X.sub.u.sup.k using the prediction coefficient tensors, {.sub.u.sup.k,f.sup.k}.sub.f.sub.k.sub.=1.sup.F, and the prediction bias tensors, {.sub.u.sup.k,f.sup.k}.sub.f.sub.k.sub.=1.sup.F obtained during the previous calibration phase:

err.sub.u.sup.k,f.sup.k=.sup.k err.sub.u1.sup.k,f.sup.k+Y.sub.u.sup.k.sup.k,f.sup.k(9)

in which .sup.k,f.sup.k is the estimate of Y.sub.u.sup.k obtained from the prediction .sub.u1.sup.k,f.sup.k and prediction bias coefficients .sub.u1.sup.k,f.sup.k from the previous calibration phase.

[0063] The rank f.sub.opt.sup.k leading to the minimum prediction error is then selected:

[00014]$\begin{array}{cc}{f}_{o.Math.p.Math.t}^k=\underset{f^k=1,.Math.\mathrm{.Math.}.Math.,F}{\arg .Math..Math.\min }.Math.\left({\mathrm{err}}_u^{k,f^k}\right)& \left(10\right)\end{array}$

and for the calibration phase, the prediction coefficients tensor and the prediction bias tensor corresponding to this optimal rank are selected, in other words the update is made at the end of the calibration phase u:

.sub.k=.sub.u.sup.k,f.sup.opt.sup.k (11-1)

.sub.k=.sub.u.sup.k,f.sup.opt.sup.k (11-2)

for each of the experts E.sub.k, k=1, . . . , K, independently. The result obtained is thus an estimate of the set of prediction parameters .sub.e making use of the REW-NPLS regression method applied to each of the experts.

[0064] Innovatively, the combination coefficients .sub.k.sup.t, k=1, . . . , K are also estimated using an REW-NPLS regression method adapted to discrete decoding, as described later.

[0065] During a calibration phase, the elements a.sup.ij of the transition matrix A are firstly updated as follows:

[00015]$\begin{array}{cc}{a}_u^{\mathrm{ij}}=.Math.a_{u-1}^{\mathrm{ij}}+\frac{v_u^{\mathrm{ij}}}{K}& \left(12\right)\end{array}$

in which a.sub.u.sup.ij (or a.sub.u1.sup.ij) are elements of the transition matrix, estimated during the calibration phase u (or u1) and is the number of transitions from state i to state j during the calibration phase u. This number of transitions is obtained starting from the matrix Z.sub.u. The lines in the transition matrix are then normalised such the sum of the probabilities of the transition from an initial state is equal to 1.

[0066] The transition matrix A thus updated is supplied by the transition matrix calculation module, 260, to the state estimating module, 240, of the HMM model.

[0067] z.sub.l defines a vector with size K for which the elements represent whether or not the machine is in state k=1, . . . , K at time t. Thus, if the machine is in state k at time t, only the kth element of z.sub.t will be equal to 1 and the other elements will be null. z.sub.t can be considered like a process with discrete values that vary in time according to a multilinear predictive model, the input variable of which is the observation tensor X.sub.t. This predictive model is trained under supervision during each calibration phase u, using the same principle as the predictive model of each expert. In other words:

{circumflex over (z)}.sub.t=BX.sub.t+b (13)

in which {circumflex over (z)}.sub.t expresses the probabilities {circumflex over (z)}.sub.k,t that the machine is in state k at time t, B is a prediction coefficients tensor with size KI.sub.1 . . . I.sub.n and b is a bias vector with size K.

[0068] The tensor B and the vector b are updated during each calibration phase u by means of an REW-NPLS regression with a forget factor i , (0<<1) based on observations X.sub.u and Z.sub.u (matrix with size KL). This forget factor may be identical to the common forget factor, when such a common forget factor is used to estimate prediction parameters for the different experts.

[0069] Working iteratively on rank f=1, . . . , F of the latent variables space starting from the covariance tensor cov(, ) in which is a centred and normalised version of X.sub.u, and the cross-covariance tensor cov(, Z.sub.u). These covariance and cross-covariance tensors are modified by:

cov(,)=.Math.cov(,)+x.sub.1(14-1)

cov(,Z.sub.u)=.Math.cov(,Z.sub.u)+x.sub.1Z.sub.u (14-2)

[0070] The REW-NPLS method provides a set of projection vectors {w.sub.u,1, . . . , w.sub.u,n}.sub.f=1.sup.F, with dimensions I.sub.1, . . . , I.sub.n respectively, a set of prediction coefficient tensors, {.sub.u.sup.f}.sub.f=1.sup.F, and a set of prediction bias vectors, {b.sub.u.sup.f}.sub.f=1.sup.F. Rank f.sub.opt leading to the minimum prediction error is selected, then the tensor B and the vector b are updated by:

B=B.sub.u.sup.f.sup.opt (15-1)

b=b.sub.u.sup.f.sub.opt (15-2)

[0071] The tensor B and the vector b are supplied to the mixture coefficients estimating module, 250.

[0072] Starting from elements {circumflex over (z)}.sub.k,t of the predicted vector {circumflex over (z)}.sub.t, the state estimating module can calculate the conditional probabilities using the softmax function:

[00016]$\begin{array}{cc}p\left(z_t=k|{\underset{}{X}}_t\right)=\frac{\exp \left({\hat{z}}_{k,t}\right)}{\underset{i=1}{\overset{K}{.Math.}}.Math.\exp \left({\hat{z}}_{i,t}\right)}& \left(16\right)\end{array}$

[0073] This expression gives the probability that the machine is in state k at time t, knowing the input data at this time.

[0074] The mixing coefficients estimating module, 250, obtains these coefficients from Bayes rule:

[00017]$\begin{array}{cc}{}_k^t=p\left(z_t=k|{\underset{\_}{X}}_{1:t},\right)=\frac{p\left(z_t=k,{\underset{\_}{X}}_{1:t}\right)}{p\left({\underset{\_}{X}}_{1:t}\right)}=\frac{p\left(z_t=k,{\underset{\_}{X}}_{1:t}\right)}{\underset{i=1}{\overset{K}{.Math.}}.Math.p\left(z_t=i,{\underset{\_}{X}}_{1:t}\right)}& \left(17\right)\end{array}$

in which p(z.sub.t=i, X.sub.1:t) are obtained using the forward algorithm, namely:

[00018]$\begin{array}{cc}p\left(z_t=i,{\underset{\_}{X}}_{1:t}\right)=p\left({\underset{\_}{X}}_t|z_t=i\right).Math.\underset{{}^{j=1}}{\overset{K}{.Math.}}.Math.a^{\mathrm{ij}}.Math.{}_k^{t-1}& \left(18\right)\end{array}$

[0075] The mixture coefficients can be obtained by recurrence by combining expressions (17) and (18):

[00019]$\begin{array}{cc}{}_k^t=\frac{p\left({\underset{\_}{X}}_t|z_t=k\right).Math.\underset{{}^{j=1}}{\overset{K}{.Math.}}.Math.a^{\mathrm{kj}}.Math.{}_k^{t-1}}{\underset{i=1}{\overset{K}{.Math.}}.Math.\left(p\left({\underset{\_}{X}}_t|z_t=i\right).Math.\underset{j=1}{\overset{K}{.Math.}}.Math.a^{\mathrm{ij}}.Math.{}_i^{t-1}\right)}& \left(19\right)\end{array}$

[0076] The conditional emission probabilities, p(X.sub.t|z.sub.t=i) can be obtained using Bayes rule starting from a posteriori conditional probabilities p(z.sub.t=i|X.sub.t) and a priori probabilities p(z.sub.l=i|X.sub.t):

[00020]$\begin{array}{cc}p\left({\underset{\_}{X}}_t|z_t=i\right)=\frac{p\left(z_t=i|{\underset{\_}{X}}_t\right).Math.p\left({\underset{\_}{X}}_t\right)}{p\left(z_t=i\right)}\frac{p\left(z_t=i|{\underset{\_}{X}}_t\right)}{p\left(z_t=i\right)}& \left(20\right)\end{array}$

the term p(X.sub.t) being a multiplication coefficient common to all states.

[0077] The probabilities of occupancy of the different states p(z.sub.t=i) at time t are calculated from the initial probabilities of occupancy and the transition matrix A.

[0078] FIG. 3 is a flowchart representing a method of iteratively calibrating a direct neural interface using a Markov mixture of experts according to one embodiment of the invention;

[0079] Said calibration method works during calibration phases planned at determined instants, for example at regular time intervals throughout the trajectory. Each calibration phase u uses a plurality of input data blocks corresponding to a plurality L of successive times in the observation times sequence. It also makes use of a plurality of output blocks giving kinematic set parameters at the same times.

[0080] The plurality of input data blocks for the calibration phase u is represented by the input tensor X.sub.u and the plurality of output blocks during this same calibration phase is represented by the output tensor Y.sub.u.

[0081] It is also assumed that the states of the HMM machine are known during the calibration phase. The different states of the HMM machine can relate to different elements to be controlled to make the trajectory, for example different members of an exoskeleton. Furthermore, for each of these elements, different states can be considered for example such as an active state when the patient controls this element, an idle state when he does not control it, and possibly a preparation state immediately preceding the active state and immediately following the idle state, during which the characteristics of the neural signals are not the same as in the idle state nor the active state for this element. Thus, when the direct neural interface must control P elements to make the trajectory, the HMM machine will comprise 2.sup.P or even 3.sup.P states (if preparatory states are envisaged). The states of the HMM machine during the calibration phase can be represented by a binary matrix Z.sub.u with size KL, all values in each column in Z.sub.u being 0 except for one element equal to 1 indicating the state of the machine at the corresponding time.

[0082] The input tensor, X.sub.u, the output tensor, Y.sub.u, the states matrix, Z.sub.u, are given in 310 for the calibration phase u.

[0083] In step 320, the matrix Z.sub.u is used as a starting point to determine the tensors X.sub.u.sup.k and Y.sub.u.sup.k formed by the input data blocks and the output data blocks related to the state k=1, . . . , K, respectively. These tensors are extracted from the input tensor, X.sub.u, and the output tensor, Y.sub.u.

[0084] In step 330, for each expert E.sub.k, k=1, . . . , K the centred and normalised tensors and are calculated and are used to deduce the covariance tensor , and the cross-covariance tensor of and . These tensors are modified by adding the covariance tensor of .sub.1 and the cross-covariance tensor respectively, obtained in the previous calibration phase, weighted by a forget factor .sup.k. These tensors thus modified will be used as covariance and cross-covariance tensors during the next calibration phase.

[0085] In step 340, the prediction parameter tensors .sub.k, .sub.k of each expert E.sub.k are updated using an REW-NPLS regression, starting from the covariance and cross-covariance tensors modified in the previous step.

[0086] After step 340, we have an updated version of the set .sub.e of expert prediction parameters K.

[0087] In step 350, the elements of the transition matrix A are advantageously updated starting from the number of transitions between successive states observed during the calibration phase u. The number of transitions between successive states is obtained starting from the matrix Z.sub.u. The matrix A is then modified by adding to it the matrix A obtained during the previous iteration, weighted by a forget factor .

[0088] In step 360, the centred and normalised tensor is added, followed by the covariance tensor of and the cross-covariance tensor of and Z.sub.u. These tensors are modified by adding the covariance tensor of and the cross-covariance tensor and Z.sub.u1, respectively, obtained in the previous calibration, weighted by the forget factor . These tensors thus modified will be used as covariance and cross-covariance tensors during the next calibration phase.

[0089] In step 370, a multi-linear predictive model is trained giving a state vector of the HMM machine as a function of the input tensor, {circumflex over (z)}.sub.t=BX.sub.t+b, the components of the state vector {circumflex over (z)}.sub.t providing probabilities that the machine is in the different states k=1, . . . , K respectively at time t. More precisely, the prediction parameter tensor B, and the bias vector b are updated using an REW-NPLS regression, starting from the covariance and cross-covariance tensors modified in the previous step.

[0090] The components of {circumflex over (z)}.sub.t calculated during an operational phase are then used at each observation time t to obtain mixture coefficients .sub.k.sup.t=p(z.sub.t=k|X.sub.1:t) of the different experts.