## 6 Results and Discussion  
A consistent random seed was applied to divide the dataset into 70% for training and 30% for testing, ensuring reproducibility in the model evaluation. In engineering applications, it is often essential to maximize predictive accuracy while keeping the resulting expressions simple; however, these objectives can be at odds. To derive the final equations, each symbolic regression algorithm was trained through multiple iterations, progressively increasing the allowed complexity. In each case, a pool of candidate formulas was generated and assessed based on a trade-off between predictive performance and simplicity. To avoid overfitting and ensure model interpretability, complex or higher-order terms were penalized or excluded. The final expressions were chosen by comprehensively evaluating each candidate for both predictive performance and structural simplicity. The results of this process were the formulas included in Table 4.  
To ensure reproducibility, the main hyperparameters adopted for each model were included in Tables 11-14. Hyperparameters not listed were those set by default in each model and therefore not modified. As much as possible, all algorithms were restricted to the use of basic operators: addition, subtraction, multiplication, division, exponentiation, square root, and logarithm, excluding more complex ones. For the application of these equations, it is important to consider the input units, which are LL, PI, w, and #0.075 mm in % and \(\gamma_{\text{d}}\) in kN/m\({}^{3}\), yielding p\({}_{\text{s}}\) in kPa. Additionally, the ranges in which the algorithms were trained must be taken into account, as predictions outside these ranges may be unrealistic. The following are the ranges:  
* Pass #0.075 mm (%): \(\geq\)28.2 and \(\leq\)100.0. * LL (%): \(\geq\)27.7 and \(\leq\)85.1. * PI (%): \(\geq\)7.9 and \(\leq\)57.8. * w (%): \(\geq\)8.9 and \(\leq\)49.1. * \(\gamma_{\text{d}}\) (kN/m\({}^{3}\)): \(\geq\)13.0 and \(\leq\)19.2. * p\({}_{\text{s}}\) (kPa): \(\geq\)17.2 and \(\leq\)514.2. The parameters selected by each SR algorithm were different, except for PySR and MOSR, which selected the same parameters (all except w). This provides an advantage, as it enables the use of the equations depending on the available parameters. Natural water content was included in only one of the four equations (PGSR), in line with the findings of the statistical analysis conducted previously. It is normal not to observe a clear impact due to the undisturbed nature of the samples, as each sample has varying levels of water content and saturation. Figure 7 presents scatterplots illustrating the relationship between the actual values of p\({}_{\text{s}}\) (\(x\)-axis) and the predicted values of p\({}_{\text{s}}\) (\(y\)-axis) for the four SR algorithms considered. The analysis demonstrated that a significant portion of the data points lay in close proximity to the identity line and were mostly located within the error range of \(\pm\)15%, indicating a strong agreement between the predicted and actual values. This suggests that the models performed well in accurately predicting p\({}_{\text{s}}\), with minimal deviation from the true values. The residuals (defined as the difference between the real and predicted values) are also shown in Figure 6, revealing that PhySO, MOSR, and PGSR slightly \begin{table}
\begin{tabular}{c c c} \hline
**SR Model** & **Equation** \\ \hline PhySO & p\({}_{\text{s}}\)(kPa) = 0.070(-PI -0.150 LL) -(00.075 mm + PI + 0.180) & (1) \\ PySR & p\({}_{\text{s}}\)(kPa) = \(\frac{\text{PI}\cdot(0075 mm)+\text{PI}^{2}}{\gamma_{\text{d}}}\) - \(\pm\)147.092 & (2) \\ MOSR & p\({}_{\text{s}}\)(kPa) = \(\frac{\text{2}\cdot(0075 mm)+\text{PI}^{2}}{\gamma_{\text{d}}}\) & (3) \\ PGSR & p\({}_{\text{s}}\)(kPa) = \(\frac{\text{PI}\cdot(0075 mm)+\text{PI}^{2}}{\gamma_{\text{d}}}\) & (4) \\ \hline \end{tabular}
\end{table} 
Table 4: Formulas obtained for each SR algorithm implemented.