As a relatively new means of transportation, the subway has become an important tool for the sustainable development of many cities. Being buried deep in soil under the weight of vital infrastructure, subway stations can be vulnerable to seismic excitations. Considering the high randomness of ground motions, it is important to research the failure probability and seismic performance of the subway station based on stochastic dynamic analysis. In this paper, a probability density evolution method (PDEM) coupled with a spectral representation random function is used to analyze the stochastic dynamic response and seismic probability of a subway station. First, according to the improved power spectral density model and the seismic design code of urban rail transit structures in China (GB 50909-2014), a set of nonstationary ground motions consistent with the code spectrum are obtained. Then, a great deal of deterministic dynamic calculations for Daikai subway station considering soil–structure interaction based on elastic–plastic methods are performed. In addition, the nonlinear stochastic response analysis and the dynamic probability analysis are obtained for the subway station by solving the PDEM equation. Finally, the probability density function (PDF) and cumulative distribution function (CDF) of the subway station under stochastic earthquake excitations are obtained based on three performance indices, including drift angle in the middle column, relative vertical displacement between floor and roof, and damage area ratio (DAR). The results show that the stochastic dynamic analysis and the probability density evolution method can analyze seismic response and evaluate seismic performance of subway stations effectively. The proposed method will serve as an effective tool for the seismic design of underground structures.