This study introduces a probabilistic framework for assessing active earth pressures in soils exhibiting spatial variability in friction angles and unit weight. These properties are modeled using random fields with log-normal distributions and spatial correlation lengths. Monte Carlo simulations (MCS) are integrated with finite element limit analysis (FELA) to evaluate the failure probability under different design safety factors. To improve computational efficiency and prediction accuracy, machine learning models, such as Multivariate Adaptive Regression Splines (MARS), are utilized to predict failure probabilities based on key spatial variability parameters. A two-phase optimization approach, combining Random Search and Adaptive Sampling, is employed to refine the hyperparameters of the machine learning model. Confidence intervals are incorporated to quantify prediction reliability, providing engineers with robust decision-making tools under uncertainty. Furthermore, adaptive finite element meshes are applied to capture irregular stochastic failure mechanisms, offering deeper insights into the impact of spatial variability. The study produces parametric results in the form of practical contour design charts, aiding engineers in optimizing safety margins while accounting for soil variability. By combining computational methods, machine learning, and uncertainty quantification, this research enhances geotechnical design practices, ensuring more reliable and cost-effective solutions.