In order to convert ellipsoidal heights obtained by the Global Navigation Satellite System (GNSS) to orthometric heights, it is necessary to know the distance between the ellipsoidal and geoid surface, called the geoid undulation. The geoid undulation can be predicted using emerging mathematics tools and algorithms. The objective of this study was to develop a model for predicting the geoid undulation using Gaussian Process Regression (GPR), one of the soft machine learning algorithms having different covariance functions. This method was then compared with the radial basis function neural network (RBFNN), generalized regression neural network (GRNN), and the interpolation method of inverse distance to a power (IDP) with the power of 1, 2, 3, 4, and 5. First, 70 % of GNSS/leveling data (422 points) were used in the training phase. The remaining 185 points were used as testing data to check the effectiveness of the constructed model. In the GPR modeling, ten covariance functions (Materniso d= 1, 3, 5; Maternard d= 1, 3, 5; SEiso; SEard; RQiso; and RQard) were tested for prediction on this dataset. The GPR based on the Materniso (d = 1) covariance function model was introduced as an effective method for predicting geoid undulation and provided the best results (RMSE = 8.32 cm, MAE = 5.51 cm, R 2 = 0.98968) when compared with the other developed GPR models. In addition, the statistical findings showed that the accuracy of all the GPR models was also better in predicting geoid undulation than the RBFNN, GRNN, and IDP with the power of 1, 2, 3, 4, and 5.