This paper focuses on the problem of nonlinear dynamic response variability resulting from stochastic system properties and random loads. An efficient and accurate method, which can be employed to analyze the dynamic responses of random finite element systems with local nonlinearity, is presented in this paper. This method, dubbed as the partition expansion method, is based on the partitioned time integration algorithm in conjunction with the Neumann expansion technique within the framework of the Monte Carlo simulation. Two numerical examples involving structural and mechanical stochastic vibration problems are employed to illustrate the advantage of the proposed method with respect to accuracy and efficiency. By comparing the results obtained by the direct Monte Carlo simulation, the dynamic response statistics can be accurately determined using the proposed method with four order expansion while the computational efforts are significantly reduced. The comparison of computing time indicates that the proposed method is efficient and practical for analyzing the statistical quantities of stochastic dynamic systems with local nonlinearity.