Impact Statement
The proposed uncertainty-incorporated approaches extend the application range of popular random forest and evolutionary polynomial regression and ensure their practicability for reliability-based design in various industries. The proposed approaches can meet diverse requirements due to their respective merits in strong fitting capability and formulable expressions and hold substantial potential in the field of computational science.
1. Introduction
The identification of mechanical indices of materials is vital in engineering design. Experiments in materials can reveal these indices as design parameters but take much effort and expense (He et al., Reference He, Zhang and Yin2024a). Alternatively, some readily measured physical characteristics are often employed to estimate these indices through certain transformation models (Phoon et al., Reference Phoon, Cao, Ji, Leung, Najjar, Shuku, Tang, Yin, Ikumasa and Ching2022). Numerous empirical relations and data-driven approaches were proposed to identify mechanical indices of materials in the past (Karimi et al., Reference Karimi, Rahmani and Jia2023; Rezania et al., Reference Rezania, Faramarzi and Javadi2011; Vázquez-Escobar et al., Reference Vázquez-Escobar, Hernández and Cárdenas-Montes2021). However, the majority of these methods neglected underlying uncertainties and merely generated point estimations (Li et al., Reference Li, Rahardjo, Satyanaga, Rangarajan and Lee2022; Nassr et al., Reference Nassr, Javadi and Faramarzi2018; Pei et al., Reference Pei, Qiu and Shen2023; Rezania et al., Reference Rezania, Faramarzi and Javadi2011; Wang et al., Reference Wang, Yin, Zhang and Jin2020; Zhang et al., Reference Zhang, Wu, Li, Wang and Samui2019, Reference Zhang, Yin, Jin and Chan2020, Reference Zhang, Yin and Jin2021b). The direct application of these methods without uncertainty quantification (UQ) may pose huge risks in engineering practice.
Uncertainty is an inevitable factor contributed by various sources in geotechnical engineering. In general, uncertainty can be categorized into epistemic and aleatoric uncertainties (Zhang et al., Reference Zhang, Jin and Yin2021a), which are also known as reducible and irreducible uncertainties (Olivier et al., Reference Olivier, Shields and Graham-Brady2021). The former arises in modeling the studied phenomena of interest and is dominated by a lack of knowledge, while the latter is induced by slightly intrinsic randomness and is hard to mitigate such as measurement error in laboratory tests. Thereafter, numerous research works focus on reducing prominently epistemic uncertainty if there is no specific demonstration. More specifically, epistemic uncertainty is primarily contributed by intrinsically spatial variability due to widely spaced materials (Ching and Phoon, Reference Ching and Phoon2019; Guan and Wang, Reference Guan and Wang2022; Lyu et al., Reference Lyu, Hu and Wang2023), transformation uncertainty in estimating design parameters (Phoon and Kulhawy, Reference Phoon and Kulhawy1999), and statistical uncertainty due to limited site investigation data in geotechnics (Ching and Phoon, Reference Ching and Phoon2017; Guan and Wang, Reference Guan and Wang2020). The current study focuses on the second source of epistemic uncertainty, that is, transformation uncertainty associated with data-driven models and aims to tailor traditionally deterministic approaches for broader applications.
To capture substantial uncertainty for reliability-based design, many scholars have been equipping traditionally deterministic approaches with UQ to improve their reliability and practicability. For example, Zhang et al. (Reference Zhang, Yin and Jin2021b) enhanced the reliability of neural networks (NNs) by incorporating Monte Carlo dropout (Gal and Ghahramani, Reference Gal and Ghahramani2016), which adopts random dropout of weights to generate mean and variance for predicting soil maximum dry density. Zhang et al. (Reference Zhang, Yin and Jin2022b) further integrated NN with variational inference (VI), which approximates the parametric posterior distribution of weights for interval estimates of soil compression index
$ {C}_c $
and creep index
$ {C}_{\unicode{x03B1}} $
. The bootstrapping technique has also been used to quantify uncertainty in ensemble learning models, such as NN (Olivier et al., Reference Olivier, Shields and Graham-Brady2021; Psaros et al., Reference Psaros, Meng, Zou, Guo and Karniadakis2023; Zhao et al., Reference Zhao, Tan, Lin, Yin and Yin2023) and support vector regression (Li et al., Reference Li, Yin and Liu2023). Apart from the preceding integrations, Wang and Zhao (Reference Wang and Zhao2017) proposed a novel Bayesian compressive sampling approach that weighs over prespecified basis functions to capture the uncertainty of the tip resistance profile. Ching et al. (Reference Ching, Phoon and Wu2022) also designed a hierarchical multivariate normal probability density functions (PDFs) to generate three-dimensional (3D) probabilistic predictions of
$ {C}_c $
.
Nevertheless, some popular approaches with unique structures, such as random forest (RF) and evolutionary polynomial regression (EPR), have not been equipped with UQ in numerous geotechnical applications (Nassr et al., Reference Nassr, Javadi and Faramarzi2018; Pei et al., Reference Pei, Qiu and Shen2023; Rezania et al., Reference Rezania, Faramarzi and Javadi2011; Zhang et al., Reference Zhang, Yin and Jin2021b). Through incorporating bootstrap aggregating (Breiman, Reference Breiman1996) and random subspace (Ho, Reference Ho1998), RF characterizes flexible tree structures built upon dichotomy and exhibits strong fitting capability (Karimi et al., Reference Karimi, Rahmani and Jia2023; Wang et al., Reference Wang, Yin, Zhang and Jin2020; Zhang et al., Reference Zhang, Yin, Jin, Chan and Gao2021d). In contrast, EPR works by combining symbolic and numerical regressions to search for a polynomial function, known as user-friendly and applicable expressions in vast scenarios (Ghorbani and Hasanzadehshooiili, Reference Ghorbani and Hasanzadehshooiili2018; Gomes et al., Reference Gomes, Gomes and Vargas2021). Based on these characteristics, RF has been widely used in soil indices prediction (Li et al., Reference Li, Rahardjo, Satyanaga, Rangarajan and Lee2022), 3D reconstruction of granular grains (Zhang et al., Reference Zhang, Yin and Chen2022a), landslide susceptibility mapping (Sun et al., Reference Sun, Xu, Wen and Wang2021), settlement estimation (Zhang et al., Reference Zhang, Wu, Li, Wang and Samui2019), and slope stability evaluation (Pei et al., Reference Pei, Qiu and Shen2023), while EPR has been applied in constitutive modeling (Nassr et al., Reference Nassr, Javadi and Faramarzi2018; Zhang et al., Reference Zhang, Yin, Jin and Liu2021e), liquefaction evaluation (Ghorbani and Eslami, Reference Ghorbani and Eslami2021; Rezania et al., Reference Rezania, Javadi and Giustolisi2010), settlement prediction (Shahin, Reference Shahin2016), and so forth. Nonetheless, absent UQ hampers their scenarios, and hence, the current study aims to tailor them for broader applications.
A key issue behind RF is that the output of each decision tree only relies on the average observation in a certain leaf and neglects the distribution information inside. The distribution information typically can be interpreted by PDF, cumulative distribution function or quantiles (Hao and Naiman, Reference Hao and Naiman2007). Owing to this neglect, uncertainty was not involved in the aforementioned studies of RF because it merely averages the generated conditional mean of each tree. Some efforts have been made to address this issue, but they only pay attention to classification tasks (Reis et al., Reference Reis, Baron and Shahaf2018; Vázquez-Escobar et al., Reference Vázquez-Escobar, Hernández and Cárdenas-Montes2021). Regarding this issue, quantile regression holds the potential to utilize the distribution information in each leaf because it can estimate arbitrary quantiles (Koenker and Bassett, Reference Koenker and Bassett1978) such as the quartile or median. Therefore, RF deserves to be integrated with quantile regression to enhance reliability for broader applications in scientific computing domains.
On the other hand, the coefficients of polynomials in EPR are scalars, simply fitted by the least-squares method. In this aspect, Jin and Yin (Reference Jin and Yin2020) discussed the potential of equipping EPR with uncertainty evaluation by replacing the least-squares method with a Bayesian approximation through the Markov chain Monte Carlo (MCMC) (Marasco et al., Reference Marasco, Marano and Cimellaro2022). Nonetheless, MCMC comes with several limitations: sampling directly from a target PDF becomes intractable when the prior and likelihood PDFs are quite different (Phoon et al., Reference Phoon, Cao, Ji, Leung, Najjar, Shuku, Tang, Yin, Ikumasa and Ching2022); it suffers from the curse of dimensionality and becomes inefficient if the number of target parameters is large (Lykkegaard et al., Reference Lykkegaard, Dodwell and Moxey2021); it tends to be stuck in local convergence for some tasks (Lee and Song, Reference Lee and Song2017). VI has recently become an effective alternative to the MCMC method and has already been incorporated into NN for probabilistic shear strength predictions of clays (Zhang et al., Reference Zhang, Yin and Jin2022b). The integration of EPR and VI is hence worth investigating for designing an efficient modeling approach with UQ.
To this end, this study proposed two uncertainty-incorporated baseline modeling approaches, that is, quantile-based RF and VI-based EPR (QRF and VIEPR) to enrich the current data-driven modeling framework and identify material indices. As a vital design parameter, the
$ {C}_c $
of clays was taken as an exemplar to develop models and verify the feasibility of the proposed approaches. Therefore, two models denoted by QRF and VIEPR with three input variables (i.e., initial void ratio e
0, liquid limit LL and plasticity index Ip) were developed and applied to predict
$ {C}_c $
values. Their performances were fairly compared in terms of accuracy and reliability and also compared with results of deterministic counterparts. The dependency of error and uncertainty on input variables was also investigated. The better-performing method was further applied to predict
$ {C}_{\unicode{x03B1}} $
for exploring its generalization ability.
2. Methodology
2.1. QRF
As a type of nonparametric algorithm, RF incorporates bagging (Breiman, Reference Breiman1996) and random subspace (Ho, Reference Ho1998) to generate trees and leaf nodes for making predictions. Specifically, bagging generates a bootstrap set by sampling with replacement for training each tree, while random features are considered for node split. In this context, each tree is an independent predictor, in which the output for a given dataset is determined by the average observation in a certain leaf. RF further averages outputs generated by each tree to control overfitting (Zhang et al., Reference Zhang, Yin and Jin2021b). However, such a strategy only considers the average observation in leaf nodes and neglects the distribution information inside, offering predictions without reliability.
Figure 1 illustrates RF using similar notations as Breiman (Reference Breiman2001), let
$ {\theta}_k $
be learned splits that determine the growth of the kth tree
$ {T}_k $
based on a bootstrap set
$ {D}_k $
. Each bootstrap set is sampled from original training datasets D
$ =\{({X}_i,{Y}_i)|{X}_i\in {\unicode{x211D}}^m,{Y}_i\in \unicode{x211D},i=1,2,\dots, n\} $
.
Figure 1. Schematic of RF.
As can be seen from Figure 1, dropping arbitrary input x from the top of
$ {T}_k $
will ultimately fall in a unique leaf
$ \lambda \left(x,{\theta}_k\right) $
. To calculate the output, observations in this leaf are assigned the same weight
$ w\left(x,{\theta}_k\right)=1/C\left[\lambda \left(x,{\theta}_k\right)\right] $
, where
$ C[\cdot ] $
counts the number of observations inside. Conversely, observations outside the leaf are assigned zero weight. Therefore, the output of a single tree
$ {T}_k $
can be calculated by weighting each observation
$ {Y}_i $
as follows:
$$ {\hat{Y}}^k=\sum \limits_{i=1}^n{w}_i\left(x,{\theta}_k\right){Y}_i $$
where
$ {w}_i\left(x,{\theta}_k\right) $
denotes the weight of ith observation in
$ {T}_k $
. When it comes to RF with K trees, the weight of each observation
$ {w}_i(x) $
becomes the average of weights over all trees as follows:
$$ {w}_i(x)=\frac{1}{K}\sum \limits_{k=1}^K{w}_i\left(x,{\theta}_k\right) $$
Consequently, the prediction of RF can be expressed by
$$ \hat{Y}=\sum \limits_{i=1}^n{w}_i(x){Y}_i $$
Note that
$ \hat{Y} $
is merely an average of the conditional mean generated by leaf nodes and neglects distribution information inside. To address this issue, the integration of RF with quantile regression (Koenker and Bassett, Reference Koenker and Bassett1978; Meinshausen and Ridgeway, Reference Meinshausen and Ridgeway2006), QRF is used to extract distribution information in leaf nodes for reliability evaluation as shown in Figure 2.
Figure 2. Schematic of QRF.
Different from RF that weights over each observation
$ {Y}_i $
in Equation (3), QRF replaces
$ {Y}_i $
with a conditional distribution
$ P\left({Y}_i\le y\right) $
to estimate the CDF of predictions by:
$$ \mathrm{CDF}(y)=\sum \limits_{i=1}^n{w}_i(x)P\left({Y}_i\le y\right) $$
where
$ P\left({Y}_i\le y\right) $
represents the probability of
$ {Y}_i\le y $
, which equals 1 if
$ {Y}_i\le y $
or 0 otherwise. As a result, a schematic of predicted CDF can be drawn in the right plot of Figure 2. By performing the inverse of Equation (4), arbitrarily specified quantile
$ {Q}_{\unicode{x03B1}} $
of predictions can be calculated as follows:
$$ {Q}_{\unicode{x03B1}}=\operatorname{inf}\left\{y:\mathrm{CDF}(y)\ge \unicode{x03B1} \right\} $$
where inf{
$ \cdot $
} denotes the infimum of a set and is applied to a set of y that allows
$ \mathrm{CDF}(y)\ge \unicode{x03B1} $
.
To sum up, QRF provides a nonparametric distribution of y as a function of x, instead of merely a conditional mean generated by RF. For QRF, the median is employed for point estimations due to its robustness toward outliers, similar to previous research works (Ching et al., Reference Ching, Phoon and Wu2022; Löfman and Korkiala-Tanttu, Reference Löfman and Korkiala-Tanttu2021). For interval estimate, the commonly used 95% credible interval (CI) can be formulated by quantiles
$ {Q}_{2.5\%}(x) $
and
$ {Q}_{97.5\%}(x) $
as the lower and upper bounds, respectively. It is noteworthy that QRF keeps the same structure as RF when absorbing distribution information over numerous trees. Compared with existing Bayesian methods (Ching et al., Reference Ching, Phoon, Li and Weng2019; Guan et al., Reference Guan, Wang and Phoon2024; Wang and Zhao, Reference Wang and Zhao2017; Yoshida et al., Reference Yoshida, Tomizawa and Otake2021; Zhang et al., Reference Zhang, Yin and Jin2022b), the prominent difference of QRF is its frequentist manner to generate interval estimates simply using votes inside trees without involving any prior information, rather than update prior distributions to posterior distributions. On the other hand, QRF only involves the simplest dichotomy for tree growth and avoids complex sampling. As a baseline algorithm, QRF is generic and can be applied in broad applications, for example, 3D reconstruction of granular grains (Zhang et al., Reference Zhang, Yin and Chen2022a) and landslide susceptibility mapping (Sun et al., Reference Sun, Xu, Wen and Wang2021), not limited to material indices identification. Given these promising characteristics, such an approach holds the potential of further fusing data with limited features to refine the generated CDF for interval estimates.
2.2. VIEPR
Inspired by numerous empirical polynomials used for estimating material indices (Nagaraj and Srinivasa Murthy, Reference Nagaraj and Srinivasa Murthy1986; Zhu and Yin, Reference Zhu and Yin2000), EPR has been proposed to automatically search for a polynomial function (Giustolisi and Savic, Reference Giustolisi and Savic2006). EPR works by combining symbolic and numerical regressions in two phases: structure identification and coefficient estimation. For structure identification, evolutionary algorithms are used to search for an exponent matrix
$ \boldsymbol{E}\in {\unicode{x211D}}^{nt\times m} $
such that m original input variables can be combined with interior exponents to form
$ nt $
new transformed variables XT as follows:
$$ {XT}_j={x}_1^{\boldsymbol{E}(j,1)}\cdot {x}_2^{\boldsymbol{E}(j,2)}\cdot {x}_i^{\boldsymbol{E}(j,i)}\cdots \cdot {x}_m^{\boldsymbol{E}(j,m)} $$
where
$ {x}_i $
= ith original input variable;
$ {XT}_j $
= jth transformed variable;
$ \boldsymbol{E}\left(j,m\right) $
= (
$ j,m $
)-entry of
$ \boldsymbol{E} $
; Next, coefficient estimation is implemented by the least-squares method to obtain a general EPR expression as follows:
$$ y=\sum_{j=1}^{nt}{w}_j\cdot {XT}_j+{w}_0 $$
where y = predicted output;
$ {w}_j $
= an adjustable coefficient of
$ {XT}_j $
;
$ {w}_0 $
= an additional coefficient as bias. Although EPR provides efficient search and explicit solutions, it only generates predictions without reliability.
To address this issue, all coefficients
$ \boldsymbol{w} $
in Equation (7) can be denoted by probability distributions, such as normal distributions instead of scalars. Suppose that they comply with a certain prior distribution
$ P $
(
$ \boldsymbol{w} $
), their posterior distribution then can be formulated by Bayes’ theorem as follows:
$$ P\left(\boldsymbol{w}|D\right)=\frac{P\left(D|\boldsymbol{w}\right)P\left(\boldsymbol{w}\right)}{P(D)} $$
where
$ P\left(D|\boldsymbol{w}\right) $
is the likelihood of observing D given
$ \boldsymbol{w} $
;
$ P(D) $
is model evidence, which can be regarded as a constant. From the view of Bayesian, Jin and Yin (Reference Jin and Yin2020) discussed the potential of equipping EPR with uncertainty by approximating the distribution of coefficients through the MCMC method (Al-Bittar and Soubra, Reference Al-Bittar and Soubra2014). However, solutions may be intractable due to large computational costs and convergence problems, especially for high-dimensional tasks. Thus, VI is introduced to address this issue.
As an approximation method, VI utilizes a parametric distribution
$ q $
to approximate target posterior
$ P\left(\boldsymbol{w}|D\right) $
. Specifically, VI learns variational parameters
$ \boldsymbol{\theta} $
in
$ q $
to make parametric probability density
$ q\left(\boldsymbol{w}|\boldsymbol{\theta} \right) $
and target posterior PDF
$ P\left(\boldsymbol{w}|D\right) $
as close as possible. This objective is mathematically achieved by minimizing the distance between two distributions, measured by Kullback–Leibler (KL) divergence (Kullback and Leibler, Reference Kullback and Leibler1951). Given two distributions p and q, the KL divergence between p and q can be formulated by
$$ K L(q\Vert p)={E}_q[\mathrm{ln}\hskip2pt q-\mathrm{ln}\hskip2pt p]=\int q(\boldsymbol{w})\mathrm{ln}\frac{q(\boldsymbol{w})}{p(\boldsymbol{w})}\mathrm{d}\boldsymbol{w} $$
where
$ q\left(\boldsymbol{w}\right) $
and
$ p\left(\boldsymbol{w}\right) $
denote the probability densities given certain
$ \boldsymbol{w} $
.
Therefore, the variational parameters
$ \boldsymbol{\theta} $
are found by minimizing the KL divergence between parameterized
$ q\left(\boldsymbol{w}|\boldsymbol{\theta} \right) $
and target posterior
$ P\left(\boldsymbol{w}|D\right) $
as follows:
$$ {\displaystyle \begin{array}{c}\hskip-9.839995em \hat{\boldsymbol{\theta}}=\mathrm{arg}\hskip2pt \underset{\boldsymbol{\theta}}{\min \hskip2pt } K L[ q(\boldsymbol{w}|\unicode{x03B8} )\Vert P(\boldsymbol{w}| D)]\hskip0.24em \\ {}\hskip-3.7em =\mathrm{arg}\hskip2pt \underset{\boldsymbol{\theta}}{\min \hskip2pt }\int q(\boldsymbol{w}|\boldsymbol{\theta} )\mathrm{ln}\frac{q(\boldsymbol{w}|\boldsymbol{\theta} )}{P(\boldsymbol{w})P(D|\boldsymbol{w})}\boldsymbol{d}\boldsymbol{w}\\ {}\hskip2.7em =\mathrm{arg}\hskip2pt \underset{\boldsymbol{\theta}}{\min \hskip2pt } K L[ q(\boldsymbol{w}|\unicode{x03B8} )\Vert P(\boldsymbol{w})]-{E}_{q(\boldsymbol{w}|\boldsymbol{\theta} )}[\mathrm{l}\mathrm{n}\hskip2pt P( D|\boldsymbol{w})]\end{array}} $$
Using MC sampling to evaluate the expectation term in Equation (10), the KL divergence is approximated and taken as the loss function as follows (Olivier et al., Reference Olivier, Shields and Graham-Brady2021; Zhang et al., Reference Zhang, Yin and Jin2022b):
$$ L\left(D,\boldsymbol{\theta} \right)\approx \frac{1}{s}\sum \limits_i^s\left(\ln q\left({\boldsymbol{w}}^{(i)}|\boldsymbol{\theta} \right)-\ln P\left({\boldsymbol{w}}^{(i)}\right)-\ln P\left(D|{\boldsymbol{w}}^{(i)}\right)\right) $$
where
$ {\boldsymbol{w}}^{(i)} $
denotes the ith MC sample drawn from the variational posterior
$ q\left({\boldsymbol{w}}^{(i)}|\boldsymbol{\theta} \right) $
and s represents the number of samples.
Overall, the development of VIEPR follows procedures shown in Figure 3. For brevity, the variational posterior of
$ \boldsymbol{\theta} $
is assumed to comply with normal distribution with only two parameters: mean and standard deviation denoted u and ρ. The training process aims to optimize u and ρ by gradient descent, instead of a scalar in EPR, meaning that the number of parameters merely doubles in comparison with EPR. Compared with other Bayesian methods such as Bayesian NNs (Zhang et al., Reference Zhang, Yin and Jin2022b), VIEPR has the same computational principle with gradient descent to estimate parametric posterior distributions of coefficients/weights. For practical applications, VIEPR facilitates the understanding of engineers by directly organizing related variables as polynomials in vast scenarios although with limited structures and fitting capability.
Figure 3. Schematic of VIEPR.
3. Data source
The database used in this study is extracted from (Zhang et al., Reference Zhang, Yin, Jin, Chan and Gao2021d) and includes 331 datasets, in which the
$ {C}_c $
of various clays was measured with three widely observed physical indices:
$ {e}_0 $
, LL, and Ip. Numerous empirical relations were built on these indices to identify
$ {C}_c $
, primarily because
$ {e}_0 $
directly affects the available space to be compressed, while LL and Ip reflect the range of soil moisture content associated with compression-induced plastic deformation. Nonetheless, these equations neglected underlying uncertainties and merely generated point estimation. Table 1 summarizes the statistics of indices in the database, which covers wide ranges of initial void ratio (
$ 0.66\le {e}_0\le 4.64 $
), liquid limit (
$ 25\%\le LL\le 166.2\% $
), and plasticity index (
$ 8\%\le Ip\le 113.9\% $
). They were then rescaled into the range (0, 1) using min–max normalization to eliminate scale effects before training.
Table 1. Statistics of soil indices
Note. STD = standard deviation.
Figure 4a presents the scatterplots of these soil indices, in which dense and sparse clusters can be easily distinguished. In this case, reliability evaluation becomes significant because accurate point estimations are elusive, due to insufficient knowledge of sparse clusters. On the other hand, the distributions of the indices are illustrated by the histograms in Figure 4b, in which they exhibit notably inhomogeneous distribution with positive values, similar to most observations (Camós et al., Reference Camós, Špačková, Straub and Molins2016; Jin and Yin, Reference Jin and Yin2020; Shi and Wang, Reference Shi and Wang2022; Tang and Phoon, Reference Tang and Phoon2019). Therefore, the natural logarithm is applied to the indices before training to ensure nonnegative predictions, similar to other published works (Ching et al., Reference Ching, Phoon and Wu2022; Löfman and Korkiala-Tanttu, Reference Löfman and Korkiala-Tanttu2021). This operation can be mathematically inverted in closed form. For instance, the median value of the predictions corresponds to the inverse logarithmic value of the mean of the normal distribution (Ching et al., Reference Ching, Phoon and Wu2022).
Figure 4. Statistics of soil indices. (a) Scatterplots. (b) Histograms.
Among the database, 80% was randomly drawn to train the model and the remaining 20% was employed to test the model. For comparison, all models adopt the same training and testing sets such that their performances can be fairly compared. To assess model performance, the commonly used mean absolute error (MAE) and the coefficient of determination (
$ {R}^2 $
) are computed as follows:
$$ \mathrm{MAE}=\frac{1}{n}{\sum}_{i=1}^n\left|{y}_i^m-{y}_i^p\right| $$
$$ {R}^2=1-\frac{\sum_{i=1}^n{\left({y}_i^m-{\overline{y}}_i^p\right)}^2}{\sum_{i=1}^n{\left({y}_i^m-{\overline{y}}_i^a\right)}^2} $$
where n is the total number of data points;
$ {y}_i^p $
,
$ {y}_i^m $
, and
$ {\overline{y}}_i^m $
are predicted, measured, and average measured quantity of interest, respectively. For probabilistic approaches, model reliability can be assessed by the ratio of the target quantity of interest falling inside the predicted 95% CI. Specifically, this ratio is expected to reach 95% for a desirable 95% CI in statistics, implying that a model would be unreliable if the ratio is less than 95% (He et al., Reference He, Zhang and Yin2024a, Reference He, Zhang, Yin and Goh2024b; Lyu et al., Reference Lyu, Hu and Wang2023; Zhao and Wang, Reference Zhao and Wang2018). The generalization refers to the ability to produce accurate and reliable predictions on unseen data and can be reflected by model performance on the testing set (Zhang et al., Reference Zhang, Yin and Jin2021c).
4. Development and evaluation of uncertainty models
4.1. Configurations for model development
The training of data-driven models was conducted using an algorithm under a given configuration. In this section, QRF and VIEPR were proposed as baseline algorithms for predicting
$ {C}_c $
, and their performances related to respective hyper-parameters. Nowadays, MAE integrated with 10-fold cross-validation (Stone, Reference Stone1974) has been commonly used as the loss function when optimizing hyperparameters in several data-driven modeling platforms (Zhang et al., Reference Zhang, Yin and Jin2021b, Reference Zhang, Yin and Jin2022b) and has also been used in the current study.
During the training process of QRF, each validation set is evaluated by trees that are trained on the remaining datasets to prevent overfitting (Breiman, Reference Breiman2001). Subsequently, the generated MAE is adopted as the loss function as shown in Equation (12). Similarly, the KL divergence with regularization term is employed as the loss function of VIEPR to optimize its sole hyperparameter nt. The hyperparameters of QRF and VIEPR are summarized in Table 2. For QRF, the numbers of trees, the candidate features, and the minimum number of samples for split (ntree, mtry, and min_sample_split) are optimized because they determine the size of forests and the growth of trees (Pei et al., Reference Pei, Qiu and Shen2023). As a result, the depth of each tree can be automatically determined once they are fixed. As a widely used evolutionary algorithm, particle swarm optimization (PSO) is employed to optimize these hyper-parameters and search for the elements of the exponent matrix in VIEPR. More details of PSO are introduced in Appendix A.
Table 2. Hyperparameters of QRF and VIEPR
The configuration of QRF and VIEPR is summarized in Table 3, where sufficient ranges of hyperparameters can be arbitrarily assigned and optimized to find the optimal one for predicting any quantity of interest. QRF determines the split by dichotomy to minimize MAE, while VIEPR is trained by gradient descent to minimize the KL divergence. The Adam optimizer is employed because of its effectiveness in handling nonstationary objectives (Kingma and Ba, Reference Kingma and Ba2014; Zhang et al., Reference Zhang, Yin and Jin2022b). QRF utilizes efficient dichotomy and only requires 500 epochs for hyperparameter optimization, whereas VIEPR updates the posterior coefficients over 20,000 epochs. Since the coefficients in VIEPR do not involve physical knowledge, their prior distributions are assigned a standard normal distribution. Finally, a total of 1000 MC samples are drawn to evaluate the KL divergence of VIEPR and are also employed for predictions.
Table 3. Configuration of QRF and VIEPR based models
4.2. Training processes of models
According to the configuration and data source described in the preceding sections, the training results of QRF and VIEPR are presented in Figure 5, illustrating the evolution of their losses. In Figure 5a, the loss value of QRF drops rapidly at the beginning and stabilizes after around 100 epochs, reaching its minimum. This corresponds to the optimized hyperparameters, that is, ntree = 180, mtry = 1, and min_sample_split = 1. Figure 6 presents one of its trees to illustrate the optimized QRF, where each node attaches the split value of the feature, the MAE value, and the number of samples inside. Starting with an initial MAE of 0.45 at the root node, the error gradually decreases as the tree branch undergoes random splits. This indicates that tree structures effectively minimize errors by distributing datasets to leaf nodes with similar features.
Figure 5. Evolution of loss value. (a) QRF. (b) VIEPR (number of transformed variables).
Figure 6. Structure of one tree of QRF.
Similarly, Figure 5b illustrates the loss values of VIEPR, using 1, 4, 7, and 10 transformed variables as examples. The VIEPR model with only one transformed variable quickly reaches a relatively higher convergence value at around 1000 epochs primarily because of limited fitting capability. By contrast, VIEPR models with multiple transformed variables attain relatively lower loss values through more epochs, due to larger search space. The loss values for models with 7 and 10 variables are close to each other, while the minimum loss value is attained by the model with four variables. This is mainly because apart from likelihood, the KL divergence also includes a regularization term to penalty model complexity, as shown in the first term of Equation (10). Consequently, the VIEPR model with four transformed variables reaches the minimum loss, and its expression is formulated by
$$ {y}_{\mathrm{VIEPR}}={w}_0+{w}_1{e}_0 L L+{w}_2\frac{LL}{I_p}+{w}_3 L L+{w}_4{e}_0 $$
where each coefficient in the VIEPR model follows a normal distribution and can be expressed by
$ {w}_0\sim $
N(−1.364, 0.062);
$ {w}_1\sim $
N(−0.001, 0.001);
$ {w}_2\sim $
N(−0.253, 0.026);
$ {w}_3\sim $
N(0.006, 0.001);
$ {w}_4\sim $
N(0.324, 0.027). Figure 7 illustrates the prior and posterior of coefficients in VIEPR, in which the posterior distribution of coefficients was derived in Equation (14). Utilizing the QRF and VIEPR models with the minimum losses, material indices can be predicted. Note that Equation (14) represents the expression in interior computations and requires the addition of a natural logarithmic base to generate nonnegative predictions, due to data preprocessing before training.
Figure 7. Prior and posteriors of coefficients in VIEPR.
4.3. Accuracy and reliability evaluation of models
Based on the above training results, developed QRF and VIEPR models were applied to predict
$ {C}_c $
values and were also compared with RF and EPR using the same training and testing sets (Zhang et al., Reference Zhang, Yin, Jin, Chan and Gao2021d). By comparison, QRF makes inferences based on the generated CDF, while RF relies on the average output of each tree. On the other hand, EPR followed the same procedure as VIEPR, but its coefficients are approximated by the least-squares method.
$ {C}_c $
is also computed by three commonly used empirical relations (Skempton and Jones, Reference Skempton and Jones1944; Sridharan and Nagaraj, Reference Sridharan and Nagaraj2000; Tiwari and Ajmera, Reference Tiwari and Ajmera2012), denoted three relations in Figure 8 as a benchmark.
Figure 8. Performance of models in terms of (a) MAE and (b) R2.
To assess model accuracy, Figure 8a and 8b presents the performance of the four models as well as three relations, in terms of MAE and R
2 values, respectively. Notably, these empirical relations showed significantly larger errors and less agreement with measurements, compared with others. For machine learning-based models, QRF and RF attained notably lower errors, reaching less than one-third of those generated by VIEPR and EPR on both training and testing sets. This difference can be attributed to model structures. Specifically, EPR-based models are fixed expressions using several variables, while RF-based models allow flexible random trees and splits to minimize errors. Compared to the results generated by RF (Zhang et al., Reference Zhang, Yin, Jin, Chan and Gao2021d), QRF observed significantly improved accuracy on both training and testing sets. This is primarily because QRF makes predictions based on the generated CDF incorporating distribution information in each leaf node, while RF merely relies on the average output of each tree. On the other hand, VIEPR showed a similar level of accuracy as EPR with a minor increase of error in the testing set, which primarily relates to the regularization term on model complexity in VIEPR. Overall, the trends observed in the R
2 values align with the discussion on the MAE values. Among the four models, QRF generated the highest R
2 and showed the greatest agreement with the measured
$ {C}_c $
values, demonstrating its superior capability in capturing correlations from limited data.
Regarding reliability evaluation, Figure 9a and 9b illustrates the results of QRF and VIEPR through a measured versus predicted
$ {C}_c $
(median) relationship, along with 95% CIs to examine model reliabilities. Specifically, the measured
$ {C}_c $
values outside 95% CIs on the training and testing sets are plotted as green triangles and circles, of which the number can reveal the reliability of models. Ideally, a 95% CI approximately covers 95% of measurements, and therefore, excessive measurements outside the interval indicate poor reliability and vice versa. Figure 9a shows that nearly all measured
$ {C}_c $
values fell in the 95% CI predicted by QRF, apart from one measurement in the testing set. Such a performance relates to the weighted sum of distribution information over an entire forest, therefore offering desirable reliability and generalization ability. In contrast, Figure 9b shows that VIEPR observes around 10% of measurements falling outside the 95% CIs, showing a slightly underestimated uncertainty due to limited model structures. Also, predictions of QRF are distributed denser around the diagonal than VIEPR, and this observation also aligns with what has been discussed earlier about R
2.
Figure 9. Comparison between measured and predicted
$ {C}_c $
with 95% CI generated by (a) QRF and (b) VIEPR.
To investigate the dependency of model accuracy and uncertainty on input variables, e
0 is taken as an example to illustrate its relationship between error and 95% CI in Figure 10. The results show that QRF obtained less error and uncertainty than explicit VIEPR along the distribution of e
0. Therefore, QRF is expected to perform more competitively than VIEPR to optimize cost-effectiveness in reliability-based design, due to reliable control of design parameters. For QRF, smaller uncertainty was found to appear in the densely distributed range of e
0, that is, 1.0–2.0, indicating that more datasets refine predicted CDF and enhance confidence level. In this regard, the understandable CDF using a simple weighted form allows engineers to easily understand the principle of interval estimates in tree structure-based models. By comparison, the uncertainty of VIEPR relies on the contribution of each variable in Equation (14), in which increased e
0 enlarges the predicted uncertainty of
$ {C}_c $
since they are positively related. Such a relationship is uniform with physics phenomena because a larger e
0 implies broader space to be compressed.
Figure 10. Relationship of predicted error and uncertainty on e0. (a) Training set. (b) Testing set.
With the regularization term, VIEPR generated a consistent error and mean standard deviation on the training set (0.054 and 0.075) than on the testing set (0.059 and 0.072). Since model accuracy and uncertainty relate to the mean and variance of polynomials, respectively, it was possible to observe such a result on randomly sampled training and testing sets. A key advantage of EPR models is their explicit expression, making them highly accessible and easy to apply in practice (Gomes et al., Reference Gomes, Gomes and Vargas2021; Nassr et al., Reference Nassr, Javadi and Faramarzi2018). Although not as accurate as QRF, the developed VIEPR can be readily formulated and even applied through Excel. For instance, Equation (14) can be used by engineers without expertise in data-driven modeling and also delivered into other scientific computing software. For small datasets, these baseline approaches can fuse data from multiple sources to generate predictions for material indices of interest, similar to other research works (Bozorgzadeh et al., Reference Bozorgzadeh, Harrison and Escobar2019; Ching and Phoon, Reference Ching and Phoon2020; He et al., Reference He, Zhang and Yin2024a; Wu et al., Reference Wu, Ching and Phoon2022).
5. Verification of generalization with application to creep index prediction
The previous results have demonstrated the capabilities of QRF and VIEPR with UQ, improving their applicability in engineering practice compared to deterministic methods. The findings indicate that QRF generates more reliable predictions, exhibiting superior accuracy and reliability using the same database. Moreover, its flexible structures and easily adjustable hyperparameters contribute to a simpler training process and application. Several platforms have already been developed to integrate RF with automatic hyperparameter optimization and model development (Zhang et al., Reference Zhang, Yin and Jin2021b, Reference Zhang, Yin and Jin2022b), ensuring its availability. Considering these highly proposing characteristics, QRF is suggested to further verify its generalization ability through additional material index prediction.
As a common design parameter,
$ {C}_{\unicode{x03B1}} $
is associated with multiple soil physical indices and involves substantial uncertainty.
$ {C}_{\unicode{x03B1}} $
holds significant importance in engineering practice as it affects the long-term performance of infrastructure such as settlement. Therefore, it is worthwhile to explore
$ {C}_{\unicode{x03B1}} $
for reliability-based design. The database was drawn from Zhang et al. (Reference Zhang, Yin, Jin and Chan2020) and includes 151 datasets, where the
$ {C}_{\unicode{x03B1}} $
of various clays was measured with three commonly observed variables:
$ {e}_0 $
, clay content CI and Ip. For these datasets, 80% was randomly drawn for model training, while the remaining 20% was reserved for testing. The same configuration of QRF used earlier is employed.
Figure 11a presents the predicted error of
$ {C}_{\unicode{x03B1}} $
along the distribution of
$ {e}_0 $
using QRF. The results show that QRF observed smaller MAE values, compared to 0.0012 and 0.0039 on the training and testing sets given by RF in the previous study (Zhang et al., Reference Zhang, Yin, Jin and Chan2020), demonstrating superior accuracy and generalization ability. Apart from point estimation, the 95% CI was provided in Figure 11a, which became broader in the sparsely distributed range of
$ {e}_0 $
, especially the range from 1.5 to 2.0 in the testing set. This means that QRF can perceive unknown information and generate a reasonable level of confidence when applied to reliability-based design. Figure 11b compares the predicted and measured
$ {C}_c $
values, in which the majority of measured
$ {C}_{\unicode{x03B1}} $
values fell in the predicted 95% CI with merely four outliers out of 31 datasets in the testing set. It is worth noting that the ratio of outliers did not strictly align with the expected 5%, which relates to insufficient datasets. These results showcase the effectiveness of QRF in enhancing accuracy and reliability on both training and testing sets, thus validating its generalization ability.
Figure 11. Performance of QRF in predicting
$ {C}_{\unicode{x03B1}} $
. (a) Accuracy. (b) Reliability.
6. Conclusions
This study has tailored two generic uncertain modeling approaches based on RF and EPR. Specifically, RF is integrated with quantile regression to extract distribution information in leaf nodes for UQ, while VI is employed to approximate parametric posterior distributions of coefficients in EPR. Thereafter, QRF and VIEPR were developed and applied to predict the soil compression index of clays.
The results indicated that QRF effectively reduced errors and covered almost all observations in the predicted 95% CI by integrating tree structures and distribution information inside. In contrast, VIEPR efficiently found parametric posterior distributions of coefficients to generate explicit solutions but observed greater errors and uncertainty due to restricted model structures. Therefore, QRF is suggested to optimize cost-effectiveness for reliability-based design due to efficient control of design accuracy and range. VIEPR is recommended to find user-friendly solutions even for those without expertise in data-driven modeling, offering enhanced reliability and wide applicability.
Compared with deterministic counterparts, QRF offered slightly improved accuracy than RF mainly because QRF generated a robust inference incorporating the distribution information in each leaf node, rather than merely relying on the average of each tree. Conversely, VIEPR showed a minor decrease in accuracy than EPR, primarily due to an additional penalty in model complexity. The dependency of model accuracy and uncertainty on input variables has been also investigated. Both QRF and VIEPR observed larger error and uncertainty in sparsely distributed ranges of input variables, meaning that they can recognize where their drawbacks exist and provide reasonable confidence in predictions.
Taking into account superior performance in predicting the compression index, QRF was further utilized to predict additional soil creep index using a new database. The accuracy and generalization of QRF were validated by comparing to RF in the previous study, while also incorporating interval estimates. It is worth noting that the proposed approach to quantify uncertainty is not limited to RF but also various tree structure-based models such as gradient boosting and extremely randomized trees. These highly promising characteristics strongly suggest that the proposed approaches have immense potential to enhance the data-driven modeling framework and can be further extended to a broader range of applications in scientific computing domains.
Data availability statement
Data are well reported in this study and will be made available from the corresponding author upon reasonable request.
Acknowledgements
The authors want to thank anonymous reviewers for their calibrations and suggestions.
Author contribution
Geng-Fu He: Writing—review and editing, Writing—original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation. Pin Zhang: Writing—review and editing, Visualization, Formal analysis, Methodology, Conceptualization. Zhen-Yu Yin: Writing—review and editing, Software, Supervision, Project administration, Funding acquisition.
Funding statement
This research was financially supported by the Research Grants Council (RGC) of Hong Kong Special Administrative Region Government (HKSARG) of China (Grant Nos.: 15220221, 15227923, 15229223, T22-607/24-N, E-PolyU501/24). The second author is supported by the Royal Society under the Newton International Fellowship.
Competing interest
The authors declare there are no competing interests.
Appendix A. Particle swarm optimization (PSO)
As an evolutionary algorithm, PSO searches for a globally optimal solution by iteratively updating a population of particles within a series of generations. Each particle represents a candidate solution and contains three components including position (X), velocity (V), and an objective function value. For the ith particle at the (k
$ + $
1)th generation, it is updated by (Zhang et al., Reference Zhang, Yin and Chen2022a).
$$ {V}_i^{k+1}=\mathrm{w}{V}_i^k+{c}_1{r}_1\left({P}_i^k-{X}_i^k\right)+{c}_2{r}_2\left({P}_g^k-{X}_i^k\right) $$
$$ {\displaystyle \begin{array}{c}{X}_i^{k+1}={X}_i^k+{V}_i^{k+1}\end{array}} $$
where
$ {c}_1 $
and
$ {c}_2 $
= acceleration coefficients;
$ w $
$ = $
an inertia weight;
$ {r}_1 $
and
$ {r}_2 $
$ = $
random numbers, distributed within the range [0, 1];
$ {P}_i $
denotes the optimal location of the ith particle, and
$ {P}_g $
denotes the globally optimal location among all particles. Table A1 summarizes the configuration of PSO used in this study, in which
$ w $
,
$ {c}_1 $
, and
$ {c}_2 $
are 1.0, 1.5, and 1.5, respectively, within commonly used ranges. PSO is assigned 20 particles and 100 generations to guarantee its exploration ability. Considering hyperparameters should be integers in this study, X is always approximated to the closest integer.
Table A1. Configuration of PSO
Note. Popsize = size of population; Maxgen = maximum number of generations.
Open access
Comments
No Comments have been published for this article.