The dataset analyzed in this article is comprised of the second version of the publicly released NinaPro (Non-invasive Adaptive Prosthetics) database (Atzori et al., 2012; Gijsberts et al., 2014). The NinaPro database contains sEMG recordings collected from 40 intact subjects while performing a large number of common hand/finger movements and grasping actions. Of particular interest is the dataset containing measurements of the sEMG and the corresponding applied forces at the fingertips for different finger and thumb movements.

Each subject had 12 electrodes attached to their forearm and upper arm in order to measure sEMG signals. The 12 electrodes were placed in such a way that some electrodes provided dense sampling while others targeted specific muscles. Eight electrodes (1–8) were equally spaced around the forearm at the height of the radio-humeral joint, one electrode each was placed on the finger extensor (9) and flexor (10) muscles, and on the biceps (11) and triceps (12), see Figure 1. Six DoF force measurements were considered: flexion of the five digits (little, ring, middle, index and thumb) along with abduction of the thumb. These force signals were collected by way of a Finger-Force Linear Sensor (FFLS) (Kõiva et al., 2012). Each subject was asked to produce a set of nine force patterns, given in Table 2, by pressing down with one or more digits of the dominant hand in response to an external stimulus. Each force pattern was repeated six times, and a rest period was enforced in between each movement and each repetition in order to prevent muscle fatigue. More information and detail regarding the setup and data collection can be found in Gijsberts et al. (2014), Nin (2014).

The data preprocessing steps employed in this article follow the steps proposed in Gijsberts et al. (2014), since a benchmark for force regression on the NinaPro data has been set. The data was first split into a training and testing set based on the repetition of movements; the second and fifth repetition for each movement were used for testing and the four remaining repetitions were used for training. All the data were then standardized to be zero mean and unit variance, using statistics calculated solely from the training set.

A common practice used in smoothing sEMG signals is to first segment the signal into overlapping windows, and features are then extracted from each window. The sEMG signals were sampled at 2,000 Hz, and following on from Gijsberts et al. (2014) windows of 400 ms (800 samples) with a sliding window increment of 10 ms (20 samples) were used. For computational feasibility the training set was subsampled by a factor of 10 (at regular intervals) resulting in approximately 3,000 data points for training. The subsample rate used here is slightly higher than that used by Gijsberts et al. (2014) for regression analysis but preliminary analysis indicated that a higher subsample rate achieved a better performance. The MoE model is trained using the extracted sEMG features as inputs to the modeling algorithm.

Preliminary results using multivariate MoE for force regression alluded to low-dimensional FE representations outperforming high dimensional FEs (Baldacchino et al., 2015). Thus commonly used low dimensional time-domain FEs for sEMG signals are considered in this work, see Table 3. All the FEs considered have dimensionality equal to the number of sEMG channels (in this case 12). All the FEs, except for FILT, were performed on the windowed data as outlined in the previous paragraph. For FILT FE, a zero-phase second-order low pass Butterworth filter with a cutoff frequency of 2 Hz was applied to the full wave rectified sEMG data, separately for the training and test datasets. In order to obtain a similar number of data points to the other FE datasets, this dataset was subsampled at regular intervals of 200 samples.

The only instance of an MoE model applied to myoprosthetic control appears in Hahne et al. (2014) for regression analysis of wrist movements. The model was limited to two experts and trained using iterative reweighted least squares separately for each DoF. In Section 2.2.1, a sophisticated MoE is defined which is capable of accommodating unlimited DoFs and achieving analytical update solutions. The MoE model structure used in this article follows the work of Ueda and Ghahramani (2002), and further developed in Baldacchino et al. (2016). This particular choice of MoE model was adopted for two reasons; tractable closed form update equations for the parameters are obtained, and it is easily extended to accommodate multivariate outputs. These two characteristics reduce the computational effort during the training phase resulting in fast training times. The model used in this article is thus referred to as the multivariate Bayesian MoE model. Since a Bayesian approach to training the model is employed, prior distributions for the random variables need to be specified and these are given in Section 2.2.2.

Given the multivariate MoE and priors distributions for the random variables described in the previous sections, a learning algorithm is needed to train the model. An approximate Bayesian framework is desired in order to find the posterior distribution of the parameters p(Θ,a|Y). The marginal likelihood p(Y) (which appears in the denominator of Bayes’ theorem) is analytically intractable due to a multi-dimensional integral over the parameter space. The choice of conjugate-exponential distributions, along with a latent variable model is elegantly accommodated by the variational Bayes expectation-maximization (VBEM) framework. In short, the VBEM algorithm forms a lower bound of the marginal likelihood which is maximized iteratively so as to obtain a tight bound. Detailed information regarding the VBEM algorithm can be found in Beal (2003) and Beal and Ghahramani (2003).

Since conjugate priors for the model parameters were used, the functional form of the variational distributions will thus be the same as the priors. The optimal variational distributions are noted below, and expressed as q*(⋅). The VBM-step update equations are given in (10)-(16) while the VBE equations for the latent variables Z are shown in equatios (17)–(19). The equations given here are similar to the ones detailed in Baldacchino et al. (2016), however, some equations have been modified to cater for the multivariate output force signals.

The joint variational posterior distribution of the gates’ mean and covariance is a Gaussian-Wishart distribution, given by (10)q∗(μi,Λi)=Nμi|mi,(βiΛi)−1W(Λi|Bi,νi), where (11)mi=β0m0+∑n=1N rnixnβi,βi=β0+Ni,Bi−1=B0−1+∑n=1N rnixnxn′+β0m0m0′−βimimi′νi=ν0+Ni,Ni= ∑n=1N rni,rni=E[zni].

The joint variational posterior distribution of the expert functions’ mean and covariance is a matrix-Normal-Wishart distribution having the following form (12)q∗(Wi,χi)=MNWi|Ŵi,χi−1,LiW(χi|λi,Qi), where (13)Ŵi=Li[X1N×1]′ViYLi=[X1N×1]′Vi[X1N×1]+ϒi−1Vi=diag(E[z1i],…,E[zni])λi=λ0+NiQi−1=Q0−1+Y′ViY−Ŵi′(X′ViX+ϒi)Ŵi.

The term ϒ_i = ℰ[A_i] is defined in equation (16). The variational distribution for the ARD parameters is (14)q∗(ai,j)=Ga(ai,j|ci,di,j), where (15)ci=c0+0.5dy,di,j=d0+0.5ξi,jξi,j=dy(Li)j,j+λiw^i,jQiw^′i,j, where (L_i)_j,j is the jth diagonal element of L_i, and w^i,j∈R1×dy is the jth row of Ŵ_i. Using the statistic of a mean from a Gamma distribution, then (16)ϒi=E[Ai]=diagcidi,1,…,cidi,dx+1.

The VBE-step consists of updating the variational distribution of Z and the relevant equations are listed below. The variational distribution for the latent indicator variables follows a multinomial distribution, such that (17)lnq∗(Z)=∑n=1N ∑i=1M znilnrniandrni=γni∑l=1M γnl, (18)γni=πiexp12lnΛ^i+lnχ^i−ϖni−ξni.

The terms required in equation (18) are, (19)lnΛ˜i=E[ln|Λi|]=∑j=1dxψνi+1−j2+dxln2+ln|Bi|,lnχ˜i=E[ln|χi|]=∑j=1dyψλi+1−j2+dyln2+ln|Qi|, ϖni=Tr(E[Λi]E[(xn−μi)(xn−μi)′])=νi(xn−mi)Bi(xn−mi)′+βi−1dx,ξni=Tr(E[χi]E[(yn−[xn1]Wi)(yn−[xn1]Wi)′])=λi(yn−[xn1]Ŵi)Qi(yn−[xn1]Ŵi)′+dy[xn1]′Li[xn1]. where ψ(⋅) is the digamma function.

The update equation for the mixing coefficients is obtained by maximizing the variational lower bound with respect to the mixing coefficients. The update equation is given by (20)πi=1N∑n=1N rni.

The equation given above was used to determine the number of mixtures in GMMs (Corduneanu and Bishop, 2001), and it was later extended to a MoE framework in Baldacchino et al. (2016). Using (20), any surplus experts will have their π_i driven to zero, effectively removing them from the model. This equation needs to be interleaved into the iterative procedure since (18) depends on π. Following on from previous authors the mixing coefficients are updated after a pass of the variational posterior distribution update equations. The equations presented in this section are summarized in Algorithm 1.

The posterior predictive distribution allows one to perform predictions of the output force signals to new unseen sEMG signals. This distribution also gives information regarding the variance of the predictions, and thus confidence bounds can be calculated. The posterior predictive distribution is expressed as p(yn′|xn′,D), where 𝒟 = [Y, X], and n′=N+1 is the new data point. For the mixture of experts described in this article, the predictive distribution is given by, (21)p(yn′|xn′,D)=∑i=1M Ωi,n′Tyn′[xn′1]Ŵi,κiΣi,κi where {Ωi,n′}i=1M take value 1 with probabilities gi,n′MAPi=1M, respectively. gi,n′MAP=gi(xn′,πi,θiMAPg) using (2) at the maximum a posteriori (MAP) estimates θMAPg=μMAP,ΛMAP obtained from the posterior distribution (10), and the final value for π. At any given time n′ only one {Ωi,n′}i=1M can be 1 (the rest are zero) corresponding to the gate with the largest probability. T([xn′1]Ŵi,κiΣi,κi) is a multivariate Student-t distribution with κ_i = λ_i − d^y + 1 degrees of freedom, mean [xn′1]Ŵi and scale matrix κ_iΣ_i, where Σi=Qi(1+[xn′1]′Li[xn′1])−1. The relevant statistics for prediction are (22)E[yn′]=∑i=1M Ωi,n′[xn′1]Ŵi,cov[yn′]=∑i=1M Ωi,n′Σi−1κi−2.

The multivariate VBEM MoE algorithm, given in Algorithm 1, was trained on the sEMG/force data from all 40 subjects, and all the 6 force signals were trained together. The gate divisions on the input space are common to all the outputs. Since 9 movements and a rest period between repetitions and movements are considered here, the algorithm is initialized with 10 experts. The algorithm was run 100 times with random initialization for γ and m₀ so as to avoid the local maxima problem. The model with the largest lower bound was chosen as being the best model to represent the data. The algorithm was run for all the FE representations listed in Table 3, and a comparison between them is considered in Section 3.

In this article, the authors apply SA techniques with the intention of assessing the importance of the number of electrodes to the performance of the modeling framework. Interested readers are referred to Saltelli et al. (2000, 2008) for an in-depth exploration of SA. Here, the regression-based global SA technique developed in Xu and Gertner (2008) is applied since this method accommodates correlated input signals. The contribution of an individual sEMG signal to the variance of the model force output, V_j, is decomposed into two parts: the uncorrelated contribution VU_j, which describes variations unique to a sEMG input x^j which are independent of the other sEMG inputs, and correlated contribution VC_j, which explains variations of an sEMG input x^j which are correlated with other sEMG signals (Xu and Gertner, 2008) (23)Vj=VUj+VCj.

This discrimination between uncorrelated and correlated contributions highlights, respectively, if the sEMG input itself dominates or if the correlated variations among the sEMG signals dominate. The three indices found in equation (23) can be calculated using the regression-based SA equations found in Xu and Gertner (2008). The final indices of interest are the first-order sensitivity indices given by the total S_j, uncorrelated SU_j and correlated SC_j contributions of sEMG input x^j such that (24)Sj=VjV,SUj=VUjV,SCj=VCjV, where V=1Nlhs∑n (yn−y¯)2 is the estimated variance of the force output, and y¯ is the mean of the force signal. A Latin hypercube sampling (LHS) procedure is used to generate N_lhs rank-correlated samples for X from which the output {y}n=1Nlhs can be generated using the final MoE model. The LHS strategy used here follows that given in Iman and Conover (1982) since it takes into account the correlations between the inputs when generating samples. Using the sensitivity indices given in equation (24), examination of spurious sEMG signals is possible since they will affect the force output only due to their strong correlations with other significant sEMG inputs (that is, they will have a very low SU_j compared to SC_j). The regression-based SA approach relies on an approximate linear relationship between the sEMG signals and the force. This is not a bad assumption for the data used in this article since in Gijsberts et al. (2014), an average R² = 60% was achieved for a linear model using RMS FE.

The performance of the final regression models is analyzed using the normalized root mean square error (NRMSE) measure, given by (25)NRMSE=(1∕N)∑n=1N (yn−y^n)2ymax−ymin, and the coefficient of determination (R²), given by (26)R2=1−∑n(yn−y^n)2∑n(yn−y¯)2. where y^n is the model predicted output at time instance n, and y¯ is the mean of y These two statistics were chosen because of their wide use by the bio-robotics and rehabilitation/prosthetics communities in assessing the quality of their models.

The average normalized root mean square error (NRMSE) for each force DoF across all subjects, is plotted for all the different FEs in Figure 3. All FEs perform similarly and low NRMSE values are recorded. The multivariate VBEM MoE algorithm performs consistently well across all subjects since the SDs associated with each DoF are small compared to their mean values (i.e., the coefficients of variation are overall less than 30%). Focusing on the FILT FE as being fairly typical, the best results are obtained for the thumb abduction force signal (3.6 ± 1.2%), while the worst result is for the thumb flexion force signal (4.96 ± 1.5%). Prediction of little, ring, middle and index flexion force signals have similar performance.

Figure 4 shows the average NRMSE for all subjects (averaged over the 6 DoF force activations) for each FE, including unit SD (given by the black vertical lines). The performance of the algorithm is fairly consistent across all subjects, with most subjects achieving a level of NRMSE below 5%. Subject 37 for MAV has the worst overall NRMSE of 9.61 ± 1.76% while Subject 36 performs the worst for the other FEs. The best results were achieved by Subject 8 for WL and Subject 33 for MAV having an NRMSE of 2.64 ± 0.46 and 2.64 ± 0.22%, respectively (Subject 33 has the worst NRMSE for RMS and FILT). The mean (over all DoF and all subjects) and corresponding unit SD for both NRMSE (given by red and black horizontal lines in Figure 4) and R² are reported in Table 4. Overall, the FILT FE has the best performance with the other FEs having very similar performance.

To investigate the performance of the models on each of the 9 force patterns, the average NRMSE per pattern is shown in Figure 5. Similar to that reported in Gijsberts et al. (2014), patterns involving the individual activation of the four fingers (F1-F4) are all characterized by slightly better performance overall (i.e., by NRMSE values around 3%). The remaining movements (F5–F9) have a slightly worse performance (between 3.5 and 5%), especially those movements involving flexion of the thumb (F6 and F9); this is consistent with the results shown in Figure 3, where thumb flexion shows the worst NRMSE result.

Classification can be seen as a natural by-product of the VBEM MoE algorithm presented in this article. The learning algorithm consists of soft competition among the experts, such that only one expert dominates in different regions of the input space. The dominant expert is probabilistically chosen by assessing the gates using (2).

Insight into the relation between the different force movements and the sEMG signals is obtained via principal component analysis (PCA). Using PCA on the sEMG signals, the first two principal components can be plotted against each other, with the different force patterns each taking on a unique color, as shown in Figure 7. In this plot, the columns indicate 4 and 10 (9 movements plus rest) force patterns, respectively, while the rows show the effect of increasing the number of subjects. Even when just one subject is considered, increasing the number of movements results in the force patterns having lots of overlap, hence making it harder to distinguish between the different movements. Similarly, increasing the number of subjects also results in a large overlap. This figure also highlights the high variability in the recorded sEMG signals between subjects, and it is expected that this variability will be reflected in large deviations in classification performance. High variability among subjects was also observed by Atzori et al. (2012) for different hand movements.

A common metric for assessing the quality of the classifier is the classification accuracy, which is the proportion of correctly classified instances to the total number of instances. The left-hand plot in Figure 8 shows the average classification accuracy for each force pattern where F0 is the rest position, for all the FEs. All FEs achieve similar performance, and multifinger movements have a lower accuracy; more information can be obtained by analyzing the confusion matrix given in the top plot of Figure 9.

The confusion matrix profiles the misclassification of movements since correct predictions would result in a distinct black diagonal (left to right), while any non-clear off-diagonal cells are indicative of misclassifications. The top plot in Figure 9 shows that the VBEM MoE framework is consistent in correctly classifying movements, since the diagonal is prominent. The most conspicuous misclassification is given by the non-clear first column indicating that all force patterns are sometimes mistaken as rest, that is, the absence of movement. This phenomenon was also observed by Kuzborskij et al. (2012), and they attribute this to: 1/the data contains rest-to-movement segments which causes ambiguity and these segments are technically neither rest nor non-rest, and 2/when windowing is used, some windows include sections of rest and non-rest samples and so these samples are again neither rest nor non-rest. For the FILT FE only point 1 is valid (since no windowing is used), and this is enough to confuse the algorithm into thinking that these regions of the data are associated with no movement. In Chan and Englehart (2003), it was reported that this type of error is treated as being acceptable since it is expected that due to mechanical inertia, a prosthesis would not be able to respond to transitory misclassifications. In the literature regarding finger movements classification, the rest position is either ignored or the sEMG data stream is segmented in such a way that regions of ambiguity are removed from the training set.

The right hand plot of Figure 8 and the bottom plot of Figure 9 show the results when only sections of the data set are considered when calculating classification accuracy, that is, each movement is divided into three equally sized segments and only labels from the center segment were retained (no retraining was performed). A higher accuracy is now attainable: single finger movements and multi-finger movements have an average accuracy of 87.89% and 62.39%, respectively.

Referring back to the confusion matrix plots (Figure 9) gives useful insight into why multi-finger movement achieves lower accuracy than single finger movement. Movement F8 (flexion of ring and middle fingers) has the lowest classification accuracy from all movements. The algorithm mostly misclassified F8 for the corresponding individual finger movements of F2 (flex ring) and F3 (flex middle). The same trend is observed with movements F7 (flexion of index and little fingers) and F9 (flexion of index finger and thumb).

The SA procedure outlined in Section 2.6 is applied to the dataset using the final MoE model obtained with the FILT FE (similar results were obtained using the other models, so only the results of one FE are shown here). All 12 sEMG signals are treated as individual inputs, and the SA results revealed that the contributions of all the sEMG input signals to the variance of the force outputs are mainly due to the correlations with other inputs, since all inputs had a low SU index compared to SC index.

The SA is performed again but this time the sEMG input signals are grouped according to their electrode location. Thus the groups comprise of electrodes 1–8 (forearm), 9–10 (finger extensor and flexor), and 10–11 (biceps and triceps)—refer to Figure 1. The three sensitivity indices, averaged over all 40 subjects, for all the output force signals for these groups of inputs are shown in Figure 10. The SU index for electrodes 1–8 is now significant (and comparable to the SC index) for most outputs. Hence, the sEMG electrodes located around the forearm are necessary for describing the relationship between force at the fingertips and the sEMG signals because this group of inputs’ individual contribution is important. The other 2 groups of inputs are almost entirely dominated by their correlated effect, and so their individual contribution is low.

The multivariate VBEM MoE algorithm was run again but this time only the 8 electrodes around the forearm were used as inputs to the model. The plot shown in Figure 11 shows the performance of the models on the reduced input set (solid lines) versus that of the full input set (dashed lines) for both regression and classification over the different force patterns. The top plot shows the quality of the models for the force regression given by the average NRMSE for the different force activations considered. A slight degradation in NRMSE is observed for the reduced input set (across all FEs), however, the models are still performing exceptionally well. The bottom plot demonstrates the quality of classifying the different finger movements including the rest position using the average classification accuracy. In this case, the reduced set and full set models have similar performance, with no clear winner.

The results and observations obtained for the force regression using the multivariate VBEM MoE are comparable to those reported in the literature, for example (Castellini and Kõiva, 2012; Gijsberts et al., 2014). Predictions of thumb movements are worse than for the other fingers due to no sEMG activity being recorded from the majority of the thumb muscles since these are located at the wrist (and hence are not usually available on the amputee’s stump) (Kõiva et al., 2012; Gijsberts et al., 2014).

In the benchmark set by Gijsberts et al. (2014), overall R² values of 91.74 and 88.93% for marginal discrete wavelet transform (mDWT) and RMS, respectively, for a exp-χ² kernel ridge regression were reported. The results of Table 4 indicate that the VBEM MoE algorithm achieves similar performance to that given in Gijsberts et al. (2014) but with the use of a much simpler FE representation for sEMG signals. The mDWT has a considerably higher feature dimensionality (36 inputs) than the FEs considered in this article (12 inputs). The use of a low dimensional FE allows for faster feature computation of the windowed sEMG signals compared to the mDWT in the preprocessing stage. Faster training times are also achieved when a smaller input dimension is used in any algorithm. It has also been shown that TD features performed better than frequency domain and WTs in real-time applications (Englehart and Hudgins, 2003).

The NRMSE values for each finger force pattern are comparable to those reported in the literature, such as (Castellini and Kõiva, 2012; Gijsberts et al., 2014). The R² values are not reported for the individual force patterns since the sum of error over each force pattern is not necessarily 0, and so the R² value can go negative although the model provides a good fit. Therefore, it is not possible to perform a direct comparison to the R² values reported in Gijsberts et al. (2014) since it would result in an erroneous interpretation. However, the NRMSE values reported here highlight what the authors in Gijsberts et al. (2014) and others have established: force signals from single finger movements are easier to predict than from multiple finger movements and thumb movements.

The misclassification of multiple finger movements reported in this article tends to be represented by the corresponding individual finger movements. This result seems to suggest that multiple finger movements can be accurately represented by their corresponding individual finger movements implying that the sEMG signals of the multiple finger movements can be decomposed into the individual sEMG signals. This outcome appears to give weight to the hypothesis presented in Castellini and Kõiva (2012) whereby the authors suggested that for example flexion of ring + middle would be statistically similar to flexion of little + ring + middle. To the best of the authors’ knowledge, this phenomenon of sEMG signals of some multiple finger movements being similar to the sEMG signals of the individual fingers has not been analyzed in the literature. Possible explanations for this could be that similar finger movements have not been tackled (some articles do not specify which multiple finger movements are being analyzed), and confusion matrices are not reported hence by-passing the opportunity to assess in detail how the classifier is misclassifying movements.

The classification accuracy reported in this article is similar to that found in the literature, refer to Table 1. The algorithm presented here is capable of classifying finger movements which is the result of a by-product from regression analysis—no extra computational expense is incurred. Another important point to note is that this is the first instance of simultaneous regression and classification of finger movements via sEMG signals using MoE models.

The initial SA investigation indicated that the contribution to the force variance was a result of correlated contributions of all the sEMG signals. This result seems to suggest that the VBEM MoE learning algorithm is not sensitive to the nominal placement of electrodes since the uncorrelated contribution of all the electrodes was low suggesting that individual contribution of the electrodes themselves are not important. This result supports the research performed in Hargrove et al. (2008), where the authors report that the pattern recognition framework is insensitive to nominal electrode placements, however, it is sensitive to electrode displacement during the training/testing phase. Hence, the same locations need to be used for training and testing in order for a classifier/regressor to have a high accuracy.

Good force regression and movement classification was achieved by the VBEM MoE model on a reduced input set which consists of a dense sampling of electrodes around the forearm. This setup in which no relevant muscles are targeted was considered by Castellini and Kõiva (2012) such that 9 electrodes were uniformly positioned around the forearm. The regression results obtained here are comparable to those reported by Castellini and Kõiva (2012). Reducing the number of electrodes reduces the overall cost in terms of both physical cost (less electrodes) and computational cost (training the model). With a training dataset of around 3,000 samples, an average training time of 114 s for the full input set and an average of 93 s for the reduced input set was achieved. Therefore, SA techniques provide a structured methodology for reducing the number of electrodes needed for clinical application.