In this paper, the probabilistic solutions of the multi-degree-of-freedom (MDOF) or large-scale stochastic dynamic systems with polynomial type of nonlinearity and excited by Gaussian white noise excitations are obtained and investigated with the subspace method proposed recently by the authors. The space of the state variables of large-scale nonlinear stochastic dynamic (NSD) system excited by white noises is separated into two subspaces. Both sides of the Fokker–Planck–Kolmogorov (FPK) equation corresponding to the NSD system is then integrated over one of the subspaces. The FPK equation for the joint probability density function of the state variables in another subspace is formulated. Therefore, the FPK equation in low dimensions is obtained from the original FPK equation in high dimensions and it makes the problem of obtaining the probabilistic solutions of large-scale NSD systems solvable with the exponential polynomial closure (EPC) method. A simple flexural beam on nonlinear elastic springs is analyzed with the subspace method to show the effectiveness of the subspace-EPC method in this case.