Nonnegative matrix factorization (NMF) is a powerful tool for parameter estimation applied in numerous robotics applications, such as path planning, motion trajectory prediction, and motion intention detection. In particular, NMF has been successfully used to extract simplified and organized movement primitives from myoelectric signal (MES) for robust control of multi-degree of freedom humanoid robots. However, MES is typically contaminated by complex noise sources. The system performance often degrades due to the simplified Gaussian assumption of the noise distribution in existing NMF methods. Furthermore, most existing NMF models are unable to automatically determine the rank of the latent matrices. To address these issues, this article presents a hybrid variational Bayesian Gaussian mixture and NMF (GMNMF) model with a finite Gaussian mixture model adopted to fit the mixed noise density function of MES. In addition, the automatic relevant determination criterion is applied to automatically infer the number of movement primitives. The coordinate descent update rules for the proposed model are formulated by mean-field variational Bayesian inference. We assess the model performance on five synthetic noise distribution functions and an experimental MES dataset to perform six wrist movements. The results demonstrate that GMNMF yields low error and high robustness in extracting the movement primitives over four competitive methods for robust cybernetic control.