Ground motion parameters prediction (peak ground acceleration, PGA and peak ground velocity, PGV) is of the essence in rescue efforts aftermath of earthquakes and seismic hazard analysis. As new developed approaches for predicting ground motion parameters, machine learning algorithms do have some advantages, but also have difficulties in estimating predictive uncertainties and interpreting machine learning models. In this study, we use the natural gradient boosting (NGBoost) algorithm to evaluate predictive uncertainties, and use the SHAP values to interpret trained machine learning models. Based on NGA-WEST2 database, we trained machine learning models which are suitable for predicting PGA and PGV in active tectonic regions. The correlation coefficients between the predicted PGA and PGV and observations in testing dataset reach up to 0.972 and 0.984, respectively. The trained machine learning models also provide reasonable probability distributions of predicted values. With the computed SHAP values, we figured out the influence of the input features (moment magnitude, M-W ; Joyner-Boore distance, R-jb; V-S over top 30 m, V-S30; rake angle, Rake; dip angle, Dip; depth to the top of fault, Z(TOR); and depth to V-S=2.5 km.s(-1), Z(2.5)) on the outputs of machine learning models. According to the SHAP values of input parameters, we find that the predicting mechanisms of trained machine learning models make sense in physics which illustrates the machine learning models are reasonable. In addition, SHAP values also revealed some facts which are ignored in previous studies: (1) The SHAP values of Z(TOR) in general are low when the depths of rupture planes are shallow (Z(TOR) < similar to 5 km), indicating that the ground motions from ruptures in the shallow part of crust may be controlled by velocity strengthening and are systematically weaker.The SHAP values of Z(TOR) decrease with Z(TOR), which indicate ground motions from ruptures in the shallow part of crust may also be affected by depth-varying geometrical attenuation; (2) When depths of ruptures are large (Z(TOR) > similar to 5 km), the SHAP values of Z(TOR) increase with Z(TOR), which indicate ground motions from ruptures in the deep part of crust may highly be impacted by depth-varying stress drops or quality factors (Q); (3) The variations of SHAP values of Z(2.5) are different for predictions of PGA and PGV when Z(2.5) are low (Z(2.5) < similar to 1 km), which indicate impacts of differences in resonance frequencies of sediments caused by variations of Z(2.5) on PGA and PGV are different, since the frequencies of velocity and acceleration are different.