## 1 Introduction  
Gaussian process regression (GPR) is a nonlinear modeling technique that has recently received interest in various fields of engineering. It provides a non-parametric and probabilistic model fully described by its mean and covariance functions. One of the advantages of the GP approach over traditional methods is its ability to provide information for the system even in off-equilibrium regions, where limited data are available. Moreover, with this method, the hyperparameters can be directly learned from the data points [1]. Furthermore, the analytic expression of the model uncertainty provided by GP models as well as the robustness against overfitting are useful for model-based control strategies. GPs were first employed in the representation of the geostatistics field (also known as kriging) and then widespread in machine learning [2; 3], robotics, controls, and system identification frameworks [4; 5]. Models for various spatial and temporal problems, such as Brownian motion, Langevin processes, and Wiener processes are examples of GP applications [6]. A consistent number of publications can be found about the use of parametric prediction methods to obtain a mathematical model of a dynamic system from observed input--output data. Considering classical approaches, such as prediction error methods (PEM), statistical properties of prediction error are used to build an optimal model by selecting a proper model structure, as in [7]. Linear structures, such as black box structures, include the FIR model (finite impulse response model), OE model (output error model), ARX model (autoregressive with an external input model), ARMAX model (autoregressive moving average with external input model), and BJ model (Box--Jenkins model), but nonlinear black box modeling techniques are also available (see [8] for an overview). Some of the examples that can be found in the literature are neural networks, fuzzy logic-based models, and wavelet expansions. They is evidence of the influence of machine learning techniques in system identification fields. A survey of different techniques can be found in [9]. Furthermore, some research is focused on the non-convex optimization problem that provides the model parameters, and on methods to improve convergence to the optimal solution. In [10], the authors provide an overview on the use of evolutionary strategies and swarm algorithms in the field of system identification. In a more recent contribution [11],a novel fractional hierarchical gradient descent algorithm was used to find the NARX model parameters. The use of the GP model in nonlinear system identification is documented in [12], where Gaussian processes are presented as alternatives to other black box structures with the advantages of having a small number of required training parameters, a facilitated structure determination, and a confidence estimation for the model. A complete discussion about system identification and control using GPs is also provided in [13]. The application of GP for system identification can be seen in [14], where this approach is used for online target tracking and smoothing, and in [15], where it allows the description of the target dynamics in which external disturbances are explicitly incorporated. Unfortunately, GP is subjected to the \"curse of dimensionality\", and it requires high computational power both for training and prediction, which scales with \(\mathcal{O}(n^{3})\) and with \(n\) the number of data points. As a consequence, classic GP is not suitable for large datasets or online computation. Several methods have been proposed to solve this problem, such as sparse GPs [16] and recursive GPs [17; 18]. The novelty of this research lies in the application of GPR in the system identification framework for modeling ground robots.