Conﬁdence Intervals for Comparing the Variances of Two Independent Birnbaum–Saunders Distributions

: Fatigue in a material occurs when it is subjected to ﬂuctuating stress and strain, which usually results in failure due to the accumulated damage. In statistics, asymmetric distribution, which is commonly used for describing the fatigue life of materials, is the Birnbaum–Saunders (BS) distribution. This distribution can be transform to the normal distribution, which is symmetrical. Furthermore, variance is used to examine the dispersion of the fatigue life data. However, comparing the variances of two independent samples that follow BS distributions has not previously been reported. To accomplish this, we propose methods for providing the conﬁdence interval for the ratio of variances of two independent BS distributions based on the generalized ﬁducial conﬁdence interval (GFCI), a Bayesian credible interval (BCI), and the highest posterior density (HPD) intervals based on a prior distribution with partial information (HPD-PI) and a proper prior with known hyperparameters (HPD-KH). A Monte Carlo simulation study was carried out to examine the efﬁcacies of the methods in terms of their coverage probabilities and average lengths. The simulation results indicate that the HPD-PI performed satisfactorily for all sample sizes investigated. To illustrate the efﬁcacies of the proposed methods with real data, they were also applied to study the conﬁdence interval for the ratio of the variances of two 6061-T6 aluminum coupon fatigue-life datasets.


Introduction
Fatigue, defined as the degradation of the mechanical properties of a material under loading that change over time, is one of the leading causes of machine and structural failure. A critical characteristic of fatigue is that the load is not sufficiently large to cause instantaneous failure. Instead, failure occurs after a particular number of load fluctuations have been encountered (i.e., after the cumulative damage has reached a critical threshold) [1]. As a result, having a good understanding of the fatigue life of materials is critical for preventing damage caused by their failure, predicting the consequences of changes in operational conditions, identifying the cause of fatigue failure, and taking effective mitigating measures. In order to evaluate the fatigue life of materials, statistical distributions of the fatigue life can be considered. These distributions are often positive asymmetry or skewness (non-normality) and start from zero, since the fatigue life is always non-negative. Therefore, the fatigue life of materials cannot be described by either the normal or symmetrical distributions. In recent decades, asymmetric distribution that has received considerable attention for describing the fatigue life of materials is the Birnbaum-Saunders (BS) distribution. It was originally developed in response to a material fatigue problem and has been extensively used in reliability and fatigue research [2]. The BS distribution explains the total amount of time that will pass until a dominant crack develops and grows to a point where the cumulative damage exceeds the threshold and causes tributions. Wongyai and Suwan [25] developed the confidence interval for the ratio of variances of bivariate non-normal distributions by using a kurtosis estimator. Recently, Maneerat et al. [26] presented the HPD interval based on the normal-gamma prior and the method of variance estimates recovery (MOVER) to compute the confidence interval for the ratio of variances of delta-lognormal distributions. Nevertheless, the construction of the confidence interval for the ratio of variances of two independent BS distributions has not yet been reported. Therefore, the goal of the present study is to propose methods for constructing the confidence interval for the ratio of the variances of two BS distributions based on the generalized fiducial confidence interval (GFCI), a Bayesian credible interval (BCI), and the HPD intervals based on a prior with partial information (HPD-PI) and a proper prior with known hyperparameters (HPD-KH).
The rest of this article is structured as follows. The background on the BS distribution and the concepts of each of the methods for constructing the confidence interval for the ratio of variances of two BS distributions are described in Section 2. The simulation studies and results are presented in Section 3. Section 4 provides an illustration of the proposed methods with real fatigue datasets from Birnbaum and Saunders [5]. The final section contains conclusions on the study.

Methods
Let X ij = (X i1 , X i2 , . . . , X in i ), i = 1, 2 and j = 1, 2, . . . , n i be non-negative random samples drawn from BS distributions denoted by X ij ∼BS(α i , β i ), where α i and β i are the shape and scale parameters, respectively. The cumulative distribution function (cdf) can be written as where x ij > 0, α i > 0, β i > 0, and Φ(·) is the standard normal cdf. Note that β i is also the median of the distribution. The corresponding probability density function (pdf) of this cdf is given by The following transformation was applied to generate samples from the BS distributions and to enable the derivation of some of their other properties, including various moments. If X ij ∼BS(α i , β i ), then Thus, By applying the above transformation, the expected value and variance of X ij are E(X ij ) = β i 1 + 1 2 α 2 i and V(X ij ) = (α i β i ) 2 1 + 5 4 α 2 i , respectively. Since X ij are independent, the ratio of the variances simply becomes

Generalized Fiducial Inference
Generalized fiducial inference can be used to transform the original data into other distributions that are known. According to the rules of that distribution, the transformed data are manipulated, and the results are transferred back to the original via an inverse transformation [27]. This idea brings us to construct the confidence interval for the ratio of the variances of two BS distributions. Let Y = G(η, U) (6) be the relationship between Y and parameter η ∈ Ξ, where G(·, ·) is a structural equation and U is a random variable for which the distribution is definitively known and independent of any parameters. For any given realization y of Y, inverse η = H(y, u) always exists for any realization u of U. Since the distribution of U is definitively known, random samplẽ u 1 ,ũ 2 , . . . ,ũ M can be generated from it. This random sample of U can be transformed into a random sample of η via the inverseη 1 = H(y,ũ 1 ),η 2 = H(y,ũ 2 ), . . . ,η M = H(y,ũ M ), such that a random sample of η (i.e., a fiducial sample) can be obtained. However, in some situations, the inverse does not exist, for which Hannig [27,28] proposed the following solutions.
For a BS distribution, the generalized fiducial distribution of (α i , β i ) derived by Li and Xu [17] is in the form where and Note that the symbol "∝" means "is proportional to." In brief, if a is proportional to b, then the only difference between a and b is a multiplicative constant. By applying Equation (11), Li and Xu [17] showed that the priors of α i and β i can be denoted as Thus, f (α i , β i |x ij ) is proper for the particular case of a prior with partial information given by the priors of α i and β i (12). Therefore,α i andβ i , which are the generalized fiducial samples of α i and β i , can be obtained from the generalized fiducial distribution in the same way as the Bayesian posterior. The adaptive rejection Metropolis sampling (ARMS), which origins to adaptive rejection sampling (ARS), was used to generate the fiducial samples (α i andβ i ) from the generalized fiducial distribution (9). The ARS was proposed by Gilks and Wild [29]. It was only suitable for log-concave target densities. In order to address the limitations of ARS, Gilks et al. [30] improved ARS to handle multivariate distributions and non-log-concave densities by permitting the proposal distribution to remain lower than the target in some regions and adding a Metropolis-Hastings step to guarantee that the accepted samples are properly distributed. This method was called ARMS, which can be easily implemented via the function arms in package dlm of R software suite (version 3.5.1). Note thatα i andβ i are treated as random variables. Therefore, α i and β i are substituted bŷ α i andβ i , respectively, resulting in the generalized fiducial estimates of θ being derived aŝ Finally, the 100 is the 100ν% percentile ofθ. Algorithm 1 summarizes the steps for constructing GFCI for θ, as seen below.

2.
Generate K samples of α i and β i by applying the arms function in the dlm package of the R software suite.

3.
Burn-in B samples (the number of remaining samples is K − B).

4.
Thin the samples by applying sampling lag L > 1 (the final number of samples is K = (K − B)/L). Note that the generated samples are not independent, and so we need to reduce the autocorrelation by thinning the samples.

Bayesian Inference
For this method, Xu and Tang [14] illustrated that the reference prior of a BS distribution (a type of Jeffreys' prior) results in an improper posterior distribution. Thus, to guarantee its propriety, proper priors with known hyperparameters are obtained by assuming that an inverse-gamma distribution with parameters a i and b i is the prior of β i and an inverse-gamma distribution with parameters c i and d i is the prior of λ i = α 2 i [16]. In accordance with Bayes' theorem, the joint posterior density function of (α i , β i ) can be written as Integrating the joint posterior density function (14) with respect to α i yields the marginal posterior distribution of β i as follows: From the joint posterior density function (14), it is clear that the fully conditional posterior distribution of λ i given β i is given by Posterior samples are drawn by adopting the Markov chain Monte Carlo technique. Since the marginal posterior distribution (15) cannot be written as if it were known, generating posterior samples of β i from this density is impossible using the usual methods. There are three common approaches, such as the random-walk Metropolis procedure, the Metropolis-Hastings algorithm and the slice sampler by introducing an auxiliary variable to simplify the sampling problem that might be considered to sample from the marginal posterior distribution (15). However, all three approaches are susceptible to serially correlated draws, indicating that a very large sample size is frequently required to produce a reasonable estimate of any desired attribute of the posterior distribution. To avoid these potential problems when generating the posterior samples, the generalized ratio-of-uniforms method of Wakefield et al. [31] is used to generate posterior samples of β i (denoted asβ i ). The concept of the generalized ratio-of-uniforms method is as follows.
Suppose that a pair of variables (u i , v i ) is uniformly distributed inside region where r i ≥ 0 is a constant term and π(·|x ij ) is specified by using the marginal posterior distribution (15). Subsequently, the pdf of β i = v i /u r i i becomes π(β i |x ij )/ π(β i |x ij )dβ i . For generating random points uniformly distributed in A(r i ), the accept-reject method from a convenient enveloping region (usually from the minimal bounding rectangle) is applied. According to Wakefield et al. [31], the minimal bounding rectangle for A(r i ) is given by where a(r i ) = sup Note that π(β i |x ij ) → 0 as β i → 0 + and π(β i |x ij ) → O(β −(a i +c i +3/2) i ) as β i → +∞. Hence, b − (r i ) = 0, a(r i ) is finite, and b + (r i ) is also finite when choosing an appropriate value for r i [16]. The generalized ratio-of-uniforms method consists of the following three steps.

1.
Compute a(r i ) and b + (r i ).

2.
Draw u i and v i from U(0, a(r i )) and U(0, b + (r i )), where U(v, w) refers to a uniform distribution with parameter v and w, and compute ρ i = v i /u r i i .
Meanwhile,λ i , which are the posterior samples of λ i , can be obtained from the conditional posterior distribution (16) by applying the LearnBayes package from the R software suite. Subsequently, the posterior samples of α i (denoted asα i ) comprise the square roots ofλ i . Note thatα i andβ i are also treated as random variables. Hence, the Bayesian estimates for θ can be written as Finally, the 100(1 − γ)% BCI for θ is [θ(γ/2),θ(1 − γ/2)], whereθ(ν) is the 100ν% percentile ofθ. In conclusion, BCI for θ can be obtained by using Algorithm 2.

1.
Set a i , b i , c i and d i , where i = 1, 2.

2.
Compute a(r i ) and b + (r i ).

3.
At the m th steps, (a) Generate u i ∼U(0, a(r i )) and v i ∼U(0, b + (r i )), and then compute Compute the Bayesian estimates for θ by applying Equation (22).

The Highest Posterior Density (HPD) Interval
The HPD interval is where the posterior density for every point within the interval is higher than the posterior densities of the points outside of it, indicating that the interval contains the more likely values of the parameter while excluding the less likely ones. According to Box and Tiao [32], the HPD interval has two main properties:

1.
Every point within the interval has a higher probability density than the points outside of it.
By applying Equation (11), Li and Xu [17] showed that J(x ij , (α i , β i )) is a special case of a prior with partial information, and the generalized fiducial estimates of α i and β i can be obtained by using the same method as for the Bayesian posterior. Therefore, at Step (6) in Algorithm 1, we applied the HDInterval package (version 0.2.2) from the R software suite to compute the HPD interval based on a prior with partial information (HPD-PI). Moreover, we also applied the HDInterval package at Step (5) in Algorithm 2 to compute the HPD interval based on a proper prior with known hyperparameters (HPD-KH).
Without loss of generality, scale parameters β 1 and β 2 were set as 1 in all scenarios. The confidence intervals were calculated at the nominal level of 0.95. All simulation results were obtained by running 1000 replications with K = 3000 and B = 1000 for GFCI and HPD-PI while M = 1000 for BCI and HPD-KH. According to Wang et al. [16], BCI and HPD-KH were considered with r 1 = r 2 = 2 and hyperparameters a i = b i = c i = d i = 10 −4 . The criteria for evaluating the performances of the proposed methods are their coverage probabilities and average lengths. The method with a coverage probability greater than or close to the nominal level 0.95 and with the narrowest average length was chosen as the best performing method for a particular scenario.
Tables 1 and 2 report the simulation results, while Figures 1 and 2 summarize the coverage probabilities and average lengths in Tables 1 and 2. The simulation results from Tables 1 and 2 indicate that the coverage probabilities of the four methods were greater than or close to 0.95 under all configurations. In addition, it was found that the differences in coverage probability among the four methods were very small. Moreover, HPD-PI provided the narrowest average lengths while BCI provided the longest ones under all circumstances. The average lengths of HPD-IP were mostly narrower than GFCI, while the average lengths of HPD-KH were mostly narrower than those of BCI. Moreover, the average length of the four methods decreased as the sample sizes (n 1 , n 2 ) increased. Table 1. Coverage probabilities and average lengths of the 95% confidence interval for the ratio of variances of two BS distributions with equal sample sizes (n 1 = n 2 ) constructed via the various methods.

Coverage Probability
Average Length (n 1 , n 2 ) (α 1 , α 2 ) GFCI BCI HPD-PI HPD-KH GFCI BCI HPD-PI HPD-KH  GFCI, the generalized fiducial confidence interval; BCI, the Bayesian credible interval; HPD-PI, the highest posterior density interval based on a prior with partial information; HPD-KH, the highest posterior density interval based on the proper prior with known hyperparameters.

Application of the Methods to Real Fatigue Life Data
To illustrate the effectiveness of the confidence interval construction methods proposed in this study in a real-life scenario, we used real datasets concerning the fatigue life of The results for the 95% confidence interval for the ratio of variances reported in Table 5 indicate that the length provided by HPD-PI was the narrowest while that of the Bayesian credible interval was the longest. These results are in accordance with those from the simulation studies where (n 1 , n 2 ) = (100, 100). In addition, the confidence intervals constructed by using the various methods did not contain 1, and so it can be concluded that there is no significant difference in terms of variance for the fatigue lifetime of 6061-T6 aluminum coupons with maximum stress per cycle of 21,000 and 26,000 psi, respectively. GFCI, the generalized fiducial confidence interval; BCI, the Bayesian credible interval; HPD-PI, the highest posterior density interval based on a prior with partial information; HPD-KH, the highest posterior density interval based on the proper prior with known hyperparameters.

Conclusions
Four methods, namely GFCI, BCI, HPD-PI, and HPD-KH, were proposed for constructing the confidence interval for the ratio of variances of two BS distributions. A Monte Carlo simulation study was conducted to assess their performances based on their coverage probabilities and average lengths. The simulation study results indicate that the coverage probabilities of all of the methods were greater than or close to the nominal level of 0.95. However, HPD-PI outperformed the others by providing the narrowest average lengths in all of the scenarios studied. In addition, the results of using fatigue lifetime data of 6061-T6 aluminum coupons coincided with those from the simulation study. Therefore, HPD-PI can be recommended for constructing the confidence interval for the ratio of variances of two BS distributions.

Data Availability Statement:
The real data sets of the fatigue life of 6061-T6 aluminum coupons are available in the work of Birnbaum and Saunder [5].