This paper presents a method to generate random samples of soil dynamic shear modulus and dynamic damping curves with full consideration of the correlations of critical soil dynamic parameters to investigate the influence of their uncertainties on the engineering site seismic response in the implementation of equivalent linear method. A one-dimensional (1D) equivalent linear site seismic response analysis program, which serves stochastic dynamic response analysis of the engineering site, has been developed in MATLAB. A 1D free field model for the typical layered engineering sites of site class II is established in this study. The target response spectra, which are defined with spectra of the outcrop corresponding to different earthquake return periods based on the acceleration response, are employed to generate artificial seismic records. These records are scaled down by ½ and referred to as input motions at the engineering bedrock for the 1D free field model. The numerical results show that the uncertainties of the critical soil dynamic parameters has significant influence on seismic response of engineering sites, which are highly related to the amplitude and the frequency aspects of the input ground motions as well as the fundamental periods of the engineering sites. The variations of the peak shear strain and the peak ground acceleration of the site, which reach 10% and 14%, respectively, increase with the amplitudes of the input ground motions. Besides, the variations of acceleration response spectra corresponding to the plateaus of the target response spectra and the fundamental periods of engineering sites exceeds 20% with the consideration of the uncertainties of soil dynamic parameters.