Confidence intervals for the difference between the coefficients of variation of Weibull distributions for analyzing wind speed dispersion

The Generalized Confidence Interval

The generalized confidence interval (GCI) was introduced by Weerahandi (1993). The most important concept of GCI is the generalized pivotal quantity (GPQ) defined as

Definition Suppose that X = (X₁, X₂, …, X_n) be a random variables from distribution, which depends on a parameter of interest φ, and a nuisance parameter γ. Furthermore, suppose that x =(x₁, x₂, …, x_n) be the observed value of X. To obtain a GCI for φ, a GPQ R(X; x, φ, γ) is also required to satisfy the following two conditions.

(i) For a fixed x, the distribution of R(X; x, φ, γ) is free of unknown parameters.

(ii) The observed value of R(X; x, φ, γ) at X = x, denoted as r(X; x, φ, γ) does not depend on nuisance parameters.

Let R_φ(α∕2) and R_φ(1 − α∕2) are the 100(α/2)-th and 100(1 − α/2)-th percentile of R_φ(X; x), respectively. The 100 (1 − α)% two-sided confidence interval based on GCI for parameter of interest is given by [R_φ(α∕2), R_φ(1 − α∕2)]. Note that 1 − α is the probability that population parameter will be in the interval. The 100 (1 − α)% confidence interval will include the true value of the population parameter with probability 1 − α, i.e., if α = 0.05, that probability is about 0.95 that the 95% confidence interval will include the true population parameter.

For Weibull distribution, Thoman, Bain & Antle (1969) showed that the distributions of $\frac{\hat{k}}{k}$ and $\hat{k}ln\left(\frac{\hat{a}}{a}\right)$ do not depend on a and k, where $\hat{a}$ and $\hat{k}$ are the MLEs of a and k, respectively. Let $\widehat{a^{\ast }}$ and $\widehat{k^{\ast }}$ be the MLEs based on a sample of size n from a Weibull(1, 1). Then $\frac{\hat{k}}{k}\sim \widehat{k^{\ast }}$ and $\hat{k}ln\left(\frac{\hat{a}}{a}\right)\sim \widehat{k^{\ast }}ln\left(\widehat{a^{\ast }}\right)$. Both distributions do not depend on any parameters. Hence, they are pivotal quantity.

Krishnamoorthy, Mukherjee & Guo (2007) presented the GPQs of scale and shape parameters from Weibull distribution. Let ${\hat{a}}_0$ and ${\hat{k}}_0$ be the observed values of $\hat{a}$ and $\hat{k}$. The GPQs of the parameters are given by (8)${\displaystyle R_k=\frac{k}{\hat{k}}{\hat{k}}_0=\frac{{\hat{k}}_0}{{\hat{k}}^{\ast }}}$ and (9)${\displaystyle R_a={\left(\frac{a}{\hat{a}}\right)}^{\frac{\hat{k}}{{\hat{k}}_0}}{\hat{a}}_0={\left(\frac{1}{{\hat{a}}^{\ast }}\right)}^{\frac{{\hat{k}}^{\ast }}{{\hat{k}}_0}}{\hat{a}}_0.}$ From the random variables X and Y, the GPQs for k_X and k_Y are given by (10)${\displaystyle R_{k_X}=\frac{{\hat{k}}_{X_0}}{{\hat{k}}_X^{\ast }}}$ and (11)${\displaystyle R_{k_Y}=\frac{{\hat{k}}_{Y_0}}{{\hat{k}}_Y^{\ast }}.}$ Thus, the GPQ for δ is defined by (12)${\displaystyle R_{\delta }=R_{{\lambda }_X}-R_{{\lambda }_Y}=\sqrt{\frac{\Gamma \left(1+\frac{2}{R_{k_X}}\right)}{{\left(\Gamma \left(1+\frac{1}{R_{k_X}}\right)\right)}^2}-1}-\sqrt{\frac{\Gamma \left(1+\frac{2}{R_{k_Y}}\right)}{{\left(\Gamma \left(1+\frac{1}{R_{k_Y}}\right)\right)}^2}-1}.}$

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on GCI is (13)${\displaystyle C{I}_{\delta \left(gci\right)}=\left[L_{\delta \left(gci\right)},U_{\delta \left(gci\right)}\right]=\left[R_{\delta }\left(\alpha ∕2\right),R_{\delta }\left(1-\alpha ∕2\right)\right],}$ where R_δ(α∕2) denotes the 100(α/2)-th percentile of R_δ.

The following algorithm is used to construct confidence interval based on GCI for the difference between the coefficients of variation of Weibull distributions:

Algorithm 2

1. Compute ${\hat{a}}_X$, ${\hat{a}}_Y$, ${\hat{k}}_X$ and ${\hat{k}}_Y$ from Algorithm 1.

2. Generate $X_1^{\ast }$, $X_2^{\ast }$, …, $X_n^{\ast }$ and $Y_1^{\ast }$, $Y_2^{\ast }$, …, $Y_m^{\ast }$ from Weibull(1, 1).

3. Compute ${\hat{a}}_X^{\ast }$, ${\hat{a}}_Y^{\ast }$, ${\hat{k}}_X^{\ast }$ and ${\hat{k}}_Y^{\ast }$ from Algorithm 1.

4. Compute R_kX and R_kY from Eqs. (10) and (11).

5. Compute R_δ from Eq. (12).

6. Repeat steps 1-5 for q times, where q is the number of generalized computation.

7. Compute 95% confidence interval based on GCI, as given in Eq. (13).

The bootstrap confidence intervals

The bootstrap method was introduced by Efron & Tibshirani (1993). This is a resampling method to determine precision measures for statistical estimation.

Let ${\mathbf{x}}_i=\left(x_1,x_2,\dots ,x_n\right)$ and ${\mathbf{y}}_i=\left(y_1,y_2,\dots ,y_m\right)$ be random samples of size n and m from Weibull distributions, and let ${\mathbf{x}}_i^{\ast b}=\left(x_1^{\ast b},x_2^{\ast b},\dots ,x_n^{\ast b}\right)$ and ${\mathbf{y}}_i^{\ast b}=\left(y_1^{\ast b},y_2^{\ast b},\dots ,y_m^{\ast b}\right)$ be bootstrapped samples drawn with replacement from the original data, using the same sample sizes. After resampling B bootstrap samples, the difference between the coefficients of variation are calculated in each bootstrap sample, as follows: ${\delta }^{\ast b}={\lambda }_X^{\ast b}-{\lambda }_Y^{\ast b},$b =1 , 2, …, B.

The percentile bootstrap confidence interval

The percentile bootstrap confidence interval is based on the percentile of the distribution of the bootstrapped replications. The value of the bootstrap statistic, δ^∗b are ordered from the smallest to the largest.

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on percentile bootstrap is given by (14)${\displaystyle C{I}_{\delta \left(pb\right)}=\left[L_{\delta \left(pb\right)},U_{\delta \left(pb\right)}\right]=\left[{\delta }^{\ast b}\left(\alpha ∕2\right),{\delta }^{\ast b}\left(1-\alpha ∕2\right)\right],}$ where δ^∗b(α∕2) denotes the 100(α/2)-th percentile of δ^∗b.

The following algorithm is used to construct confidence interval based on percentile bootstrap for the difference between the coefficients of variation of Weibull distributions:

Algorithm 3

1. Generate X₁, X₂, …, X_n from Weibull(a_X, k_X) and Y₁, Y₂, …, Y_m from Weibull(a_Y, k_Y).

2. For b = 1.

3. Resampling bootstrap samples $X_1^{\ast },X_2^{\ast },\dots ,X_n^{\ast }$ from X₁, X₂, …, X_n and compute ${\lambda }_X^{\ast b}$.

4. Resampling bootstrap samples $Y_1^{\ast },Y_2^{\ast },\dots ,Y_m^{\ast }$ from Y₁, Y₂, …, Y_m and compute ${\lambda }_Y^{\ast b}$.

5. Compute ${\delta }^{\ast b}={\lambda }_X^{\ast b}-{\lambda }_Y^{\ast b}$.

6. Repeat steps 3-5 for B times, where B is the number of bootstrap sample.

7. Sort δ^∗b from the smallest to the largest.

8. Compute 95% confidence interval based on percentile bootstrap, as given in Eq. (14).

The bootstrap confidence interval with standard errors

From the B bootstrap statistic denoted as δ^∗b, b =1 , 2, …, B, we can calculate the standard error of a statistic. They can be estimated using the standard deviation of the bootstrap distribution.

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on bootstrap confidence interval with standard errors can be written as (15)${\displaystyle C{I}_{\delta \left(bs\right)}=\left[L_{\delta \left(bs\right)},U_{\delta \left(bs\right)}\right]=\left[\hat{\delta }-Z_{\left(\alpha ∕2\right)}SE,\hat{\delta }+Z_{\left(\alpha ∕2\right)}SE\right],}$ where SE is the standard error of a statistic.

The following algorithm is used to construct confidence interval based on bootstrap confidence interval with standard errors for the difference between the coefficients of variation of Weibull distributions:

Algorithm 4

1. Generate X₁, X₂, …, X_n from Weibull(a_X, k_X) and Y₁, Y₂, …, Y_m from Weibull(a_Y, k_Y).

2. For b = 1.

3. Resampling bootstrap samples $X_1^{\ast },X_2^{\ast },\dots ,X_n^{\ast }$ from X₁, X₂, …, X_n and compute ${\lambda }_X^{\ast b}$.

4. Resampling bootstrap samples $Y_1^{\ast },Y_2^{\ast },\dots ,Y_m^{\ast }$ from Y₁, Y₂, …, Y_m and compute ${\lambda }_Y^{\ast b}$.

5. Compute ${\delta }^{\ast b}={\lambda }_X^{\ast b}-{\lambda }_Y^{\ast b}$.

6. Repeat steps 3-5 for B times, where B is the number of bootstrap sample.

7. Compute SE.

8. Compute 95% confidence interval based on bootstrap confidence interval with standard errors, as given in Eq. (15).

The method of variance estimates recovery

According to Donner & Zou (2010), this approach can be used to construct a confidence interval for a function of two parameters, λ_X − λ_Y, defined as (16)${\displaystyle C{I}_m=\left[L_m,U_m\right],}$ where the lower limit and upper limit for ${\hat{\lambda }}_X$ - ${\hat{\lambda }}_Y$ are given by (17)${\displaystyle L_m=\left({\hat{\lambda }}_X-{\hat{\lambda }}_Y\right)-\sqrt{{\left({\hat{\lambda }}_X-l_X\right)}^2+{\left(u_Y-{\hat{\lambda }}_Y\right)}^2}}$ (18)${\displaystyle U_m=\left({\hat{\lambda }}_X-{\hat{\lambda }}_Y\right)+\sqrt{{\left(u_X-{\hat{\lambda }}_X\right)}^2+{\left({\hat{\lambda }}_Y-l_Y\right)}^2}.}$ Hendricks & Robey (1936) presented confidence intervals for λ_X and λ_Y defined as (19)${\displaystyle \left(l_{X.HR},u_{X.HR}\right)=\left({\hat{\lambda }}_X-t_{\left(\alpha ∕2,n-1\right)}\frac{{\hat{\lambda }}_X}{\sqrt{2n}},{\hat{\lambda }}_X+t_{\left(\alpha ∕2,n-1\right)}\frac{{\hat{\lambda }}_X}{\sqrt{2n}}\right)}$ and (20)${\displaystyle \left(l_{Y.HR},u_{Y.HR}\right)=\left({\hat{\lambda }}_Y-t_{\left(\alpha ∕2,m-1\right)}\frac{{\hat{\lambda }}_Y}{\sqrt{2m}},{\hat{\lambda }}_Y+t_{\left(\alpha ∕2,m-1\right)}\frac{{\hat{\lambda }}_Y}{\sqrt{2m}}\right),}$ where t_{(α∕2,n−1)} and t_{(α∕2,m−1)} denote the 100(α/2)-th percentile of t-distribution with n − 1 and m − 1 degrees of freedom, respectively. To construct a confidence interval for the difference between the coefficients of variation of Weibull distributions based on MOVER with Hendricks and Robey’s confidence interval, we substitute l_X.HR, u_X.HR, l_Y.HR and u_Y.HR for l_X, u_X, l_Y and u_Y, respectively.

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on MOVER with Hendricks and Robey’s confidence interval becomes (21)${\displaystyle C{I}_{\delta \left(m.HR\right)}=\left[L_{\delta \left(m.HR\right)},U_{\delta \left(m.HR\right)}\right],}$ where (22)${\displaystyle L_{\delta \left(m.HR\right)}=\left({\hat{\lambda }}_X-{\hat{\lambda }}_Y\right)-\sqrt{{\left({\hat{\lambda }}_X-l_{X.HR}\right)}^2+{\left(u_{Y.HR}-{\hat{\lambda }}_Y\right)}^2}}$ (23)${\displaystyle U_{\delta \left(m.HR\right)}=\left({\hat{\lambda }}_X-{\hat{\lambda }}_Y\right)+\sqrt{{\left(u_{X.HR}-{\hat{\lambda }}_X\right)}^2+{\left({\hat{\lambda }}_Y-l_{Y.HR}\right)}^2}.}$

The following algorithm is used to construct confidence interval based on MOVER for the difference between the coefficients of variation of Weibull distributions:

Algorithm 5

1. Generate X₁, X₂, …, X_n from Weibull(a_X, k_X) and Y₁, Y₂, …, Y_m from Weibull(a_Y, k_Y).

2. Compute the intervals for λ_X and λ_Y from Eqs. (19) and (20).

3. Compute 95% confidence interval based on MOVER, as given in Eq. (21).

The Bayesian confidence intervals

Here, we derive Bayesian estimates of the parameters of a Weibull distribution. Once again, for random variable X = (X₁, X₂, …, X_n) from a Weibull distribution, we transform scale parameter $a^{\prime }={\left(\frac{1}{a}\right)}^k$, and so we can provide its probability density function in another form as follows: (24)${\displaystyle f\left(x;a^{\prime },k\right)=a^{\prime }k{x}^{k-1}exp\left[-a^{\prime }{x}^k\right],x>0.}$

It is assumed that the prior distributions for the scale and shape parameters are gamma with hyperparameters (v₁, v₂, z₁, z₂): (25)${\displaystyle \pi \left(a^{\prime }\right)\sim gamma\left(v_1,z_1\right)}$ and (26)${\displaystyle \pi \left(k\right)\sim gamma\left(v_2,z_2\right).}$ Hence, using Bayes’ theorem, the joint posterior density function of a′ and k is given by (27)${\displaystyle \pi \left(a^{\prime },k|x\right)\propto L\left(a^{\prime },k|x\right)\pi \left(a^{\prime }\right)\pi \left(k\right),}$ where L(a′, k|x) is a likelihood function.

However, the posterior distributions cannot be computed, and so we used the MCMC method with Gibbs sampling to provide them for the parameters. We used the Gibbs sampling procedure to generate samples from the joint posterior distribution, as given in Eq. (27). The conditional posterior distributions of parameters are (28)${\displaystyle \pi \left(a^{\prime }|k,x\right)\sim gamma\left(v_2+n,z_2+\sum x^k\right)}$ and (29)${\displaystyle \pi \left(k|a^{\prime },x\right)\propto k^{v_1+n-1}exp\left[-v_{1}k-a^{\prime }\sum x^k\right].}$ For Eq. (29), we used the Metropolis-Hasting (MH) algorithm to update shape parameter k. Moreover, the random walk Metropolis (RWM) and Gibbs sampling methods are applied as follows (Geman & Geman, 1984):

Algorithm 6 The Gibbs algorithm

1. Take the initial value of parameter (a′⁽⁰⁾, k⁽⁰⁾).

2. Generate a′^(t) ∼ gamma(v₂ + n, z₂ + ∑x^k(t−1)).

3. Using Random walk Metropolis (RWM) algorithm for calculate k^(t).

4. Repeat step 2–3 for T times, where T is the number of replications of MCMC.

5. Burn in 1000 samples and compute the parameter of interest.

Algorithm 7 The random walk Metropolis (RWM)

1. Start with (a′^(t), k^(t−1)).

2. Generate $\epsilon \sim N\left(0,{\sigma }_k^2\right)$.

3. Compute k^∗ = k^(t−1) + ε.

4. Compute $A_k=\frac{L\left(k^{\ast },a^{\prime }|x\right)\pi \left(k^{\ast }\right)}{L\left(k,a^{\prime }|x\right)\pi \left(k\right)}$.

5. Generate u ∼ U(0, 1).

6. If u ≤ min(1, A_k) set k^(t) = k^∗ and If u > min(1, A_k) set k^(t) = k^(t−1).

The Bayesian-MCMC

Once again, let X and Y be random variables from Weibull distributions. First, we used Algorithms 6 and 7 to calculate the difference between the coefficients of variation denoted as δ^t, t =1 , 2, …, T, based on the Bayesian-MCMC.

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on the Bayesian-MCMC is given by (30)${\displaystyle C{I}_{\delta \left(MCMC\right)}=\left[L_{\delta \left(MCMC\right)},U_{\delta \left(MCMC\right)}\right]=\left[{\delta }^t\left(\alpha ∕2\right),{\delta }^t\left(1-\alpha ∕2\right)\right],}$ where ${\delta }^t\left(\alpha ∕2\right)$ denotes the 100(α/2)-th percentile of δ^t.

The following algorithm is used to construct confidence interval based on Bayesian-MCMC for the difference between the coefficients of variation of Weibull distributions:

Algorithm 8

1. Calculate the initial value of parameter (a′⁽⁰⁾, k⁽⁰⁾).

2. Calculate δ^t from Algorithms 6 and 7.

3. Compute 95% confidence interval based on Bayesian-MCMC, as given in Eq. (30).

The Bayesian-Highest Posterior Density (HPD) Interval

Here, the assumption is that the density of every point inside the interval is greater than that of every point outside the interval, and it is also the shortest interval (Box & Tiao, 2011). In this article, the HPD interval was computed by using the R package with HDInterval.

Therefore, the 100 (1 − α)% two-sided confidence interval for the difference between the coefficients of variation of Weibull distributions based on the Bayesian-HPD interval can be written as (31)${\displaystyle C{I}_{\delta \left(HPD\right)}=\left[L_{\delta \left(HPD\right)},U_{\delta \left(HPD\right)}\right].}$

The following algorithm is used to construct confidence interval based on Bayesian-HPD for the difference between the coefficients of variation of Weibull distributions:

Algorithm 9

1. Calculate the initial value of parameter (a′⁽⁰⁾, k⁽⁰⁾).

2. Calculate δ^t from Algorithms 6 and 7.

3. Compute 95% confidence interval based on Bayesian-HPD interval, as given in Eq. (31).

Example 1

Wind speed data set were collected at 90-meter wind energy potential stations in Trad and Chonburi provinces in 2016 by the Department of Alternative Energy Development and Efficiency, Ministry of Energy (https://www.dede.go.th/ewt_news.php?nid=47706). These are coastal provinces in eastern Thailand that have been identified as potential areas for harvesting wind energy. Figure 3 shows Q-Q plots of the Weibull distributions of the data from the two provinces. Furthermore, the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values reported in Table 3 indicate that these two datasets fit Weibull distributions because the latter have the smallest values for both criteria.

The statistical summary of the estimates of the parameters for these datasets is (n, m) = (12,11), ${\hat{a}}_X=1.7621$, ${\hat{a}}_Y=4.5839$, ${\hat{k}}_X=1.4832$, ${\hat{k}}_Y=9.0794$, ${\hat{\lambda }}_X=0.6860$ and ${\hat{\lambda }}_Y=0.1317$, while the actual difference between the coefficients of variation of the two Weibull distributions $\hat{\delta }=0.5543$. The confidence intervals based on the proposed methods are given in Table 4. The results show that the Bayesian-HPD interval performed well in terms of the coverage probability and expected length when the sample sizes are small, which supports the simulation results. Finally, Fig. 4 shows a trace plot of the generated δ value.

Example 2

Data on wind speeds measured at 90-meter wind energy potential stations in the southern and northeastern regions of Thailand were collected in April–May, 2019 by the Department of Alternative Energy Development and Efficiency, Ministry of Energy (https://www.dede.go.th/more_news.php?cid=501). The summary statistics for the southern region are n = 12, ${\hat{a}}_X=2.5334$, ${\hat{k}}_X=2.3956$ and ${\hat{\lambda }}_X=0.4441$ and the northeastern region are m = 20, ${\hat{a}}_Y=4.4608$, ${\hat{k}}_Y=3.8457$ and ${\hat{\lambda }}_Y=0.2906$. The actual difference between the coefficients of variation of the Weibull distributions of these datasets $\hat{\delta }=0.1538$. Weibull Q-Q plots of these data are shown in Fig. 5 and the correctness of fitting the data to Weibull distributions in terms of the smallest AIC and BIC values are reported in Table 5. The 95% confidence intervals for δ by applying the six methods (Table 6) indicate that the Bayesian-HPD interval performed well in terms of the coverage probability and expected length when the difference between the coefficients of variation was small, which is once again consistent with the simulation study results. Finally, the trace plot of the generated δ value are presented in Fig. 6.