Assessing seismic site response with absolute certainty remains challenging due to inherent soil variabilities. Soil variabilities primarily stem from natural geologic processes and lie beyond human control. Soil characteristics, including shear wave velocity, shear modulus, unit weight, and plasticity, exhibit inherent randomness and variability. These factors significantly influence seismic site responses, making it crucial to account for parameter uncertainty in seismic behavior characterization. Precisely characterizing this uncertainty is essential for reliable seismic hazard assessment and the design of earthquake-resistant structures. However, this task faces challenges due to the computational expense associated with many model realizations. In this study, a computationally efficient and user-friendly model for estimating peak ground acceleration (PGA), displacement (PGD), and quantifying uncertainty was developed. The methodology comprises two main steps: (1) 2D equivalent linear seismic site response analyses were simulated using MIDAS GTS NX commercial software. These simulations incorporated randomly generated properties of clay soil, including maximum shear modulus (G0), unit weight (γ), and plasticity index (PI). (2) Leveraging data from the site response analyses, a Bayesian regression model was developed using the R programming language. The accuracy and reliability of the developed model were validated using a new data set, and the results closely aligned with finite element method (FEM) outcomes. By accounting for soil inherent variabilities, the model effectively characterizes the uncertainty of PGA and PGD using mean and coefficient of variation (CoV). Remarkably, the Bayesian approach yielded CoV of response parameters up to 2.33%, a substantial 94.37% relative difference compared to the FEM. Notably, this improvement in uncertainty was achieved while maintaining computational efficiency.