This paper studies the effect of VS and layer thickness uncertainty on ground response analysis (GRA), considering statistical randomization. Initially, the effect of VS uncertainty on GRA was studied, and then combined VS and layer thickness uncertainty were simulated. In the study of VS uncertainty, the parametric variations of log standard deviations (SDs) of VS ( \({\sigma}_{\rm{ln}}V_{S}\) ), interlayer correlation coefficient (CC), and shaking level of ground motions on GRA have been studied. In the study of combined uncertainty, the parametric studies of coefficient of variations (CoVs) of layer thickness uncertainty and intensity of input ground motions are simulated. In this regard, a MATLAB code for equivalent linear GRA in the frequency domain combined with VS and thickness randomization procedure has been developed. The simulation considers 300 randomized profiles for each case with a total number of simulations ~6000. The study clearly exhibits significant variations in GRA due to the uncertainty in VS. The GRA outcomes in terms of SDs of amplification spectra (AS), response spectra (RS), and peak ground acceleration (PGA) variations are found to increase with \({\sigma }_{\rm{ln}}V_{S}\) . Another source of uncertainty that significantly affects the GRA outcome is the intensity of input ground motion. The SDs of AS, RS, and PGA variations are observed to rise as the intensity of input motion increases. This is due to the pronounced soil nonlinearity exhibited during the strong ground shaking. Little higher SDs in GRA outcomes are observed in combined uncertainty cases than VS uncertainty alone. Although the differences are not that significant for the current considered base profile, this might vary from one site to another. Finally, the cumulative distribution functions (CDFs) of the surface PGA variations for different cases have been estimated for the site.