Extreme Value Analysis of Statistically Independent Stochastic Variables

Yongho Choi; Seong Mo Yeon; Hyunjoe Kim; Dongyeon Lee

doi:10.26748/KSOE.2018.093

Abstract

An extreme value analysis (EVA) is essential to obtain a design value for highly nonlinear variables such as long-term environmental data for wind and waves, and slamming or sloshing impact pressures. According to the extreme value theory (EVT), the extreme value distribution is derived by multiplying the initial cumulative distribution functions for independent and identically distributed (IID) random variables. However, in the position mooring of DNVGL, the sampled global maxima of the mooring line tension are assumed to be IID stochastic variables without checking their independence. The ITTC Recommended Procedures and Guidelines for Sloshing Model Tests never deal with the independence of the sampling data. Hence, a design value estimated without the IID check would be under- or over-estimated because of considering observations far away from a Weibull or generalized Pareto distribution (GPD) as outliers. In this study, the IID sampling data are first checked in an EVA. With no IID random variables, an automatic resampling scheme is recommended using the block maxima approach for a generalized extreme value (GEV) distribution and peaks-over-threshold (POT) approach for a GPD. A partial autocorrelation function (PACF) is used to check the IID variables. In this study, only one 5 h sample of sloshing test results was used for a feasibility study of the resampling IID variables approach. Based on this study, the resampling IID variables may reduce the number of outliers, and the statistically more appropriate design value could be achieved with independent samples.

Abbreviations

ACF: autocorrelation function

EVA: extreme value analysis

GPD: generalized Pareto distr.

GEV: generalized extreme value

CI: confidence interval

POT: peaks-over-threshold

MPM: most probable maximum

MOM: method of moments

MLE: maximum likelihood estimation method

PPC: Pearson correlation coefficient = R

IID: independent, identically distributed

PACF: partial ACF

EVT: extreme value theory

MEF: mean excess function

IQR: interquartile range

SE: standard error

P-P: probability-probability

Q-Q: quartile-quartile

LS: least squares methods

1. Introduction

The extreme value theory (EVT) based on the Fisher–Tippett theorem (Fisher and Tippett, 1928) describes the distribution of extremes and rare events, especially in highly nonlinear slamming or sloshing impact data and long-term environmental applications. An inference about extremes is usually difficult to achieve because of the lack of data. In this case, the tail behavior of the distribution is of importance. There are mainly two popular extreme value models in the EVT: the block maximum model and peaks-over-threshold (POT) method. According to the 1st EVT theorem, which is referred to as the Fisher–Tippett–Gnedenko or Fisher–Tippett theorem, a distribution of the block maxima yields a generalized extreme value distribution. The 2nd EVT theorem, which is referred to as the Pickands–Balkema–de Haan theorem (Balkema and Haan, 1974; Pickands, 1975), is used in the tail-fitting for the POT. The POT data also consist of two distributions: a Poisson distribution for the number of events per given time period, and a generalized Pareto distribution (GPD) for the size of the exceedances.

This paper will focus only on the EVT for the design value for the sloshing impact pressures obtained from a sloshing model test. It can be understandably applied to obtain a design environmental condition. In order to be able to apply the EVT, the observations need to be independent and identically distributed (IID). Hence, IID is to be checked for measured sloshing impact pressures. But DNVGL (2015) even describes that sampled global maxima of mooring line tension are assumed to be independent stochastic variables without checking IID. ITTC (2017) never covers the independence of sampling pressures in a sloshing model test. Consequently, the assumption that consecutive data are IID may not be theoretically well justified, and the common positive correlation of successive data may lead to conservative estimates by simply regarding observations far away from a Weibull or generalized Pareto distribution as outliers.

This study first checked the IID of the sampled data. If random variables do not pass the IID check, resampling schemes that use the block maxima approach for the GEV (Generalized extreme value) and POT approach for the GPD are used. A partial autocorrelation function (PACF) is used to check automatically the IID variables. Finally, the resampled IID data are used to estimate the extreme design values. For the feasibility study of the resampling IID variables approach, only one 5 h sample of sloshing test results was used in this study.

2. Theoretical Background

2.1 Extreme value theory

Fisher and Tippett (1928) developed and organized a GEV using three asymptotic limit distributions: Gumbel (ξ=0), Fréchet (ξ>0), and Weibull (ξ<0) distributions for IID variables (Ochi, 1990). Gnedenko (1943) and Gumbel (1954) gave the credit for the GEV or convergence to the type theorem (Choi 2016).

(1)

GGEV(x;μ,σ,ξ)={exp[−(1+ξ(x−μ)/σ)−1/ξ]ξ≠0exp[−exp(−(x−μ)/σ)]ξ=0

in which μ, σ, and ξ represent the location (“Shift” in figures), scale, and shape parameters, respectively. The r-return period value x_r is obtained by,

(2)

xr={μ^+σ^ξ^−1[−(ln(1−r−1))−ξ^−1]ξ≠0μ^−σ^ln(−ln(1−r−1))ξ=0

The symbol ^ represents the estimation. For the sloshing test, the r– return period in model scale for the Hr-hour in real scale is obtained by r=Hr×3600/dt/λ=Hr/D×Np and dt=D×3600/Np/λ, where D, N_p, dt, and λ are the test duration in hours, number of peak pressures, average time interval of impacts in seconds, and real-to-model scale ratio, respectively.

The maximum likelihood estimation (MLE, see sec. 2.4) method is used to estimate the GEV parameters by maximizing the likelihood function (L) or by minimizing the negative log-likelihood function (-l).

(3)

L(μ,σ,ξ;x)=∏ingGEV(x)=∏in∂GGEV(x)∂x=1σn∏in[1+ξxi−μσ]−(1+1ξ)∏inexp[−(1+ξxi−μσ)−1ξ]

(4)

−l(μ,σ,ξ;x)=−lnL(μ,σ,ξ;x)=nlnσ+(1+1ξ)∑inln[1+ξxi−μσ]+∑in(1+ξxi−μσ)−1ξ

A Hessian matrix H is defined to obtain the SE (standard error = stdev of sample mean = stdev/N) for the CI (confidence interval), in which stdev is a standard deviation. It is a square matrix of 2nd order partial derivatives of a scalar-valued negative log-likelihood function. The estimation symbol ^ is dropped hereafter. An inverse matrix of H yields the variance-covariance matrix V. The SEs of the distribution parameters are obtained by the square root of the diagonal terms of V.

(5)

V=H(−l)−1,H(−l)=(−∂2l∂μ2−∂2l∂σ∂μ−∂2l∂ξ∂μ−∂2l∂σ∂μ−∂2l∂σ2−∂2l∂ξ∂σ−∂2l∂ξ∂μ−∂2l∂ξ∂σ−∂2l∂ξ2)

Because the SEs of the r-return period extreme values are more practical than the SEs of the distribution, the Delta (∇) method (Xu and Long, 2005) for the GEV is used. Hence, the SE is the square root of the variance of the r-return period extreme value x^r.

(6)

SEx^r=(Var(x^r)≈∇xrT⋅V⋅∇xr, yr=−ln(1−r−1)

(7)

∇xr=(∂xr∂μ∂xr∂σ∂xr∂ξ)=(1−ξ−1(1−yr−ξ)σξ−2(1−yr−ξ)−σξ−1yr−ξlnyr)

Next, Balkema and Haan (1974) and Pickands (1975) verified that the exceedances over a sufficiently large threshold value u approximately yield the GPD when the IID variables follow a GEV distribution. The GPD is derived using the conditional exceedance distribution function for IID random variables.

(8)

GGPD(y)={1−(1+ξy/σu)−1/ξξ≠01−exp(−y/σu)ξ=0, σu=σ+ξ(u−μ)

The r-return period value x_r for the GPD is shown in the equation below. In this equation, n_y is the number of measured points per unit time (1 h for the sloshing test, 1 year for the environment) and λ^u=Pr(X>u)=1-F(u) is the ‘Rate’ (empirical threshold exceedance rate). In the case of sloshing tests, r=r'ny=Hr×3600/dt/λ. In the case of the environment, r=r'ny=r'×365.25×24/∆, where r^′ and Δ are the number of years and measuring interval in hours, respectively.

(9)

xr={u+σ^ξ^−1[(r′nyλ^u)−ξ^−1]ξ≠0u−σ^ln(r′nyλ^u)ξ=0

The negative log-likelihood function (-l) for GPD is,

(10)

−l(μ,σ,ξ;y)=nlnσ+(1+ξ−1)⋅∑in(−ξyi/σ)

The Delta method for the GPD of the SEs for the r-return period extreme values is shown in the equation below. Hereafter, t←r'ny, λ←λ^u.

(11)

∇xr=(∂xr∂μ∂xr∂σ∂xr∂ξ)=(σt(tλ)ξ−1ξ−1[(tλ)ξ−1]−σξ−2[(tλ)ξ−1]+σξ−1(tλ)ξln(tλ))

The determination of threshold u is essential for the GPD because the POT approach provides a model for IID exceedances over a high threshold. There are various capable methods for determining u, such as the MEF (mean excess function or MRL (mean residual life), most popular method (Kim and Song, 2012)) plot, parameter stability plot, dispersion index plot, rule of thumb, and multiple-threshold model method (Bommier, 2014). Because the MEF is a linear function of u, the starting point of linearity found in the MEF plot is to be a proper value as the threshold value.

(12)

MEF=e(u)=E[X−u|X>u]=(σ+ξu)/(1−ξ)

2.2 Other distributions

In the conventional analysis procedures of the sloshing test to date, a 3-parameter Weibull distribution, which is one of the GEV distributions when ξ<0, or GPD is used to estimate the extreme design value. The 3-parameter Weibull for the highest 10% of the data as a popular industry practice was also studied. It is a hybrid approach with the POT and GEV instead of the GPD.

2.3 Autocorrelation and partial autocorrelation

The data for EVT should be IID stochastic variables. Because zero correlation shows that the variables are linearly independent, the independency is checked by the correlation. A coefficient of the autocorrelation function (ACF) r_x(τ) shows the correlation with lag τ between x(t) and x(t-τ) in a given time series. The PACF p_x(τ, τ) is defined by the autocorrelation between x(t) and x(t-τ) in a given time series, with the linear dependence of x(t) on x(t-1), x(t-2), ..., x(t-τ-1) removed.

(13)

rx(τ)=∑(t=1)(n−τ)(xt−x¯)(xt+τ−x¯)/(∑t=1n(xt−x¯)2

(14)

px(τ,τ)=cov(xi−xi−τ|xi−1,⋯,xi−τ+1)var(xi|xi−1,⋯,xi−τ+1)⋅var(xi−τ|xi−1,⋯,xi−τ+1)

Where

(15)

{px(1,1)=rx(1)for τ=1px(τ,τ)=rx(τ)−∑(j=1)(τ−1)px(τ−1,j)rx(τ−j)1−∑(j=1)(τ−1)px(τ−1,j)rx(j)for τ=2,3,⋯,npx(τ,j)=px(τ−1,j)−px(τ,τ)⋅px(τ−1,τ−j)

A program for all the statistical analyses except for CI for sloshing impact pressures was realized in the Python programming language in this study, rather than the ‘R’ programming language and in-house code used in Choi (2016) for the extreme environmental condition.

2.4 Parameter estimation methods of distribution

When fitting probability distributions to data, DNVGL (2017) recommends several fitting techniques such as the method of moments (MOM), least squares methods (LS), and maximum likelihood estimation (MLE). MLE is used to fit the data on GEV and GPD. Because MLE for a Weibull is applicable only when the shape parameter is greater than one (Smith, 1985), and the shape parameter for the sloshing peak pressures is usually smaller than one, the MOM was applied for a Weibull in this study. Because the purpose of this study was to check the IID for statistical distribution models, a comparison of the fitting techniques was not carried out in this study.

3. Sloshing Impact Pressures and Conventional Fitting

Special care is needed to extract the sloshing impact pressures directly from the sloshing model test because of various causes of interference such as electronic noise and meaningless small pressures in the measured pressure time record. An initial POT method with an initial time window that is usually used in GTT is applied to cope with it (Fillon et al., 2011; Kim, 2017). Because the sloshing pressures sampled with this screening approach does not include the checking IID, an additional POT for GPD and block maxima for GEV for IID random variables were introduced in this study.

A sloshing model test with a model scale ratio λ = 50 for the No. 2 tank of a 160 k (160,000 m³) conventional LNG carrier was conducted by filling it to 30% with water and air under 1 year North Atlantic beam sea condition. This study dealt with the sample pressure record from a single Kistler piezo-type pressure sensor at only one location. Fig. 1 shows the extracted 764 N_p sloshing peak pressures with an initial threshold of 2.5 kPa and time window of 0.2 s, which is hereafter defined as one block for the block maxima approach in the model scale. The measured maximum pressure was 244.7 kPa in model scale. The total test duration was 2545.6 s in model scale (D = 5h in real scale); dt was 3.33 s. Fig. 2 shows the results of the conventional procedure for the 3-parameter Weibull fitting with all 764 data points. The estimated 3 h MPM (Q = 1/r = 0.002) as a return level for the 3 h return period is 143.6 kPa, and it may be underestimated considering the last two large pressures, where r = 3h/5h×764 = 458.4. A value of 143.6 kPa (75.4 bar) for water yields 69.4 kPa (35.4 bar) for LNG in real scale. The maximum value of the abscissa is obtained by multiplying the maximum measured data by 1.5.

The Tukey (1977) method using quartiles is used to check the outliers. An outlier from the fitting curve is determined based on a distribution of residuals obtained by subtracting the theoretical values of a model’s extreme distribution from samples with the assumption of normally distributed residuals. Then, the outsiders of Q₁-γ×IQR to Q₃+γ×IQR are defined as the outliers, where Q₁, Q₃, and IQR represent the first and third quartiles and a interquartile range Q₃-Q₁ of residuals, respectively. γ is a constant value traditionally given as 1.5 by Tukey. However, in this paper, an appropriate value will be found, because smaller residuals yield a larger γ. In the case of Fig. 2, the last two large pressures are considered as outliers with IQR = 4.769 and γ = 5.

Because the conventional procedure does not check the IID, no one makes sure whether or not the 3 h MPM of 143.6 kPa is correct. To avoid this, the linear independence is checked by PACF, as shown in Fig. 3. Two graphical diagnostic plots such as the P-P plot and Q-Q plot, as shown in Fig. 4, are also used for the goodness of fit. In the Q-Q plot, R is the Pearson correlation coefficient (PCC), with a range of -1.0 to 1.0.

Because one point at time lag τ = 18 of PACF in Fig. 3 exceeds the 95% CI, and the plots in Fig. 4 look bad (R = 0.982 in the Q-Q plot), one can say that all 764 data points are not IID. Hence, the 3 h MPM of 143.6 kPa is unconvincing.

Next, the other conventional procedure is GPD fitting for the largest 8% of the data as a rule of thumb, as shown in Fig. 5 and Fig. 6, with a GPD threshold u = 24.3 kPa and N_p= 62. The estimated 3 h MPM is 164.7 kPa, which is larger than 143.6 kPa, as shown in Fig. 2. The fitting results look better than Fig. 2. In fact, considering IQR = 4.172 and γ = 7.57, Fig. 6 has no outlier. Furthermore, the data follows a Fréchet distribution rather than a Weibull distribution because positive shape parameter ξ= 0.434. Because this high threshold reduces the bias and the number of data points but increases the variance, it is necessary to choose a lower threshold as far as possible to increase the amount of data and decrease the variance.

Fig. 7 and Fig. 8 are used to check the IID and goodness of fit, and show somewhat nice fitting with R = 0.993 in Fig. 8.

The MEF plot shown in Fig. 9 indicates that selecting a threshold of 24.5 kPa is not bad, because the MEF shows a linear tendency above the threshold (see section 2.1) except for a high threshold due to the lack of data. Hence, it can be said that a proper threshold is around 24.5 kPa.

Finally, considering an industry practice of 3-parameter Weibull fitting for the largest 10% of the data, the 3 h MPM of 150.7 kPa is obtained with u = 21.0 kPa and N_p= 77, as shown in Fig. 10. This is slightly larger than 143.6 kPa for the 3-parameter Weibull fitting with all 764 data points in Fig. 2. The last two large pressures are considered the outliers, with IRQ = 14.738 and γ = 3.58. Considering ξ > 0 in Fig. 6, the Weibull with MOM is not an appropriate distribution in this case.

4. GEV and GPD for IID Variables

Because the data for GEV and GPD should be IID random variables, the block maxima for GEV and the exceedances over the threshold for GPD are used to check IID using the PACF and MEF in this study, respectively. For MEF in Fig. 9, it is, however, somewhat challenging to interpret for a linear region to choose the threshold and disadvantageous to eliminate the manual work for determining the threshold. Hence, PACF is also used in checking IID for GPD in this case.

At first, for GPD, Fig. 7 and Fig. 11 show ACF and PACF for the largest 8%, 10%, and 11% of the data, respectively. All of the PACFs satisfy the 95% CI,. Hence, all of the cases pass the IID checks. The 11% case has the largest number of data points (N_p= 85) among them, and data larger than 11% cannot pass the IID check. Therefore, u = 19.8 kPa for the 11% case is finally chosen for the GPD threshold.

Fig. 12 shows the GPD fitting with the largest 11% of the data (N_p=85). The estimated 3 h MPM for GPD is 167.14 kPa. This figure also has no outlier with IRQ = 4.504 and γ = 5.85.

Next, for GEV, Fig. 13 and Fig. 14 show the PACFs and Q-Q plots for block sizes of 4, 6, and 10, respectively. From PACF, all of the cases pass the IID checks. However, the Q-Q plots show that the case of a block size of 10 is the most appropriate choice with the largest R = 0.995. Hence, the final block size of this study is 10, with a time window of 2.0 s.

Fig. 15 shows the GEV fitting with 77 tail data points for a block size of 10. The estimated 3 h MPM for GEV is 184.8 kPa. This figure also has no outlier with IRQ = 0.087 and γ = 200.

From the GPD and GEV fitting figures, it can be said that these fitting methods with IID checking may give statistically more appropriate design pressures. The positive GEV shape parameter (ξ>0.489 in Fig. 15) yields the Fréchet distribution.

5. Concluding Remarks

In this study, GEV and GPD were used to fit IID random variables from a sloshing model test, as well as the conventional analysis procedures. In order to apply GEV and GPD, the IID check was carried out in advance for highly nonlinear data by MEF and graphical diagnostic plots. For the feasibility study of the resampling IID variables approach, only one 5 h sample of sloshing test results was used in this study. Hence, further study with various samples with test durations of 5 h or longer (e.g., 30 h) and other filling cases (e.g., 15% to 70% of tank height) may be needed to verify this approach. The following concluding remarks can be made for this study.

(1) In the case of impact pressures obtained by the sloshing model test in this study, the 3-parameter Weibull with MOM may not be an appropriate distribution for all of the sampling data using the conventional procedure with an initial threshold and time window.

(2) The GPD fitting for the largest 8% of the data, as one of the conventional procedures, gives a relatively reasonable design pressure. Because this high threshold reduces the bias and the number of data points but increases the variance, it is necessary to choose a lower threshold as far as possible to increase the number of data points and decrease the variance. In addition, the obtained positive GPD and GEV shape parameters give evidence that the data to be analyzed follows a Fréchet distribution rather than a Weibull distribution.

(3) The 3-parameter Weibull fitting for the largest 10% of the data as an industry practice showed a bad performance because of the positive GPD shape parameter for the sample in this study.

(4) Resampling of the IID variables was carried out using ACF, PACF, MEF, and graphical diagnostic plots for GEV and GPD. From this study, resampling IID variables may reduce the number of outliers, and yield a statistically more appropriate design value because of the independent sample.