Quantifying Statistical Uncertainty in Site
Investigation
Jianye CHING a , Shih-Hsuan WU b and Kok-Kwang PHOON c a
Professor, Dept of Civil Engineering, National Taiwan University, Taipei, Taiwan.
b
Postdoctoral fellow, Dept of Civil Engineering, National Taiwan University, Taipei, Taiwan.
c
Professor, Dept of Civil & Environmental Engineering, National University of Singapore, Singapore.
Abstract. Due to limited information in site investigation, it is not possible to obtain the actual values for the mean, standard deviation, and scale of fluctuation of a soil property of interest. The deviation between the estimated values and actual values is called the statistical uncertainty. There are two schools of thoughts on how to model the statistical uncertainty: frequentist thought and Bayesian thought. The purpose of this paper is to discuss their philosophical difference, to show how to quantify the statistical uncertainty based on these two distinct schools of thoughts, and to compare their performances. For the frequentist school of thought, the confidence interval will be used to quantify the statistical uncertainty, whereas the posterior probability distribution will be used for the Bayesian school of thought. Examples will be presented to compare the performances of these two schools of thoughts in terms of their consistencies. The results show that in general the Bayesian school of thought performs better in terms of consistency. In particular, the Markov chain Monte Carlo method is recommended when the information amount is very limited.
Keywords. statistical uncertainty, site investigation, frequentist, Bayesian, reliability
1. Introduction
One of the purposes of site investigation is to obtain information on the spatial distribution of geotechnical design parameters. Typically, the spatial distribution is expressed as a trend function and a zero-mean oscillating component about the trend. This oscillating component is called the spatial variability. The spatial variability has profound impact to the behavior of a geotechnical system, as illustrated in literature. Due to the importance of spatial variability, it is important to gather information on its statistical properties through site investigation. The spatial variability is typically characterized by its standard deviation (V) and scale of fluctuation (G). There are several techniques that can be employed to estimate V and G (in particular G), such as the method of moments (Uzielli et al. 2005; Dasaka and Zhang 2012; Firouzianbandpey et al. 2014; Lloret-Cabot et al. 2014), the fluctuation function method (Wickremesinghe and Campanella 1993; Cafaro and Cherubini 2002), the maximum likelihood (ML) method (DeGroot and Baecher 1993), and the Bayesian method (Wang et al.
2010). A common issue of these techniques is that the resulting V and G estimates are not identical to their actual values. The deviation between the estimated and actual values is called the statistical uncertainty.
Statistical uncertainty is present, because the amount of information collected in a site investigation is always limited (Phoon and Kulhawy 1999). The only way to eliminate the statistical uncertainty is to investigate the entire population. For site investigation, this means that the entire subsurface domain of interest must be explored, which is obviously not practical. Because of the statistical uncertainty, the V and G estimates based on different borehole/sounding data within the same soil layer will be different, although the actual V and G values are fixed numbers. The upshot is that the actual values of V and G cannot be known precisely in real world applications. In literature, case histories illustrating the statistical uncertainty can be found. One possible example is the clay site in Taranto, Italy (Cafaro and Cherubini 2002). The site consists of two fairly uniform clay layers (upper and lower layers). However, the G estimates for the lower layer obtained from five © 2015 The authors and IOS Press.
This article is published online with Open Access by IOS Press and distributed under the terms of the Creative Commons Attribution Non-Commercial License.
cone penetration test (CPT) soundings range from 0.185 m to 0.536 m. The actual G value cannot be known precisely.
Although the statistical uncertainty is well known, there are only a very limited number of previous studies focusing on this subject. Honjo and Setiawan (2007) studied the statistical uncertainty for the average of a property field along a certain depth. They were not concerned with the statistical uncertainties in V and G: V and G are prescribed in their analysis. Instead, they are concerned with the unknown constant trend (or the mean value P). The rationale for prescribing V and G is based on the argument that the number of data points in site investigation is typically insufficient to estimate these second-order statistics. To circumvent this practical difficulty, Honjo and Setiawan (2007) suggested that conservative values for V and G can be assumed based on previous studies and experiences. In other previous studies, the statistical uncertainty may have been treated implicitly by lumping with other sources of uncertainties. For instance, Wang et al. (2010) proposed a Bayesian method of characterizing the uncertainties in (P,V,G). The resulting posterior probability distribution actually incorporates the statistical uncertainties in (P,V,G).
Although it is possible to prescribe (V,G) based on previous studies and experiences, as suggested by Honjo and Setiawan (2007), there are practical difficulties for doing so. First of all, (V,G) values in literature vary in a wide range. For instance, the coefficient of variation (COV = V/P) for the undrained shear strength varies from 0.1 to 0.6 (Phoon and Kulhawy 1999). Its vertical G varies from 0.8 to 6.2 m, whereas its horizontal G is known in a very limited way (only 3 studies collected by Phoon and Kulhawy 1999). It is not simple to select suitable conservative values for (V,G) based on the above wide ranges. Moreover, G may depend on the problem scale (Fenton 1999) and (V,G) may depend on the adopted trend function and the sampling interval as well (Cafaro and Cherubini 2002). The scale considered in previous studies may not be similar to the scale applicable for the geotechnical project at hand. The trend function and sampling interval studied in the literature may not be
applicable to the conditions in the project at hand. There are strong reasons, either practical or analytical, to estimate (V,G) based on the local site investigation data.
The purpose of this paper is to quantify the statistical uncertainties for the 3 parameters (P,V,G) which are required for the second-moment characterization of a random field. To achieve this, we will first demonstrate the statistical uncertainties in (P,V,G) using simulated spatial variability examples. Second, two schools of thoughts, frequentist and Bayesian, will be adopted to interpret and quantify the statistical uncertainties. Analytical tools for quantifying the statistical uncertainties will be also introduced. For the frequentist school of thought, the tool is the confidence interval. For the Bayesian school of thought, the tool is the posterior probability distribution. Finally, the two schools of thoughts will be compared on the basis of consistency: which school of thoughts can provide confidence interval or posterior probability distribution that is more consistent to the actual values of (P,V,G).
2. Random Field Simulation
Vanmarcke (1977) proposed that the spatial variability can be modeled as a random field. Among random field models, stationary random fields are widely used due to their simplicity and possibly the only practical version that can be characterized from limited data. Three parameters are required to characterize a second-order stationary random field model: (a) mean (P), (b) standard deviation (V), and (c) auto-correlation function (ACF). For a one-dimensional random field W(z), where W is the property field and z is the location (e.g., depth), the auto-correlation function is defined to the correlation between two locations with 'z apart:
CV W z ,W z+z z = Var W z Var W z+z ª º ¬ ¼ ª º ª º ¬ ¼ ¬ ¼ (1)
where Var(.) denotes variance; CV(.,.) denotes covariance. The most popular ACF is the single exponential model (Vanmarcke 1977):
z exp ' G2 z (2)
where G is the scale of fluctuation (SOF).
A one-dimensional (1D) stationary normal random field with mean = P, standard deviation = V, and scale of fluctuation = G can be simulated by the following equation:
W P 1 V L U (3)
where W = [W(z1) W(z2) …W(zn)]T contains the point processes of W(z) at the n locations of interest; U = [U1 U2 …Un]T contains independent standard normal samples; 1 = [1 1 … 1]T; L is the Cholesky decomposition of the correlation matrix R. In this study, W = [W(z1) W(z2) …W(zn)]T represents the vertical profile of a CPT sounding in a soil layer. dz = zi+1 – zi is the vertical sampling interval, and D = zn – z1 is the total depth of the CPT. It is assumed that the CPT data has been transformed (e.g., detrended, or converted to the normalized cone tip resistance, and possibly taken logarithm, etc.) so that W is a realization from the underlying stationary normal random field. Without loss of generality, the actual values of (P,V,G) are taken to be (0,1,1). Two normalized parameters are used to characterize the total depth D and sampling interval dz: nD = D/G and ' = dz/G.
3. Frequentist and Bayesian Schools of Thoughts
3.1. Frequentist school of thoughts
From the frequentist point of view, the W data is a random realization from the population. Central to the frequentist thought are the estimators of (P,V,G), which are functions of W. (V,G) will be taken logarithm throughout this study. Because W is random, the estimators of (P,ln(V),ln(G)) are random variables. There are various estimators. Among them, the maximum likelihood (ML) estimator has several desirable asymptotic properties, including asymptotically unbiased, asymptotically normally distributed, asymptotically minimum variance, invariant, sufficient, and consistent (DeGroot and Baecher
1993). The ML estimators of (P,ln(V),ln(G)) maximize the following likelihood function:
T 1 2 1 2 n 2 2 f , ln , ln 1 1 e 2 P P V P V G S V W 1 R W 1 W R (4)
The ML estimators of (P,ln(V),ln(G)) are denoted by (PML,ln(VML),ln(GML)).
Frequentist asymptotic approximation method
(PML,ln(VML),ln(GML)) are random variables because they depend on W and W is random. The multivariate PDF of (PML,ln(VML),ln(GML)) is called the sampling distribution. According to Mardia and Marshall (1984), (PML,ln(VML), ln(GML)) are asymptotically (multivariate) normally distributed with mean equal to the actual values of (P,ln(V),ln(G)) and the covariance matrix C equal to the inverse of the information matrix, which is the expectation for the negative of the Hessian matrix for the log-likelihood. In principle, the sampling distribution has mean equal to the actual values of (P,ln(V),ln(G)) and C must be evaluated at the actual values of (P,ln(V),ln(G)). However, the actual values of (P,ln(V),ln(G)) are unknown. In practice, the actual values of (P,ln(V),ln(G)) are replaced by (PML,ln(VML),ln(GML)).
As mentioned earlier, the confidence interval/region plays a central role in the frequentist school of thought. The 95% confidence region for the actual values of (P, ln(V),ln(G)) can be constructed based on the sampling distribution of (PML,ln(VML),ln(GML)). This 95% confidence region is a three-dimensional (3D) ellipse defined as
T ML ML 1 ML ML ML ML ln ln ln ln 7.815 ln ln ln ln P P P P ª º ª º « V V » « V V »d « » « » « G G » « G G » ¬ ¼ ¬ ¼ C (5)
where 7.815 is the 0.95 quantile of the Chi-squared variable with 3 degrees of freedom.
Limited number of data points (bootstrapping) The above multivariate normal sampling distribution and the resulting confidence region are valid for an asymptotic case, which means nD is large. However, for site investigation, it is common to have small nD. With this limited amount of data, the above asymptotic approximation is not proper. The bootstrapping method (Efron and Tibshirani 1993) is a general framework that can be applied to non-asymptotic cases to draw random samples from the sampling distribution of (PML,ln(VML),ln(GML)). The steps for bootstrapping are as follows:
1. Re-sampling the data indices (1, 2, …, n). Note that after the re-sampling, there may be repetitive indices. Suppose the re-sampled indices are (5, 1, 9, …, 1). The re-sampled z data points are therefore z5, z1, z9, …, z1 and the re-sampled W data points are W(z5), W(z1), W(z9), …, W(z1).
2. Evaluate (PML,ln(VML),ln(GML)) using the re-sampled pairwise z-W data points.
3. Repeat steps 1 and 2 to obtain nB samples of (PML,ln(VML),ln(GML)). Note that nB is distinctive from n.
These nB samples are approximately from the sampling distribution of (PML,ln(VML), ln(GML)).
3.2. Bayesian school of thoughts
The main philosophical difference between the frequentist and Bayesian schools of thoughts is that nothing is random in the Bayesian thought – a quantity that cannot be determined precisely is treated as uncertain (lack of precise information) rather than random. From the Bayesian point of view, the parameters (P,ln(V),ln(G)) are fixed numbers but they are uncertain as well. A prior PDF, denoted by f(P,ln(V),ln(G)), is used to describe the a priori knowledge for (P,ln(V), ln(G)). It specifies the plausibility of various combinations of (P,ln(V),ln(G)) before the W data is obtained. By specifying the prior PDF f(P,ln(V),ln(G)), (P,ln(V),ln(G)) are not random in Bayesian point of view. They are uncertain. The likelihood function in Eq. (4) describes the plausibility of W, given the full knowledge of (P,ln(V),ln(G)). The posterior (updated) PDF can then be determined using the Bayes rule:
f , ln , ln f , ln , ln f , ln , ln f P V G P V G P V G W W W (6)
where f(W), the multivariate PDF of W, is a normalizing constant that does not depend on (P,ln(V),ln(G)). The posterior PDF f(P,ln(V),ln(G) |W) describes the plausibility of various combinations of (P,ln(V),ln(G)) after the W data is incorporated. The (P,ln(V),ln(G)) value that maximizes the posterior PDF f(P,ln(V),ln(G)|W) is called the maximum a posteriori (MAP) estimate, denoted by (PMAP,ln(VMAP),ln(GMAP)).
In the frequentist point of view, it does not make sense to speak about the PDF of (P,ln(V), ln(G)), because they are not random. Nonetheless, in the Bayesian thoughts, it makes sense to speak about the PDF of (P,ln(V),ln(G)), because they are uncertain. In particular, the posterior PDF describe the plausibility of various combinations of (P,ln(V),ln(G)) after W is obtained. The posterior PDF is central to the Bayesian thought. Asymptotic approximation method (Laplace) There is an asymptotic normal approximation for f(P,ln(V),ln(G)|W): the Laplace approximation (Bleistein and Handelsman 1986). The normal PDF that approximates f(P,ln(V),ln(G)|W) has mean = (PMAP,ln(VMAP),ln(GMAP)) and the covariance matrix = CMAP = the inverse for the negative of the Hessian for the log-likelihood, evaluated at (PMAP,ln(VMAP),ln(GMAP)). In the Bayesian framework, there is a term called the “Bayesian confidence region” analogous to the frequentist 95% confidence region (Jaynes 1976). For the asymptotic normal posterior PDF, the Bayesian 95% confidence region is an 3D ellipse described by an inequality similar to Eq. (5), with center replaced by (PMAP,ln(VMAP),ln(GMAP)), and covariance matrix
replaced by CMAP.
Limited number of data points
The Laplace approximation for the posterior PDF f(P,ln(V),ln(G)|W) is valid only for asymptotic cases with large nD. The Markov chain Monte Carlo (MCMC) method (Gilks et al. 1995) is a general framework of obtaining samples from
f(P,ln(V),ln(G)|W) regardless of the amount of information. The steps for a simplified MCMC algorithm (Metropolis et al. 1953) are as follows: 1. The initial time step t = 0, and the initial
samples are (P0,ln(V0),ln(G0
)). These sample values can be arbitrary. In this study, the MAP estimates (PMAP,ln(VMAP),ln(GMAP)) are taken to be the initial samples.
2. Draw the candidate sample (Pc,ln(Vc),ln(Gc)) from a joint normal PDF centering at (P0
, ln(V0),ln(G0)) and with a covariance matrix 66. This multivariate normal PDF is called the proposal PDF. Theoretically, the covariance matrix 6 can be arbitrary. But in reality, the choice of 6 will affect the efficiency of MCMC. It is found that the covariance matrix CMAP is a satisfactory choice for 6. This is adopted in the current study.
3. Compute the r ratio:
c c c c c c 0 0 0 0 0 0 f , ln , ln f , ln , ln r f , ln , ln f , ln , ln P V G P V G P V G P V G W W ( ) 4. Accept (P1,ln(V1),ln(G1)) = (Pc,ln(Vc),ln(Gc )) with probability min(1, r). If (Pc,ln(Vc),ln(Gc
)) is not accepted, repeat the previous sample: (P1,ln(V1),ln(G1)) = (P0,ln(V0),ln(G0
)).
5. Incrementally cycle steps 2~4 for t = 1, 2, …, nT times to obtain (P0,ln(V0),ln(G0)) … (PnT,ln(VnT),ln(GnT
)) samples.
When t is sufficiently far away from the initial time t = 0, the above process will reach a stationary state. When this happens, the samples (Pt,ln(Vt),ln(Gt)) are distributed as f(P,ln(V),ln(G)|
W). However, when t is small, the stationary
state has not been reached. (Pt,ln(Vt),ln(Gt)) can be affected by the initial position (P0,ln(V0),ln(G0
)). This initial non-stationary period is called the burn-in period. The samples within the burn-in period must be discarded.
4. Comparison of Consistencies
In this section, we compare the consistencies of the four methods (two frequentist and two Bayesian methods) introduced above. For all methods, the consistency is quantified based on how often the actual values of (P,ln(V),ln(G)),
namely (0,0,0), are within the 95% confidence region. The confidence regions for the four methods are summarized below:
1. (Frequentist asymptotic) The 95% confidence region based on the asymptotic approximation is adopted. This confidence region is the 3D ellipse described in Eq. (5). 2. (Frequentist bootstrapping) The 95%
confidence region based on bootstrapping samples is adopted. This confidence region is a 3D ellipse described by an inequality similar to Eq. (5), but the center of the ellipse is replaced by the sample mean of the bootstrapping samples, and the covariance matrix replaced by the sample covariance of the bootstrapping samples. It must be understood that the “cloud” of the bootstrapping samples may not resemble an ellipse. The elliptical form is used for simplicity and also for consistency to the asymptotic confidence region.
3. (Bayesian asymptotic) The 95% Bayesian confidence region based on the asymptotic posterior PDF (Laplace approximation) is adopted. This confidence region is the 3D ellipse described by an inequality similar to Eq. (5), but the center of the ellipse is replaced by the MAP estimates, and the covariance matrix replaced by CMAP.
4. (Bayesian MCMC) The 95% Bayesian confidence region based on MCMC samples is adopted. This confidence region is a 3D ellipse described by an inequality similar to Eq. (5), but the center of the ellipse is replaced by the sample mean of the MCMC samples, and the covariance matrix replaced by the sample covariance of the MCMC samples. It must be understood that the cloud of the MCMC samples may not resemble an ellipse. The elliptical form is used for simplicity and also for consistency to the asymptotic confidence region.
A method is considered to perform ideally if the chance that the ellipsoid does not contain (0,0,0) is exactly 0.05 (1 – 0.95 = 0.05). Random simulations are performed to understand how this chance changes with respect to nD. One hundred realizations of random fields with (P,V,G) = (0,1,1) and with ' = 0.25 (dz = 0.25) are simulated. Each random field realization is sampled with total depth = D, and the sampled 7
random field is used to construct the above four ellipsoids. Whether the four ellipsoids contain (0,0,0) can be easily by the above inequalities. The chance that the ellipsoid does not contain (0,0,0) is simply [number of realizations where the ellipsoid does not contain (0,0,0)]/100.
Figure 1 shows how this chance varies with respect to nD for the above four methods. It is evident the Bayesian MCMC method outperforms the other three methods because the chance for the MCMC method is the closest to 0.05. The two asymptotic methods (frequentist and Bayesian) have similar performances. The bootstrapping method has the poorest performance. When nD is large (nD t 20), the two asymptotic methods have performances comparable to the MCMC method. However, when nD is small, the two asymptotic methods are obviously outperformed by the MCMC method. This is significant: as mentioned earlier, it is quite common in geotechnical engineering practice to encounter cases with small nD. This makes the MCMC method the most consistent method for geotechnical engineering practice.
Figure 1. Chance that the ellipsoid does not contain (0,0,0).
5. Conclusions
The MCMC method (Bayesian) is the most consistent method for the estimation of (P,V,G), in the sense that the actual values of (P,V,G) are encompassed by the “cloud” of the MCMC samples with a high chance. This conclusion holds regardless of nD. In particular, for cases with small nD (thin layers), the MCMC method is the only consistent method. For cases with
large nD (nD t 20) (thick soil layers), both asymptotic approximation methods (frequentist and Bayesian) perform equally well to the MCMC method. The bootstrapping method (frequentist) is not consistent, even for cases with large nD.
References
Bleistein, N., Handelsman, R. (1986). Asymptotic Expansions of Integrals. Dover, New York.
Cafaro, F., Cherubini, C. (2002). Large sample spacing in evaluation of vertical strength variability of clayey soil.
ASCE Journal of Geotechnical and Geoenvironmental Engineering, 128(7), 558-568.
Dasaka, S. M., Zhang. L. M. (2012). Spatial variability of in situ weathered soil. Geotechnique, 62(5), 375-384. DeGroot, D. J., Baecher, G. B. (1993). Estimating
autocovariances of in-situ soil properties. ASCE Journal
of Geotechnical Engineering, 119(1), 147-166.
Efron, B., Tibshirani, R. (1993). An Introduction to the
Bootstrap. Boca Raton, FL. Chapman and Hall/CRC.
Firouzianbandpey, S., Griffiths, D.V., Ibsen, L.B. Anderson, L.V. (2014). Spatial correlation length of normalized cone data in sand: Case study in the north of Denmark.
Canadian Geotechnical Journal, 51(8), 844-857.
Gilks, W.R., Richardson, S. Spiegelhalter, D. (1995). Markov Chain Monte Carlo in Practice. Boca Raton, FL. Chapman and Hall/CRC.
Honjo, Y., Setiawan, B. (2007). General and local estimation of local average and their application in geotechnical parameter estimations. Georisk, 1(3), 167-176. Jaynes, E. T. (1976). Confidence intervals versus Bayesian
intervals. In Foundations of Probability Theory, Statistical Inference, and Statistical Theories of Science,
(Eds.: Harper and Hooker), Dordrecht, D. Reidel. Lloret-Cabot, M., Fenton, G.A. Hicks, M.A. (2014). On the
estimation of scale of fluctuation in geostatistics.
Georisk, 8(2), 129-140.
Mardia, K.V., Marshall, R.J. (1984). Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 71(1), 135-146.
Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21, 1087-1092.
Uzielli, M., Vannucchi, G. Phoon, K.K. (2005). Random field characterisation of stress-normalised cone penetration testing parameters. Geotechnique, 55(1), 3-20. Vanmarcke, E. H. (1977). Probabilistic modeling of soil
profiles. ASCE Journal of Geotechnical Engineering, GT11, 1227-1246.
Wang, Y., Au, S.K. Cao, Z. (2010). Bayesian approach for probabilistic characterization of sand friction angles.
Engineering Geology, 114, 354-363.
Wickremesinghe, D., Campanella, R. G. (1993). Scale of
fluctuation as a descriptor of soil variability. Probabilistic Methods in Geotechnical Engineering
(Eds.: Li and Lo), Balkema, Rotterdam, The Neterlands, 233-239.
