Running trains on a long-span bridge under earthquake is an inevitable problem for the high-speed railway in seismic regions. This paper aims to investigate the seismically random dynamic behavior of train–cable stayed bridges interaction system using the probability density evolution method (PDEM). The system simultaneously involves the multi-effect of random multi-point earthquake excitations with traveling waves and random track irregularity. The multi-point earthquake excitations, composed of lateral and vertical positions, are generated by stochastic harmful functions (SHFs) and uniformly modulated into nonstationary random processes. The motion equation of such a system is established by coupling the 38-degree vehicles and bridges through the refined random wheel/rail contact relationship and accounting for the phase-lags of seismic travelling waves between pier excitations. PDEM is proven to be applicable to such time-dependent systems and is then used to transform the random excitations into a series of deterministic representative excitations with initial probability. By solving for the corresponding deterministic probability responses, various nonstationary random responses, including the time-dependent probability density evolution functions of random responses, mean values curves and standard deviations, can be obtained easily. A case study based on the cable-stayed bridge is then presented, and the results show the effectiveness and accuracy of the proposed method by comparison with Monte Carlo Simulation. Additionally, the influence of train speed and seismic intensity on the safety and reliability of trains is discussed.