Yeser Meziani

Yann Morère

Amine Hadj-Abdelkader

Mohammed Benmansour

Guy Bourhis

Kalman Extended Smoother

Towards Adaptive and Finer Rehabilitation Assessment: a Learning Framework for Kinematic Evaluation of Upper Limb Rehabilitation on an Armeo Spring Exoskeleton

Keywords:

Providing specialized rehabilitation and tailoring the training process for patient's needs and according to recovery potentials has gained importance. To satisfy this need, a dynamic assessment of the performance of the recovery process is required. Assessing rehabilitation for the upper limb is often carried out with clinical subjective scales that do not satisfy these requirements. The use of technologies introduced several sensors into the devices used for rehabilitation and permitted the rise of kinematic assessments.

Kinematic measures provide an objective scale to follow up recovery during upper limb rehabilitation. The kinematics are still raw evaluations since they present insignificant effects if studied over short periods or on heterogeneous samples.

We propose a framework for modeling the trajectories as a means of encoding the specificity of the movement at every stage. The new technique permits detecting significant differences as soon as three training sessions became available.

We adopt an expectation-maximization algorithm and an optimization technique to encode the trajectories and the transition model from the acquired data. The framework enables us to encode in a Bayesian sense the observations from the patient and define six metrics to follow up on the progress of the movement quality. Statistical analysis of the results proved that these metrics are effective in tracking the evolution of the recovery. The results also established a strong discriminative property.

The proposed framework promises a finer scale of evaluation and extends the knowledge about kinematic assessment. This study's findings suggest that adopting these new metrics can help achieve more individualized patient care. It additionally promises to limit the amount of data needed to detect a significant change.

Introduction

The rehabilitation exercises are increasingly served by the latest technologies in the form of serious games in virtual reality environments, robotic exoskeletons, or assistive training robots which are proving to be highly beneficial rehabilitation tools [START_REF] Kwakkel | Effects of robot-assisted therapy on upper limb recovery after stroke: a systematic review[END_REF][START_REF] Aminov | What do randomized controlled trials say about virtual rehabilitation in stroke? a systematic literature review and metaanalysis of upper-limb and cognitive outcomes[END_REF].

This presence has come to reply to the numerous guidelines and recommendations that were set for the design of rehabilitation exercises [START_REF] Carr | Stroke Rehabilitation: Guidelines for Exercise and Training to Optimize Motor Skill[END_REF]. These recommendations aim to maximize the efficiency of the treatment while compensating for the shortage of qualified personnel and therapists delivering these services by opening doors for automation. This shortage becomes noticeable given the growing demographics of the elderly and the motor-impaired patients victims of certain pathologies such as stroke.

European Robot Road-map of 2010 also emphasized the urgency to provide robotic solutions for healthcare assistance that are intelligently capable of automating a large portion of the caregiving procedure [START_REF]Strategic Research Agenda for Robotics in Europe[END_REF].

One of the recurrent challenges for researchers tackling these requirements is the ability to individualize and adapt yeser.meziani@univ-tlemcen.dz (Y. Meziani) ORCID(s): 0000-0001-9787-8757 (Y. Meziani) the exercise to the evolution of the patient's recovery and his current performance. The performance is often affected by different factors and does not always reflect the actual fitness level of the patient.

Although these guidelines are recommendations, the use of ubiquitous technologies for rehabilitation training stresses the need to automate the process of assessment, diagnosis, and monitoring. This consists of developing a more inclusive patient-robot-therapist loop, where continuous feedback is present between all elements. Achieving this would enable the full exploitation of these tools while liberating the therapist to be able to provide care of more patients at once.

While surveying the literature, replying to these challenges has often been approached by modeling the patient's behavior during the rehabilitation training. This has been approached in different ways depending on the tools and the envisioned application context.

Many studies presented an intelligent agent to handle game difficulty for rehabilitation exercises based on virtual reality [START_REF] Brezany | Towards precision brain disorder rehabilitation[END_REF][START_REF] Alamri | A decision model of stroke patient rehabilitation with augmented reality-based games[END_REF][START_REF] Tsiakas | A multimodal adaptive session manager for physical rehabilitation exercising[END_REF][START_REF] Borghese | A cloud-based platform for effective supervision of autonomous home rehabilitation through exer-games[END_REF]. Both [START_REF] Taub | Autocite: automated delivery of ci therapy with reduced effort by therapists[END_REF][START_REF] Kan | The development of an adaptive upper-limb stroke rehabilitation robotic system[END_REF] are examples of approaches involving an assistive robot instead of a virtual environment. Their theoretical frameworks range between Bayesian framework, fuzzy inference, and Markov Decision Processes with Reinforcement learning. To retrieve the para-meters of the exercise they often rely on multi-modal sensing using RGB-D cameras. Otherwise recovering sensors' readings from the assistive robot were adopted for robot-based systems. Applications of telerehabilitation focus on building such assessment and monitoring models for online feedback for therapists and users as seen in [START_REF] Pirovano | Self-adaptive games for rehabilitation at home[END_REF][START_REF] Borghese | A cloud-based platform for effective supervision of autonomous home rehabilitation through exer-games[END_REF].

The authors of [START_REF] Kleine Deters | Therapeutic exercise assessment automation, a hidden Markov model approach[END_REF] employed a Hidden Markov Model (HMM) to follow and decide dynamically on the rehabilitation exercise parameters concerning the actual performance of the patient.

In [START_REF] Capecci | A hidden semi-markov model based approach for rehabilitation exercise assessment[END_REF], HMMs were used to follow on the execution of a rehabilitation exercise captured by an RGB-D camera and to provide feedback on the correctness of the execution.

A shared limitation to these approaches is that, while accounting for the environment's specific actions set i.e. game controls, they fail to generalize easily across rehabilitation systems and remain thus system-dependent.

Learning the kinematic model of a task appears to be a promising endeavor since trajectories encode for both performance and strategies, while also being easily portable between systems. To this end, we can find the process of learning the kinematic models of the patient during a specific task. These methods are often conceived for path planning for active exoskeleton control. We refer to [START_REF] Zhou | Learning control in robot-assisted rehabilitation of motor skills-a review[END_REF] for a detailed overview of these methods and techniques for rehabilitation applications.

Particularly, of interest in our context is the emphasis the authors put on the fact that the potential final model of the patient's performance remains a hidden outcome. They point out that the proposed modeling techniques need to account for motion limitations due to fatigue and stiffness among other factors influencing the performance of the patient. The fact remains that an active research endeavor is modeling the recovery process whereas the few approaches proposed within the rehabilitation context still lag in terms of clinical validity and general utility.

In [START_REF] Reinbolt | Determination of patient-specific multi-joint kinematic models through two-level optimization[END_REF], a musculoskeletal model of the patient is built by approaching the human body with a multi-joint kinematic model and estimating afterwards the parameters from the data. An extended review paper in [START_REF] Leardini | Kinematic models of lower limb joints for musculo-skeletal modelling and optimization in gait analysis[END_REF] provides more details on studies implementing similar models. These kinematic models define the joint structure of the body during a rehabilitation exercise from a morphological point of view to assess the movement conducted during training or estimate the interventions' effects.

Another work proposed in [START_REF] Mendoza-Crespo | An adaptable human-like gait pattern generator derived from a lower limb exoskeleton[END_REF] presents a learning technique for the kinematic model in the context of gait trajectories identifying centroids in the trajectories and interpolating them to model a step from multiple healthy users' demonstrations. Learning the Kinematic model as a discriminant representation of the user without considering the variability in his performance would consist an undermining factor to the generalizability of similar approaches.

Meanwhile, in [START_REF] Lauretti | Learning by demonstration for motion planning of upper-limb exoskeletons[END_REF] the authors used the Dynamic Movement Primitives (DMPs) to encode for the trajectories accomplished during the Activities of Daily Living (ADLs). The challenge that the DMP framework presents is the as-sumption of availability and knowledge of the model of the system which is indeed the case with robotic arms. The DMP-generated trajectories are also smooth splines that would fail to capture the imperfections that persist in the patient's behavior.

A Locally weighted regression has been used by [START_REF] Atkeson | Using locally weighted regression for robot learning[END_REF] to learn robot arm control by approximating the local model at each position. This approach while time-efficient for online learning does not provide enough generalization in the output model, it also may present a challenge for continuous monitoring applications since it is a memory-based approach and may present overflow challenges.

In [START_REF] Wisneski | Quantifying kinematics of purposeful movements to real, imagined, or absent functional objects: implications for modelling trajectories for robot-assisted adl tasks[END_REF] the study of the subject robot interactions by building a model of daily life activities using the standard zero jerk method was emphasized as being insufficient to capture all the settings that the subjects demonstrated. The authors also noted some significant differences in task execution and strategies which would be more apparent in pathological subjects. A serious note was put on the necessity to establish a more sensitive model to capture the curvatures that were demonstrated during the reach movements and which were reproduced for all the panel members.

To tackle these limitations, evolutive or iterative approaches have been proposed. For instance, in [START_REF] Chebotar | Combining model-based and model-free updates for trajectory-centric reinforcement learning[END_REF] the authors presented an evaluation of a hybrid of both modelbased and model-free reinforcement learning approaches and the respective algorithms used for training in the context of trajectory learning. A persistent challenge with these propositions is that they still assume fully known Markov Decision Processes besides to some knowledge of the system model.

For trajectory tracking, Iterative Learning Control [START_REF] Ahn | Iterative learning control: Brief survey and categorization[END_REF] for motor learning has been proposed as a theoretical framework to approach the evolution of the learning. These methods rely on the iterative approximation of the injected controls to ensure the convergence towards the final model. However, this can only match perfectly repetitive tasks such as in robotic manipulation. Besides, the assumption of knowledge of the patient model trajectory is unavailable in the context of rehabilitation training. The methodology assumes invariance in the dynamics governing the system evolution, whereas in rehabilitation, the dynamics are a hidden model to estimate. Moreover, the dynamics have the potential to evolve frequently considering the recovery process taking place.

In [START_REF] Jain | Learning trajectory preferences for manipulators via iterative improvement[END_REF] the authors presented an approach to learning an optimal task demonstrated by a user. The resulting model aims to control robot interactions through the reproduction of the movement.

In an assessment oriented application, rule-based techniques rely on a predefined set of rules to score and give feedback to the patient and therapists. As an example, a kinematic rule-based modeling technique was proposed by [START_REF] Zhao | Realtime motion assessment for rehabilitation exercises: Integration of kinematic modeling with fuzzy inference[END_REF] which defines an encoding of the exercise execution rules. Exercise executions are then compared against the established base truth to give feedback and assessment. Fuzzy inference is used to decide on quality assessment. This ap-proach lacks a dynamic update of the rule definition. The major benefit of the method is the ability to incorporate useful advice into the feedback given to the patient by specifying the execution error committed during the exercise.

In [START_REF] Capecci | A hidden semi-markov model based approach for rehabilitation exercise assessment[END_REF] the authors used an HMM to detect a set of features of patient's movements during the exercise using an RGB-D camera. They then provide an assessment based on a set of predefined rules that were concerted with a therapist. The definition of these rules is a non-trivial task. It involves a considerable amount of work to handwrite these rules and use them later for assessment. The defined rules are also specific to the task and are not generalizable.

In opting for a methodology to learn the patient model we were motivated by addressing some challenges:

• Continuous assessment and follow-up of the exercises;

• A tool to provide personalized feedback;

• A lightweight implementation for online deployment;

• The ability to generalize to similar rehabilitation systems.

We suppose that achieving these capabilities should permit a holistic parameterization of the patient's recovery, thus, providing current performance, rehabilitation rate, and potential recovery projections. We are adopting a statistical trajectory modeling approach to establish a means of assessing the performance of the patient. We seek to find the underlying model of the patient's hand movement during the exercise on an orthotic exoskeleton by capturing the motion of the end-effector using an HMM and imitation learning technique. Further analysis of the technique of imitation learning in the context of robotic trajectory modeling is referenced in [START_REF] Osa | An algorithmic perspective on imitation learning[END_REF].

Our main contribution is to provide an instrument to assess the longitudinal evolution of patients during the instrumented rehabilitation training sessions. Hence, a framework would be developed to model the trajectories recorded during patient movements to permit portability of this methodology to other exoskeleton devices, motion capture, or telerehabilitation systems. The framework relies on a dataestimated dynamic model. To measure the evolution of the rehabilitation new metrics will be proposed and tested statistically. We hypothesize that the resulting model will detect significant changes early on in the course of rehabilitation. We also hypothesize the significance of the findings will be pathology agonistic. Finally, we aim to evaluate the framework on an operational dataset to ensure these captured changes are detectable in the deployment environment.

Our starting point is the work that has been conducted by Coates et al. [START_REF] Coates | Learning for control from multiple demonstrations[END_REF] which proposes an HMM trajectory representation and an algorithm to infer the ideal trajectory of a given task. Our reasoning to opt for this methodology could be articulated in the following manner:

• Firstly, the HMMs representation permits us to capture the fine details of the trajectories considering that we use an upsampled chain to model the optimal execution. Meanwhile, in robotics, the trajectory modeling techniques often outcome a smooth trajectory. The fact of smoothing the resulting trajectories inadvertently contradicts our attempt to capture non-optimal and errored executions demonstrated by the patient;

• Secondly, the statistical representation would likely fit the stochasticity of the human Patient Controller and the Orthotic Exoskeleton system that we are studying, hence referred to as PCOE;

• Thirdly, the iterative procedure used permits constant update of the outcome model;

• Fourthly, the assumption of noisy demonstrations and the smoothing procedure would permit us to alleviate the rigor constraint on the data quality. This is important as quality is often imposed by the equipment; for instance, the sensors used on the Armeo Spring exoskeleton;

• Lastly, the portability of this methodology, learning state trajectories would permit the extrapolation towards a control task for any active-assist intervention using motorized exoskeletons as well as being applicable in the virtual environment systems.

To be able to capture the patient model, we first start by encoding the system dynamics. To this end, we implement a learning algorithm as in [START_REF] Abbeel | Learning vehicular dynamics, with application to modeling helicopters[END_REF] which permits us to approximate the state transition model. We present a newer iteration and cost function that can perform as well with better consistency. We show results presenting the ability of the model to keep satisfactory error rates while predicting state changes with both the original and the proposed implementations.

To model the patient ideal trajectory we implemented a variation of [START_REF] Coates | Learning for control from multiple demonstrations[END_REF] where we used the Expectation-Maximization (EM) algorithm to infer the underlying model. The resulting HMM model serves as the basis for defining a set of six new metrics that are studied statistically to investigate their properties.

In section 2 we start by introducing the design of the study and the population involved. The equipment used as well as the task studied is listed afterward. We then present the framework and detail its components. We end with detailing the definition of the proposed metrics, the statistical analysis used, and the data preprocessing methodology. In section 3 we listed results by first detailing the learning of the coefficient of the dynamic model. Secondly, we present the trajectory modeling results. Finally, the results of the statistical analysis are detailed for each of the tests that we carried out. In section 4, we discuss the results in light of potential use cases and utility in the domain of rehabilitation. We then conclude with the principal takeaways from this study and potential future research. The detailed population included in the case study. In a feasibility study we assess the proposed methodology and metrics by comparing the results of the two patients with those of the healthy user. We then investigate for significant differences between the measured metrics.

Methods

Study Design and Population

We conducted a longitudinal observational study on rehabilitation exercise data. The study was conducted on records of exoskeleton upper limb rehabilitation exercises retrieved from the Physical Rehabilitation Center at the University Hospital Center CHU of Tlemcen, Algeria. Both patients followed the standard therapy routines in parallel to the complementary exoskeleton sessions. The study aims to parameterize the recovery process through repeated measures in time. Outcome measures are defined based on the kinematically evaluated trajectories. The case study results are included as a feasibility study for the proposed framework. The subjects included are listed in Table 1. As a control, we compare the results of the patients to those of a healthy user.

Equipment and Task

The equipment used during the exercise session is an orthotic passive exoskeleton, the Armeo Spring as depicted in Fig. 1. The patients use this device to conduct vertical plane assessment exercise. The device logs the joint angles and the end effector positions during each session while permitting the interaction with a virtual environment through the resulting movement of a cursor therein. We conduct our learning phase on a dataset collected from the Physical Rehabilitation Center at the University Hospital Center in Tlemcen, Algeria where the device is presented to the patients as a complementary training exercise.

The exercise presented to the patients is to follow vertical plane movements to catch a target that appears on the screen by pointing at it with the cursor. Moving the orthotic exoskeleton fixed to the patient's upper limb in theplane moves the cursor on the screen.

We define an elementary trajectory (ET) to be the trajectory of the end effector to reach the next target catch position starting from the previous target position. The definition of ETs is depicted in Fig. 2.

The Framework Presentation

We propose to model the trajectories executed by the patient during the exercises. The outcome model serves as a tool to define a finer measure of the level of recovery. The algorithm is an EM aiming to infer an ideal model from a demonstration of the task presented at the input. The framework composition is presented in Fig. 3 main components:

and involves two

• The dynamic model which approximates the transition between states at times and + 1. The model follows the exhibited system properties by learning from the data. The system is the complex entity we note PCOE given the fact that the controller and actuator, i.e., the patient is changing dynamically during training. From a robotics standpoint, the only known parameters are the geometric configuration of the 6DoF robot. The patient acts as a controller that drives the muscles representing the actuators. The contractions of the muscles result in the final movement of the device;

• The model trajectory learning algorithm permits the inference of the hidden model trajectory of the patient for a given task from the demonstrated exercise.

We are thus attempting to encode the dynamics and the kinematic trajectory model of the patient during this elementary task.

State Transition Model Learning Algorithm

Piecewise function approximation Definition:

We consider that the PCOE system is observed over discrete periods ∈ Φ. Φ = { ∈ (1, … , )} where is the length of the particular task we are studying in time stamps. The system state is defined as = [ 1 , 2 ] ⊺ coordinates of the end effector in the reference frame attached to the tip of the gantry located at the device's shoulder sensor as appears on Fig. 1.

Assuming the evolution of the state is in the cartesian space and varies smoothly locally, we approach these variations using a third-order polynomial. The approximation can alternatively be seen as assuming the 3 derivative of the position, namely the jerk is constant which entails a smooth movement. The problem is formulated as optimizing the coefficient of the polynomial transition model to minimize the prediction error over an H-wide prediction window.

If we consider the ET's resampled frequency at approximately 51 , we can consider the transition between state and +1 is governed by the following equations:

̂ 1 ( + 1) = 1 3 1 ( ) + 1 2 1 ( ) + 1 1 ( ) + 1 ( ), (1)
̂ 2 ( + 1) = 2 3 2 ( ) + 2 2 2 ( ) + 2 2 ( ) + 2 ( ), (2)
(1) & (2) ⟹ ̂ +1 = 3 + 2 + + (3) with = [ 1 , 2 ] ⊺ , (4) = ([ 1 , 2 ]), and (5)
= ([ 1 , 2 ]), and (6)
= ([ 1 , 2 ]) (7)
We employ the optimization lagged error criterion algorithm [START_REF] Abbeel | Learning vehicular dynamics, with application to modeling helicopters[END_REF] to minimize the cost of successive predictions of state trajectories. The assumed local transition model is thus the approximation of the form:

̂ +1 = 3 + 2 + + (8)
Where , and are diagonal coefficient matrices in ℝ 2×2 as defined in Eq. (3). To run the optimization algorithm, we use exercise files with measurement logs of the position of the end effector during the predefined assessment task.

The Algorithm

We evaluate our coefficient matrices directly on the Cartesian coordinates namely = [ 1 , 2 ] ⊺ , ∈ ℝ 2 . We attempt to approximate the parameters of the cost value function by conducting estimates improvement using a black-box optimization technique on the exercise data. For this minimization, we define a cost value function as the cumulative incurred error during a prediction window of width H. The prediction errors are calculated starting from each timestamp t and sliding through the entire dataset. Predictions are calculated using the model in Eq. ( 8).

Algorithm 1 Compass : State Transition Model Learning

Initialize the 0 , 0 , 0 matrices by minimizing the one step prediction error:

OSPE = +1 --0 3 -0 2 - 0 repeat for = 0, … , -1 do Run ̂ +1 = ̂ 3 + ̂ 2 + ̂ + ̂ end for for = 0, … , -1 do for ℎ = 0, … , -1 do for = 0, … , ℎ -1 do = + ̂ 3 + | + ̂ 2 + | + ̂ + | end for = + ‖ +ℎ --‖ 2 2 end for end for ( ̂ , ̂ , ̂ ) = arg min , , +1 = (1 -) + ̂ +1 = (1 -) + ̂ +1 = (1 -) + ̂ until ‖ +1 -‖ + ‖ +1 -‖ + ‖ +1 -‖ <
In practice we run the algorithms on empirically determined hyperparameters, a learning rate set at 0.3 and a stopping criteria defined by a tolerance threshold of 10 -5 resulted in satisfactory accuracy at 68 iterations. The learning rate is set such that we change the estimated parameters slowly. By imposing this constraint, we ensure our learned coefficients are not majorly influenced given that the training data-set contains sharp turns and closed loops. The optimization was conducted using a Powell solver [START_REF] Powell | An efficient method for finding the minimum of a function of several variables without calculating derivatives[END_REF]. This optimization problem is defined as a multidimensional illconditioned, due to the cost function that is nonlinear in the parameters. The approximation that we are taking following the referenced work is that the coefficient matrices used in evaluating the sliding sub-sum in Algorithm 1 does not depend on previous estimates, which permits us to have a tractable problem formulation. The function is defined as the integral squared norm of the error incurred on successive predictions with complexity ( . . 2 ), where is the number of calls of the minimization solver, is the model length and is the window width. The Powell method does not depend on the differentiation of the cost function. This method makes comparatively less cost function calls during optimization iterations that make it a robust and efficient choice.

Patient Model Learning Algorithm

Hidden Markov Model Definition:

We consider each ET as a Hidden Markov Model where state represent a stochastic process taking a countable set of values ∈ [-, + ], = 1, 2 where represents the maximum range of motion for the ℎ coordinate. This value is commonly fixed for the assessment exercise to permit comparability between patients. The stochasticity of the state is our means to encode for the variability of the PCOE's states [START_REF] Bartolucci | Latent Markov models for longitudinal data[END_REF][START_REF] Ching | Markov chains. Models, algorithms and applications[END_REF].

Following [START_REF] Coates | Learning for control from multiple demonstrations[END_REF], we note as the hidden variable representing the ideal state at time instant . The hidden state vector is defined as = { , ∈ ℝ 2 , ∈ Φ }.

A naive transition model is assumed to be governed by the equation:

+1 ∼  ( 1 ( ), ) (9)
The function 1 represents the state transition model characterizing the PCOE system learned with Algorithm 1. The model noise is a zero-mean Gaussian distribution parameterized by covariance .

Each ℎ ET we observe represents an emission of the hidden vector and is noted

= { , = 1, … , , ∈ ℝ 2 }
, where the measurement state is evaluated through the observation model :

+1 ∼  (ℎ( ), ) (10)
The length of the HMM is fixed to twice the median of the observations' lengths. The resulting model is obtained at a higher frequency. The higher frequency HMM permits mainly to accommodate the acceleration effects that can be noted as increased spacing or "jumps" in the time series of each coordinate. This effect is notable considering the fixed sampling rate of the device. The naive predictor in Eq [START_REF] Brezany | Towards precision brain disorder rehabilitation[END_REF] does not natively accommodate this change and we can benefit from additional steps to compensate for the error.

The function ℎ represents the observation model and is assumed to follow:

+1 ∼  ( , ) Γ ( ), (11)
Where represents the measurement noise covariance and Γ is an indicator function for ∈ Γ , where

Γ = { , = 0 , … , -1 }.
The time alignment follows the probability law:

( +1 | , ) = , 3 ∑ =1 = 1 (12) with = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if +1 -= 1 2 if +1 -= 2 3 if +1 -= 3
The relationship of the hidden state vector with the observations vector is illustrated in the Fig. 4. Following the trellis diagram, the dynamic model returns, at timestamp , the naive prediction of the using the equation [START_REF] Brezany | Towards precision brain disorder rehabilitation[END_REF]. On the other hand, the observation model ℎ returns, at timestamp + 1, the last aligned smoothed hidden state with an added Gaussian noise defined by the covariance matrix .

To increase the accuracy of the model we augment the dynamic model of the HMM to be:

+1 ∼  ( 1 ( ) + , ) (13)
while is defined as: = ( -1 ( -1 )) and 0 ≡ 0. The extended vector is defined as follows:

= if Γ ( ) = 1 0 if ℎ
is a dampening coefficient determined empirically 1 . The bias permits to correct for the prediction errors of the model evaluated on the observations . We re-inject the bias into the model to compensate for the "observed" error in the model while predicting future states.

To estimate the parameters of these models namely, the distributions of the HMM state vector , the error covariance matrices , and the alignment distribution an optimization problem is formulated as:

, , , ( , , ; , , ) (14)
An Expectation-Maximization Algorithm 2 is used to estimate these parameters. Firstly, an extended Kalman smoother is used to infer the . Note that while is a learned parameter the problem cannot be properly defined without fixing it since the EKS iterations depend on deterministic alignment. We use the proposed uniform multinomial as an initialization value for the probability distribution . We also assume that the state's components vary independently. Therefore we defined the covariance matrices as diagonal initialized at , = ( ), , ∈ ℝ 2×2 . The smoothed distributions at the end of the Kalman iterations are used to update the covariance matrices, following the procedure:

̄ = ̄ +1| -1 -( ̄ | -1 ) (15)
1 The was set to 2∕3 in our experiment.

=  ( ̄ | -1 ) (16) = Σ | ⊺ Σ -1 +1| (17) = Σ +1| -1 -Σ +1| -1 ⊺ ⊺ - Σ +1| -1 (18) = 1 -1 ∑ =0 ̄ ̄ ⊺ + Σ | -1 ⊺ + (19)
With ̄ as the center of the Gaussian at position in the HMM chain.  represents the Jacobian matrix of the function 2 defined in Eq. ( 13).

= -ℎ( ̄ | -1 ) (20) = 1 -1 ∑ =0 ⊺ + ( ̄ | -1 )Σ | -1 ( ̄ | -1 ) ⊺ (21)
 represents the Jacobian matrix of the ℎ function 3 defined in Eq. [START_REF] Capecci | A hidden semi-markov model based approach for rehabilitation exercise assessment[END_REF].

We then proceed to use the estimated model to infer such that, given the distribution in Eq. ( 12) and fixing the covariance to the previously updated matrices, we have the following:

( , , , | , ) (22)
( , ; , , , )

( ; , ) ( ) (23)
[ ( , | ) + ( )] (24)
for each observation vector we determine the vector that best aligns the observed to the hidden. This process is performed using the subsequent formula:

0, = ( 0 ; 0 , ) ( 0 = ) ( 26
) , = ( ; , = ). ( 27
) max ′ [ ( ; -1 = ′ ). -1, ′ ] with ′ ∈ { -3, -2, -1} (28)
and

′ ∈ [2 -3, 2 + 3] (29)

The Algorithm

The implementation follows the suggestion in [START_REF] Coates | Learning for control from multiple demonstrations[END_REF]. The pseudo-code is listed in Algorithm 2 and we provide implementation specifications in what follows.

We start the iteration assuming a uniform multinomial distribution for values and an = 10 -5 for the covariance matrices. We start by running the standard Extended Kalman Filter using observations as measurements of the hidden. After reaching step we run backward pass that permits the estimation of the smoothed distributions' parameters,  ( | -1 , Σ | -1 ) [START_REF] Fiesler | Handbook of neural computation[END_REF]. We execute the Dynamic Time Warping algorithm [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] by maximizing the likelihood 2 The Jacobian matrix of the measurement model is evaluated to:  = ([(1), (2)]), (1) = 3

1 2 1 + 2 1 1 + 1 , (2) = 3 2 2 2 + 2 2 2 +
2 3 The Jacobian matrix of the observation model is evaluated to:  = (2) an identity matrix in ℝ 2×2

Metric

Symbol Definition Equation Hypothesis

Probability

The total probability defined as integral probability at each state of the hidden chain model.

( ) = ∫ -inf  ( ∶ ( ) , ( ) ) ∶ Φ → Γ Increase Total Variance
The sum of the covariance matrices along the hidden chain. = ∑ , ∈ Φ Decrease M Distance  M distance is defined as the sum of the squared root of the standardized differences. 26)- [START_REF] Romero | Convergence of the expectation-maximization algorithm through discrete-time lyapunov stability theory[END_REF]. The inferred solution is the time alignment vector that maximizes the likelihood of the observation given the alignment and the updated distribution. We update to the vector with the higher likelihood values and uses the Kullback-Leiber divergence as a termination test for convergence with a threshold = 10 -2 . The convergence can be safely assumed although not optimality, given the stability of the problem statement of maximizing the likelihood function which acts as a Lyapunov function for the EM iterations [START_REF] Romero | Convergence of the expectation-maximization algorithm through discrete-time lyapunov stability theory[END_REF].

 = ∑ √ ( ) , = -( ) , ∈ (0,

Assessment Metrics

Once the HMM model is evaluated for each ET during the session, we use it to define a set of metrics. The detailed definitions and formulation of the metrics are listed in Table 2.

Statistical Analysis

Statistical testing was applied to evaluate the evidence from the results of the modeling phase. First, the metrics did not prove normally distributed as a result of the Shapiro Wilks test. Consequently, we employed non-parametric tests instead of their parametric counterparts. The median and the Interquartile Range (IQR) were used as descriptive statistics.

The statistical significance thresholds is set to .05.

Evolution property

Firstly, we study the utility of each metric as an evolution metric. Hence, its ability to measure the evolution of the recovery process in a significant way. To this end, we used a within-subject dependent test for metrics samples at the beginning and the end of the treatment noted pre and post respectively. The test used is the two-tailed Wilcoxon Signed Rank test which is equivalent to the one-sample ttest. It assumes that the samples are dependent and with similar shapes to their distributions. The dependent variable is assumed to be continuous. Results of p-values are reported to eleven degrees of freedom = 11. The effect size was defined as = ∕ where W is the Wilcoxon Signed Rank statistic and n is the maximum sum of ranks. The statistic is the sum of the positive ranks and is reported as zero whenever there is none. A zero effect size follows a zero statistic by definition. Considering the fixed variability of the values of each metric and the continuous scale, the assumptions of the Wilcoxon Signed Ranks test are satisfied.

Discriminative property

Besides, we study the characteristic of a metric as a discriminative measure which refers to its ability to differentiate between different subjects at different stages of recovery. To this end, we used a between-subjects sample test for significance comparisons. The Kruskal Wallis test is used to check the differences between three or more samples. The same assumptions about the shapes of the distributions apply. Results of p-values are reported to = 2. Mutual differences between subjects were done via a post-hoc analysis. An independent sample two-tailed Mann Whitney Ranks Sum test was used to assert the mutual differences. Results of p-values are reported to = 11 and a Bonferroni rule of ( ∕3)% corrected for significance threshold. The effect size was defined as = ∕ 1 2 where U is the Mann Whitney test statistic representing the score for the second sample and 1 and 2 are the sample sizes.

Predictive property

Afterward, we assess the possibility of using the metric as predictive measure, representing the possibility of detecting trends of recovery status as a correlation relationship. To this end, we evaluated the Spearman's . The coefficient represents the strength of the monotonic relationship between the number of ETs conducted, as a form of training repetition, and the given metric. Results of p-values were reported to = 21.

Strength of evidence

To measure the likelihood of the alternative hypothesis of the non-parametric test we employ the Bayes factor bound (BFB) metric which defines the maximum odds of the alternative being true given the p-value of the test. This measure helps us ascertain the robustness of the conclusion and provide further evidence for the sake of reproducibility. It is defined as follows:

= -1∕ ln . All statistical tests were calculated using custom Python scripts using implementations of Scipy v1.3.2.

Data Preparation

The software has been developed using Python v3.7 with basic Anaconda modules. Data have been normalized using the Min-Max normalization to have the same scale for the position readings.

The elementary trajectories after segmentation were selected based on the successful completion of the task. We then re-sampled the results to the same rate to have a similar length per ET which is the median of the observations' lengths in a given session. This is necessary for the alignment phase where we seek to find the perfect fit with the ideal model which has twice the length. The resulting data are represented in Fig. 5.

The resulting trajectories are then translated to the origin and rotated such that each ends on the x-axis. This will help us unify the directions for all the trajectories and provide a baseline for comparison. An example of transformed trajectories is showed in Fig. 6 and thus permitting us to transform the learning into encoding absolute displacement instead of absolute position.

Results

In the following, we detail the results of the state transition model learning phase for both implementations comparatively and against the original work. Afterward, we detail the patient task model inference results. The statistical analyses are then presented while emphasizing the notable findings.

The Dynamic Model Learning Results

We start by retrieving an exercise file for each patient for the task defined in Fig. 2. We train the learning Algorithm 1 on the entire normalized exercise dataset (approximately 4000 timestamps) with different combinations of local models. We then chose the best three models to detail. Firstly, for the predictions are performed based The data are also unit normalized to ensure the same scale for all observations. ETs are then translated to the origin and rotated to define the movement as absolute displacements in the same direction.

on the transitions defined in Eq. ( 3) while for the a second-order variation of Eq.( 3) was used as the prediction equation. Finally, for the predictions were made using similar transition functions to the one defined in [START_REF] Abbeel | Learning vehicular dynamics, with application to modeling helicopters[END_REF], Algorithm 1 iterations were then updated accordingly to define the appropriate cost functions and rerun the algorithm.

We then proceeded to validate using the held-out method on an unseen exercise file to evaluate the accuracy of the predictions using error mean and standard deviations for the 20 timestamps starting at each time step in the data.

According to the definition in Eq. ( 12), the most possible number for blind predictions, i.e. without corrections by measurements , during EKS iterations is three timestamps. Accordingly, we detail the performance on the first three steps predictions for the trained models. In terms of execution time, the dynamic model is learned in 15 min up to 1 hour depending on the imposed tolerance. Fig. 7 shows that the defined as local cubic spline centered in achieves better accuracy results. The which implemented a parabolic model resulted in the second-best accuracy, while the original did not present any further accuracy improvements.

A detailed evolution of the mean error rate for the first three predictions is shown in Fig. 7, where we can note that both and achieve comparable accuracy levels while reaching at the third step with slightly better accuracy and consistency for the later with mean values of 1.53 × 10 -3 (2.25 × 10 -3 ) and 1.13 × 10 -3 (1.85 × 10 -3 ) respectively. We note that achieves lower standard deviation overall which entails better prediction accuracy with better consistency, justifying the proposition.

The Task Model Learning Results

We then proceed to establish a task model using the modified dynamic model defined in Eq. ( 13) augmented with a bias term .

In Fig. 8a we note the beginning of the iterations where the inferred model is too errored and the alignment is not properly maintained with a higher accumulation of error along the trajectory. The iterative process gives us a proper estimation of the noise covariances and better alignment sequences which are observed in Fig. 9 where alignment is achieved by the end of the iterations. The maximum likelihood alignment vector is chosen, this is paralleled with a decrease in the errors accumulated in the model.

A closer examination of the aligned demonstrations and the resulting model is showed in zoomed Fig. 8a where we see the beginning of the iterations with higher errors at the observations . The final iteration in Fig. 8b shows us how the converged better fits the demonstrations and the model all while keeping lower error values compared to early iterations. The 95% confidence interval around the smoothed trajectories does not show any particular accumulation of error, neither does the blind predictions in the trajectory. The smoothing procedure helps towards better and consistent modeling and thus prevents convergence to biased models at the end of learning iterations. The bias vector shifts in position with the same amplitude in coherence with its definition and has a maximum amplitude of 4% of the amplitude of the normalized measurement. Generation of a single report file takes 1 minute 33 seconds which is fairly reasonable for end of session report. The very nature of the metrics definition involving measurements from a full episode (i.e. an ET), limits to an extent a fully online use of the modeling framework without further modifications and experimentations.

The statistical analysis results

Modeling the patient's behavior as a means of extracting the underlying quality and performance of the rehabilitation exercise was presented in different forms and many studies proposed methodologies to proceed to that end.

Our objective in this study was to characterize the evolution of the patient during the rehabilitation exercises through the study of the evolution of the model of their trajectories. The modeling approach listed above evaluates a Bayesian model of each trajectory execution and thus permits us to define a set of new measures that we studied throughout the training.

The properties tested, namely the evolutive and predictive properties of the metrics, are for testing the evolution of the recovery over time and this is not the case of a healthy user whom we used as a comparative ideal for attained recovery. As such, the tests for both properties did not include results for subject C1.

Wilcoxon Ranks Sum Test

Observing the results of Table 3 shows that increase in median values for metrics , ,   and the decrease in the median of  ,  were consistent between the two subjects.

values were of the same order of magnitude for the same metric between pre and post samples and between samples of the two subjects.

To assess the hypothesis that the metrics measurements evaluated based on ETs from pre and post-treatment samples evolve significantly, and are thus capable of detecting differences between the performance levels of the patients, we used the Wilcoxon Signed Rank test. Differences between values of the entropy  and the total likelihood  were significant for subject S1. An odds ratio of 27:1 in favor of the alternative reinforced this result. Despite having nonsignificant statistics, the other metrics reported a moderate effect size of .2.

For subject S2 all metrics showed statistically significant differences between values at pre and post-treatment assessments. Significant differences had a 27:1 odds ratio in favor of the alternative hypothesis. The difference between the values of alignment likelihood  was insignificant, with a moderate effect size of .32. The results suggest that the metrics did not reach significant results for subject S2 simply because the process of recovery is still underway for that particular patient. This can be further justified if we examine the magnitude of differences between the median values pre and post-treatment for subject S1 that are twice the magnitude of the measure for subject S2. A further argument for this statement is the evaluation of the effect size that shows a small effect size for the metrics despite not reaching the statistically significant threshold. The results for the alignment likelihood were consistently insignificant, a result that can be related to the high values of which suggests the high dispersion in

.026 27

Table 3 A comparison between descriptive statistics for pre and post-treatment data shows similar trends between subjects. A Wilcoxon Ranks Sum Tests for pre and post-treatment differences in metrics' values detected many significant results. Nonsignificance is reported for both subjects for likelihood alignment  , this is due to the high value of the reported for all its samples. High odds ratio always favors the alternative and low effect sizes are still reported for the non-significant differences.

values. This dispersion makes the detection of a clear difference unfeasible and thus results in insignificant statistics.

Kruskal Wallis test

Evaluation of the between patients' differences can help distinguish notable differences between their performances. We studied the differences as observed through the Kruskal Wallis test and reported its results in Table 4. The test shows that the differences between the metrics for pre and post-treatment samples between the subjects were significant. This was reinforced with extreme odds reported in favor of the alternative.

Mann Whitney Ranks Sum test

Further analysis of these results using mutual tests of Mann Whitney are presented in Table 5. Comparing results for subject S2 to healthy control C1 gives significant results at the Bonferroni corrected threshold for all metrics except for  , , and  1 all being insignificant for pre-treatment samples comparisons. Despite the important p values, moderate to large effect sizes were reported .74 , .44 and .73 respectively.

Similarly comparing subject S1 to the healthy control showed significant results for all metrics except the alignment likelihood  . Nevertheless, the metric presented a large effect size of .5.

The results show that the proposed metrics can differentiate between the patients at the onset and the end of the treatment. Besides, it proves that the metrics are capable of characterizing the differences between the subjects at all stages of recovery. A result that can be very promising as it can be linked to the differences in treatment response that each patient will exhibit and that are otherwise undetectable. Whilst some metrics were reported to have insignificant statistics, the same conclusions may still be pronounced considering the reported high to moderate effect size.

Spearman's rho correlation

Further examining the evolution of each metric between the start and the end of the training can be seen in Table 6

Table 4

Investigating the differences between the three subjects at the start and end of treatment was conducted using the Kruskal Wallis test. The significant results are shown for all metrics and are reported with high odd ratios as evidence for the alternative. The results show that the metrics are significantly different between at least two of the subjects.

which detailed the correlation coefficients. These coefficients relay the strength of the monotonic relationship between the metric and the exercise repetitions, represented by the number of ETs accomplished.

Examining the results we can note that for the metrics proposed the coefficients reported for patient S1 are large effect sizes relating to the clear trend between the values of the metric at the beginning and the end of the training. These results are in line with the findings of the Wilcoxon tests. This also provides further details to the nature of the incremental relationship between the metric and the evolution of training repetitions.

The examination of the signs of the coefficient shows that similar trends are present for all the metrics and for both subjects despite not reaching significance for patient S2.

In addition to examining the trends of the metrics, the underlying structure of the model is represented by the different parameters that we estimate through the optimization techniques of the trajectory modeling Algorithm 2.

The observations that we can note on the value of the model parameters are a significant decrease in the total probability of the model and a decrease in the total variance.

The model noise variance is also increasing with medium effect size. We are thus able to detect changes at the level of the model that are characteristic of the changes in the structure of the trajectories that the patient exhibits throughout his training. This is a very compelling result as the intuition of the proposition of the model is the ability to capture the changes in the executions of the patient at different stages of recovery.

Detecting a clear trend on the parameters of the model indicates the ability of the proposed framework to encode for those differences in these parameters through the optimization performed using Algorithm 2 to estimate their values.

To assess the strength in the monotonic evolution of the metric between pre and post phases we use Spearman's rank. Results are reported in Table 6. For subject S2 we note three metrics  2 ,  ,  that significantly increased with moderate to large effect sizes and with moderate odds for the alternative.  showed a significant decreasing trend with large effect size and low odds ratio. A single model parameter  1 showed a significant increasing trend with moderate effect size and low odds ratio.

For subject S1, metrics ,  ,  showed a significant increasing trend with large effect sizes and extremely likely odds in favor of the alternative. Metric  showed an increasing tendency with large effect sizes and high odds as well.

Model parameters all decreased significantly at the exception of  1 that showed instead an increasing tendency. Results were associated with large effect size and extreme to likely odds in favor of the alternative.

Responsivity of the metrics

To further support this claim we carried a longitudinal analysis of the evolution of the metric values through the course of training for subject S1. We plotted the results of the p-value of the Wilcoxon test between the metric values at the beginning of the training and the next session completed by the patient that resulted in Fig. 10.

From these results, we can say that the metrics proposed can characterize differences in the status of recovery as early as after three sessions. These results can be highly beneficial in applications where the need for close follow up is strictly imposed, for instance in the telerehabilitation applications. On the other hand, the discriminative property is also present although the fluctuations may still leave several open questions to research as we can see in Fig. 11. The measure of the entropy of the model presented the most interesting results Table 5 Mann Whitney Test results assert that the metrics defined are capable of discriminating against the different subjects both at the onset of training as well as at the end. This can be seen in the significant results reported for comparisons between the three subjects. Large effect sizes and odd ratios also accompanied the results testifying the strength of the evidence reported.

Patients

between the metric as it is capable of significantly comparing the differences between the samples from different subjects regardless of the stage of recovery. The  presents another interesting behavior where, by the end of the training, it can detect no differences between the healthy user and the recovered patient. In the meanwhile, it still showed differences present between the recovered and the non recovered patient. At the beginning of the training though it shows differences between all three patients. Examining the results of Fig. 10, we can note that for all metrics we have a decreasing tendency for the p-values towards the significance threshold. The metrics on the exception of the alignment likelihood are detecting significant differences between the samples as soon as the third session since the start of the rehabilitation training.

Closely following the evolution of the p-values for the Mann Whitney test as depicted in Fig. 11 demonstrates that the entropy measure can distinguish subjects throughout the course of the training. The likelihood alignment does not have significant differences between the subject and the healthy user, while it presents some significant results at the onset of the rehabilitation between the two patients. The  measures are fluctuating throughout the training. The results by the end of treatment show that it is unable to detect significant differences between the patient and the healthy user by the end of the training. At the same phase differences with the other patient are significant which denote the recovery has progressed significantly for the subject S1 compared to S2. The  distance also presents fluctuations and is switching between significant and nonsignificant results

Patients

Table 6

The evolution of the metrics is tested using Spearman's rank. The test quantifies the existence and the strength of the monotonic variation of the metric with additional training. The significance of results was not consistent between the two subjects, whereas, coherent signs of correlations were reported. The model parameter  1 showed a significant increase for subjects and with moderate effect sizes asserting that the model adapts to the changes in the demonstrated dynamics of the movement while encoding for the trajectories.

throughout the training. The same remarks can be observed on the plots of the total likelihood values which are significant at the start and the end while fluctuating throughout the training.

Discussion

The utility

The rehabilitation exercises can be more effectively and accurately assessed using our proposed framework and the defined metrics. The new resulting metrics present strong evolutive, discriminative, and predictive properties that are key to assessment adoption for clinical validation and practice.

We modeled the trajectories of the end effector during rehabilitation exercises using an Armeo Spring exoskeleton using an EM procedure. The resulting model permitted the represent statistical significance thresholds. The results show that except the  all the metrics are capable of detecting significant differences between the performance for the post-stroke subject after one or two additional sessions. The results are moderately to strongly in favor of the alternative considering the BFB associated with the p values below 10 -3 .

definition of several metrics. The new metrics proposed to serve as a finer assessment measure of the patient's recovery. The new methodology permits detecting changes after 3 sessions with significant results for a number of the proposed metrics.

A major concern for any assessment technique is its sensitivity to the measured construct i.e. the recovery level. By showing a finer responsiveness property, the proposed framework promises to enable researchers to establish a much finer objectively evaluated kinematic based scale that correlates with the recovery process.

A strong discriminative nature for many of the resulting metrics would permit a proper classification of the patients into performance groups. By capitalizing on the correlation properties of the metrics, a narrow horizon future prediction of the outcome of the recovery might be feasible.

These findings promise to extend the kinematic approach that has been commonly used by researchers to evaluate the recovery of the upper limb using exoskeleton devices specifically. Some studies have stretched similar approaches further by proposing training to therapists using such models to ensure safe execution and reproducibility of correct exercises [START_REF] Salas | Using simulation-based training to improve patient safety: what does it take?[END_REF]. Another interesting use case has been offered in the context of modeling for diagnosis of the Parkinson [START_REF] Zhang | Fuzzy logic based assessment on the adaptive level of rehabilitation exergames for the elderly[END_REF].

A challenge that researchers often face is the inability to detect significance in recovery as measured by the traditional kinematic properties as we have shown in our clinimetric study on the same dataset [START_REF] Meziani | Evaluation clinimétrique des mesures cinématiques de rééducation du membre supérieur[END_REF]. Similar issues can be due to the small sample size and the heterogeneity of the population among others. The detection of longitudinal evolution by providing a finer assessment scale would enable researchers to study more precisely the morphological or neurological processes associated with the recovery while represent statistical significance thresholds. The discriminative power of each metric is viewed in terms of the evolution of the p-value of the between-subject Mann Whitney test. The results show that the entropy measure reports significant differences consistently. The p-value fluctuates near the significance threshold for the other metrics between the start and end of the rehabilitation.  metric seems to reach insignificance of differences between subject S1 and healthy control, suggesting a good metric for end of rehabilitation application.

enabling the possibility of parameterizing many subtle differences. Thus, enabling the possibility to accurately associate the phenomenon happening at the different levels and draw conclusions about causality or the underlying mechanisms of recovery.

In considering the underlying experimental unit as defined by the elementary trajectory, potential portability is open to any other application where the analysis of trajectories is of interest, for classification purposes for example. Meanwhile, exploration of other trajectories capturing devices might prove useful especially as there have been different input devices reported in research addressing the kinematic rehabilitation assessment.

Comparative studies that investigate the different approaches to physical rehabilitation interventions can benefit from the use of quantitative measures that are capable of detecting the small effects more reliably.

In future research, we would aim to conduct the study of the clinimetric properties of the metrics that we derived. This study will provide us with the properties and characteristics of each metric and its best use cases. Clinimetric characterization of the metric should also help move the metrics upfront towards a clinical adoption by asserting its validity.

The shortcomings

A case study of three subjects aimed to explore the feasibility of the tools proposed is far from providing enough evidence to establish the proposed framework. Nevertheless, the random nature of the framework and the heterogeneity of the data used can stand arguments in favor of the results provided.

The updated transition function while presenting relatively better accuracy has a higher iteration execution time. In the case of online deployment for a step by step prediction the higher accuracy model is encouraged since this training is to be done . As for any learning technique, the model derived can only be as good as the data used to train it and thus a good quality sensory reading from the exoskeleton is our baseline for all the models built and we are bound to the accuracy that it provides while detailing results of this study.

The data derived represent exactly the exoskeleton motions and does not necessarily represent the actual arm joints movements. Although there is not a direct relationship between the arm and exoskeleton angular displacement, the resultant movement still mirror some of the degrees of freedom of the arm joints and the end effector is submitted to the full control of the human controller. Arguably, any noticeable discrepancy in the recorded movements of the exoskeleton would have a physiological source at the functional or cognitive level and we are thus capturing a meaningful event.

The re-sampling of the data into similar lengths results in unwanted data changes and thus providing a more independ-ent way to update the model would be necessary to guarantee reliability without the need to introduce artificial measures into the data.

The normalization of the data while contributing to the stability of the algorithm might present the challenge when testing on newer data where coordinates maximum values are not known beforehand and thus might affect the final inferred model.

Conclusions

In conclusion, to address the accurate assessment of rehabilitation training from kinematic data we have developed a framework and listed a set of new evaluation metrics. The learning technique based on the Expectation-Maximization algorithm provides a means of modeling the trajectories of the end effector during the exercise. The resulting models are then used to derive a set of metrics that are used to study the recovery progress.

The results in this paper showed that the metrics derived by our methodology can provide finer measures over the progress of the recovery, demonstrating high evolutive, discriminative, and predictive characteristics.

These results while constituting a feasibility study for the framework, state also the effectiveness and sensitivity of the derived metrics in parameterizing the recovery process.

These findings can be highly beneficial in closely monitoring the rehabilitation status for patients, by providing minute follow-up while adapting to performance changes.

In future work, we aim to establish a validation study of the preliminary results that we presented. This would be accomplished by studying the results from a broader panel of subjects.

Figure 1 :

 1 Figure 1: The Armeo Spring device, a passive orthotic exoskeleton for upper limb rehabilitation training. The reference frame attached to the shoulder sensor defines all the positional measurements recorded in the dataset used for our study. The image is courtesy of Hocoma Inc. Switzerland.

Figure 2 :

 2 Figure2: A description of the vertical assessment exercise task where patients are required to point at targets appearing randomly on the screen by moving a cursor using the Armeo device. An elementary trajectory is thus the trajectory followed by the end effector to achieve this elementary task.

Figure 3 :

 3 Figure 3: The overall framework permitting the modeling of the patient behavior from exercise demonstrations by learning from data.

Figure 4 :

 4 Figure 4: The trellis diagram representing the hidden Markov model with the aligned observations vector. The alignment exemplifies the probabilities defined in Eq.(12).

Figure 5 :

 5 Figure 5: The training data was a sub-sample of the original exercise recording where we selected for the successful catch trajectories among the trajectories recorded represented in dotted lines against the whole session trajectory.

Figure 6 :

 6 Figure 6: The selected ETs are then resampled to the same rate to have the same number of timestamps per observation.The data are also unit normalized to ensure the same scale for all observations. ETs are then translated to the origin and rotated to define the movement as absolute displacements in the same direction.

Figure 7 :

 7 Figure 7: The broken bars plot represents the mean squared error evolution through the first three prediction steps. The minimum edge of the rectangle coincides with the mean error value of the model, while the length represents the variation of the measure as defined by the standard deviation . The graphs show that in addition to better accuracy the cubic model had lower standard deviations.

Figure 8 :

 8 Figure8: An example iterations over a single ET where we note that the algorithm starts with the standard alignment and error rates. The algorithm then reiterates over the observation and estimates a better alignment vector. The new alignment can be noted on the zoomed plots where observations represented in dots shifted positions between the start and end of iterations. The error on the resulting smoothed trajectories becomes less apparent as an effect of better matching of the observations to the hidden model.

Figure 9 :

 9 Figure9: The evolution of the total likelihood as estimated by the alignment EM iteration, the two extremities of each box represents the min and max values for all ETs. The values evolution between different iteration steps is depicted. We can note the clear increase in likelihood values for the observations given the alignment after the first iteration.

Figure 10 :

 10 Figure 10: The evolution of the p-values for the statistical Wilcoxon test for baseline metric measurements and those after each additional session of training for patient S1.represent statistical significance thresholds. The results show that except the  all the metrics are capable of detecting significant differences between the performance for the post-stroke subject after one or two additional sessions. The results are moderately to strongly in favor of the alternative considering the BFB associated with the p values below 10 -3 .

Figure 11 :

 11 Figure11: Plots of the log p-value of the Mann Whitney tests for mutual differences between subjects C1, S1, and S2. represent statistical significance thresholds. The discriminative power of each metric is viewed in terms of the evolution of the p-value of the between-subject Mann Whitney test. The results show that the entropy measure reports significant differences consistently. The p-value fluctuates near the significance threshold for the other metrics between the start and end of the rehabilitation.  metric seems to reach insignificance of differences between subject S1 and healthy control, suggesting a good metric for end of rehabilitation application.

Table 1

 1

	Patient Age Gender Affected Arm Pathology
	C1	27	Male	Right	Control
	S1	54	Female	Right	Post stroke
	S2	9	Female	Right	Infantile hemiplegia

Table 2

 2 The Metrics defined based on the HMM model of the ET. The metrics are measures of standardized distance and of probability.

	… , )	Decrease

Acknowledgment

The authors would like to thank Mrs. Meryem Karaouzene (PT) of the Rehabilitation Center, University Hospital Center of Tlemcen, Algeria for helping in providing access to the data record, patients' history, and subsequent explanations. This work has been supported by the PHC Tassili program agreement 19MDU210 and the DGRSDT Algeria.

Conflict of interest

None.

Graphical Abstract

Towards Adaptive and Finer Rehabilitation Assessment: a Learning Framework for Kinematic Evaluation of Upper Limb Rehabilitation on an Armeo Spring Exoskeleton

Yeser Meziani, Yann Morère, Amine Hadj-Abdelkader, Mohammed Benmansour, Guy Bourhis

Highlights

Towards Adaptive and Finer Rehabilitation Assessment: a Learning Framework for Kinematic Evaluation of Upper Limb Rehabilitation on an Armeo Spring Exoskeleton

Yeser Meziani, Yann Morère, Amine Hadj-Abdelkader, Mohammed Benmansour, Guy Bourhis • An individualized modeling approach of the trajectories of patients using the orthotic exoskeleton.

• Incremental model update to closely follow-up the recovery evolution.

• Intra-patient assessment for quick diagnosis of performance drops or training's negative effects.

• Introduce six new metrics to assess rehabilitation exercises.