In order to reliably evaluate the safety and stability of running train during operation, the distribution of maximum acceleration response of the running train under random track irregularity excitation is analyzed, based on the theory of multi-body dynamics and the theory of probability statistics. To this end, the train-track coupling model is established using the multi-body dynamics software Simpack. As the internal source excitation of the train-track coupling model, the track irregularity is simulated through trigonometric series method. In addition, the Monte-Carlo method is used to obtain multiple samples of the train’s acceleration response under random track irregularity. Subsequently, the maximum value of each sample (time history of train’s acceleration responses) is obtained and treated as a random variable for statistical analysis. The distribution of the train’s maximum acceleration response is fitted by the Gaussian Mixture Model, in which the fitting parameters are obtained through expected maximum algorithm and maximum likelihood estimation. Meanwhile, the reasonable number of samples to derive appropriate distribution model is also discussed. The results show that the Gaussian Mixture Model can accurately model the distribution of the train’s maximum acceleration response. Additionally, it is found that the variation of the train’s maximum acceleration response increases with higher vehicle speed.