This research demonstrates the results of an investigation into the California bearing ratio (CBR) of granular soils from Qassim region, Saudi Arabia, using multilinear regression (MLR), pure quadratic (PQ) models, and gene expression programming (GEP) methods utilized to develop mathematical models for estimating the CBR based on basic soil index properties. In this study, samples were collected from different borrowing pits in the Qassim area. Forty-three samples of soil were taken and transferred to a laboratory for examination. Seven multilinear regressions and seven PQ models were investigated, while four GEP models were made. The selection of each model variable depends on soil indices, grouping into grain size distribution, Atterberg limits, and compaction parameters. The results of this analysis showed that the PQ model had a higher accuracy [coefficient of determination (R2) = 0.89, root mean square error (RMSE) = 16.006, uncertainty (U95) = 16.17, and reliability = 57%] than the multilinear regression model, which has a lower accuracy model [R2 = 0.811, RMSE = 20.791, U95 = 15.569, and reliability = 51%]. The best GEP model yields [R2 = 0.776, RMSE = 22.552, U95 = 15.787, and reliability = 53%]. Furthermore, sensitivity analysis was conducted to distinguish the influences of different input variables on CBR; it was found that fines percentage (F200), maximum dry density (MDD), and optimum moisture content (OMC) are the most influential variables.