Many real-world multivariable systems need to be modeled to capture the interconnected behavior of their physical variables and to understand how uncertainty in actuators and sensors affects the system dynamics. In system identification, some estimation algorithms are formulated as multivariate data problems by assuming symmetric noise distributions, yielding deterministic system models. Nevertheless, modern multivariable systems must incorporate the uncertainty behavior as a part of the system model structure, taking advantage of asymmetric distributions to model the uncertainty. This paper addresses the uncertainty modeling and identification of a class of multivariable linear dynamic systems, adopting a Stochastic Embedding approach. We consider a nominal system model and a Gaussian mixture distributed error-model driven by an exogenous input signal. The error-model parameters are treated as latent variables and a Maximum Likelihood algorithm that functions by marginalizing the latent variables is obtained. An Expectation-Maximization algorithm that jointly uses the measurements from multiple independent experiments is developed, yielding closed-form expressions for the Gaussian mixture estimators and the noise variance. Numerical simulations demonstrate that our approach yields accurate estimates of both the multivariable nominal system model parameters and the noise variance, even when the error-model non-Gaussian distribution does not correspond to a Gaussian mixture model.