Bayesian Estimation for the Coefﬁcients of Variation of Birnbaum–Saunders Distributions

: The Birnbaum–Saunders (BS) distribution, which is asymmetric with non-negative support, can be transformed to a normal distribution, which is symmetric. Therefore, the BS distribution is useful for describing data comprising values greater than zero. The coefﬁcient of variation (CV), which is an important descriptive statistic for explaining variation within a dataset, has not previously been used for statistical inference on a BS distribution. The aim of this study is to present four methods for constructing conﬁdence intervals for the CV, and the difference between the CVs of BS distributions. The proposed methods are based on the generalized conﬁdence interval (GCI), a bootstrapped conﬁdence interval (BCI), a Bayesian credible interval (BayCI), and the highest posterior density (HPD) interval. A Monte Carlo simulation study was conducted to evaluate their performances in terms of coverage probability and average length. The results indicate that the HPD interval was the best-performing method overall. PM 2.5 concentration data for Chiang Mai, Thailand, collected in March and April 2019, were used to illustrate the efﬁcacies of the proposed methods, the results of which were in good agreement with the simulation study ﬁndings.


Introduction
When several random variables comprise non-negative values, their dissemination will fit neither a normal nor other symmetrical distributions, and so asymmetrical distributions must be considered instead. Of these, the Birnbaum-Saunders (BS) distribution is receiving considerable attention because of its useful properties and close relationship to the normal distribution. The BS distribution was originally developed to model failure due to the development and growth of cracks in a material subjected to repeated stress cycles [1]. It has two positive parameters: α, the shape parameter, and β, which represents both the scale parameter and the median of the BS distribution. Suppose that random variable X follows a BS distribution with parameters α and β (denoted as X∼BS(α, β)), then its cumulative distribution function (cdf) is given by where Φ(·) is the standard normal cdf. Subsequently, the probability density function (pdf) of the BS distribution can be written as The BS distribution has several interesting properties. For example, for X∼BS(α, β), (i) α −1 ( X/β − β/X)∼N(0, 1), where N(0, 1) denotes a standard normal distribution with exact generalized p-value approach and the delta method, to compare the mean, quantile, and reliability functions of several BS distributions. Li and Xu [26] considered two fiducial methods to estimate the parameters of BS distributions, based on the inverse of the structural equation and Hannig's method. Recently, Guo et al. [27] presented confidence intervals and hypothesis testing for the common mean of several BS distributions. Despite the diverse theoretical and methodological developments for constructing confidence intervals from the functions of parameters of BS distributions, there have not yet been any studies on the single CV and the difference between the CVs of BS distributions. To fill this gap, we propose confidence intervals for these two scenarios by applying the concepts of GCI, the bootstrap confidence interval (BCI), the Bayesian credible interval (BayCI), and the highest posterior density (HPD) interval. Moreover, we applied the proposed methods to datasets of PM2.5 (particulate matter 2.5 µm) concentration in Chiang Mai, Thailand, collected in March and April 2019, to illustrate their efficacies. The rest of this paper is organized as follows. The methodology for the construction of confidence intervals for the CV and the difference between the CVs of BS distributions are described in Sections 2 and 3, respectively. A simulation study and results are presented in Section 4. In Section 5, real datasets are used to illustrate the efficacies of the proposed confidence intervals. Finally, conclusions are provided in Section 6.

The Confidence Interval for the CV of a BS Distribution
Suppose X = (X 1 , X 2 , ..., X n ) is a vector of random samples from a BS distribution with parameters α and β (denoted as X∼BS(α, β)). The cdf of the BS distribution, as given in Equation (1), provides where N(µ, σ 2 ) refers to a normal distribution with mean µ and variance σ 2 . Hence, Therefore, the BS random variable X can be generated from the normal random variable based on this relationship. The expected value and variance of X are given by E(X) = β(1 + 1 2 α 2 ) and Var(X) = (αβ) 2 (1 + 5 4 α 2 ), respectively. Hence, the CV denoted by τ can be defined as

Generalized Confidence Interval
Weerahandi [28] proposed GCI based on the concept of the generalized pivotal quantity (GPQ). This is a generalization of the usual pivot quantity but different in that GPQ is a function of the nuisance parameter, whereas the usual pivotal quantity can only be a function of the sample and the parameter of interest. GPQ has an unknown parameter that is distribution-free and an observed value that does not depend on the nuisance parameter.
Suppose X = (X 1 , X 2 , ..., X n ) is a random sample from the BS distribution in Equation (1) with sample size n, then After denotingZ = n −1 ∑ n i=1 Z i and S 2 = (n − 1) −1 ∑ n i=1 (Z i −Z) 2 , it follows thatZ and S 2 are respectively independently distributed asZ∼N(0, α 2 /n) and (n − 1)S 2 /α 2 ∼χ 2 (n − 1), where χ 2 (n − 1) denotes a Chi-squared distribution with n − 1 degrees of freedom. Hence, R(X; β) = √ n ×Z S follows a t-distribution with n − 1 degrees of freedom that is free from unknown parameter β. Moreover, R(X; β) is only a function of β and is not related to α. Sun [29] proved that R(X; β) is a strictly decreasing function of β, thereby leading to a unique solution of the equation R(X; β) = T as follows: where T∼t(n − 1) refers to a t-distribution with n − 1 degrees of freedom and β 1 and β 2 are the two solutions of quadratic equation where Therefore, R β is a GPQ for β. Subsequently, the GPQ of α [24] is obtained as where is defined in Equation (7). By applying Equation (5), the pivotal quantity for τ can be expressed by Therefore, the 100(1 − γ)% confidence interval for τ based on GCI is where R τ (ν) denotes the 100ν% percentile of R τ .

Bootstrap Confidence Interval
The bootstrap method was originally introduced by Efron [30] as a re-sampling technique based on the random selection of new samples from the original sample. It uses information from the sample to infer certain characteristics from the distribution of a particular statistical estimator. Since τ is a function of α, an estimator of α can be considered. Ng et al. [31] showed that the maximum likelihood estimator (MLE) of α is biased. The constant-bias-correcting (CBC) parametric bootstrap recommended by Lemonte et al. [32] can be applied to reduce the bias of α, and so we used it to construct the BCI for τ. Let x = (x 1 , x 2 , ..., x n ) be an original random sample of size n drawn from BS(α, β), which has distribution function F = F α,β (x).α andβ, the MLEs of α and β, respectively, can be obtained by maximizing the log-likelihood function by applying the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton nonlinear optimization algorithm. Thus, the log-likelihood function of a BS distribution without the additive constant becomes Assume that x * = (x * 1 , x * 2 , ..., x * n ) is a sequence of bootstrap samples of size n obtained from BS(α,β). Thus, the bias of estimatorα is given by Suppose that B bootstrap samples are available, then theirα series (denoted byα * 1 ,α * 2 , ...,α * B ) can be calculated, whereα * j is a sequence of the bootstrap MLEs of α, for j = 1, 2, ..., B. The BFGS quasi-Newton nonlinear optimization algorithm is also used to compute the bootstrap MLEs of α. Following this, the bootstrap expectation E(α) can be approximated by using the meanα * (·) = 1/B ∑ B j=1α * j . Therefore, the bootstrap bias estimate based on B replications ofα is calculated asB (α, α) =α * (·) −α. (14) By using the bootstrap bias estimate, the correct estimate forα * (denoted asα) is given bỹ Note that we subtractB(α, α) twice because the bootstrap samples are generated by using the bias estimate (α,β) and theα * s are biased estimates ofα [33]. By applying Equation (5), the bootstrap estimator of τ can be expressed aŝ Therefore, the 100(1 − γ)% confidence interval for τ based on BCI becomes whereτ(ν) denotes the 100ν% percentile ofτ.

Bayesian Credible Interval
To guarantee proper posterior distributions, inverse gamma (IG) priors with known hyperparameters denoted as IG(β|a 1 , b 1 ) and IG(α 2 |a 2 , b 2 ) are utilized as priors for β and α 2 , respectively [34]. Since showing that β and α 2 are independent of each other is straightforward, then their respective priors can be written as and Let x = (x 1 , x 2 , ..., x n ) be a sample from BS(α, β), then the likelihood function is given by From Equations (18) to (20), the joint posterior density of β and α 2 can be obtained as Hence, it is clear that and The availability of the marginal posterior distribution of β in Equation (22) and the conditional posterior distribution of α 2 given β in Equation (23) allow the application of a Markov chain-Monte Carlo method to generate the posterior samples. Therefore, the Rpackage of LearnBayes was applied to obtain the samples by applying Equation (23) while the generalized ratio-of-uniforms method of Wakefield et al. [35] was used to generate the samples by applying Equation (22) [34]. The concept of the generalized ratio-of-uniforms method is covered here.
, where r ≥ 0 is a constant and π(·|x) is given in Equation (22). If (u, v) are random variables uniformly distributed in C(r), then the ratio β = v/u r has density π(β|x)/ π(β|x)dβ. The accept-reject method from the minimal bounding Note that a(r) and b + (r) are finite, while b − (r) = 0 [34]. Hence, the procedure for using the confidence intervals constructed with BayCI for the CV is as follows. First, set the values for a 1 , b 1 , a 2 , b 2 and r, then compute a(r) and b + (r) by using Equations (24) and (26), respectively. At the ith step, generate the values for u and v from U(0, a(r)) and U(0, b + (r)), respectively, where U(g, h) is a uniform distribution with parameters g and h. Subsequently, compute ρ = v/u r . If u ≤ [π(ρ|x)] 1/(r+1) , then set β (i) = ρ; otherwise, regenerate u∼U(0, a(r)) and v∼U(0, b + (r)). Next, generate a new value for α 2 (i) from the IG distribution given in Equation (23) depending on β (i) . Hence, the value of α (i) can be obtained after a simple algebraic transformation. Following this, the Bayesian estimator of τ becomeŝ The preceding process is repeated N times. Therefore, the 100(1 − γ)% confidence interval for τ based on BayCI becomes whereτ * (ν) denotes the 100ν% percentile ofτ * . The HPD interval was also considered for constructing the confidence interval for τ. The HPD interval, explained by Box and Tiao [36], is the shortest interval containing (1 − γ)100% of the posterior probability. Each point inside the interval has a higher probability density than any point outside of it. Therefore, the package HDInterval version 0.2.2 from R version 3.5.1 was utilized to compute the HPD interval for τ in the simulations and computations.

Confidence Intervals for the Difference between the CVs of BS Distributions
In this section, we present the confidence interval for the difference between the CVs of BS distributions based on the concepts of GCI, BCI, BayCI, and the HPD interval.
Suppose that X = (X 1 , X 2 , ..., X n ) and Y = (Y 1 , Y 2 , ..., Y m ) are two random samples of size n and m from a BS distribution (denoted as BS(α, β) and BS(α y , β y ), respectively). The CV of X can be calculated by applying Equation (5). By using random variable Y, the expected value and variance are given by E(Y) = β y (1 + 1 2 α 2 y ) and Var(Y) = (α y β y ) 2 (1 + 5 4 α 2 y ), respectively. Thus, the CV of Y (denoted as τ y ) becomes Since X and Y are independent, the difference between the CVs (denoted as η) can be written as

Generalized Confidence Interval
The GPQ of η is considered for constructing the confidence interval based on GCI. Let R α and R α y be the GPQs of α and α y , respectively. Following this, R α is calculated by applying Equation (9) and R α y can be obtained from where , T y ∼t(m − 1) and R β y (y; T y ) is the GPQ of β y , which can subsequently be defined as Meanwhile, β y,1 and β y,2 are the two solutions for where By applying Equation (30), the GPQ for η can be written as Therefore, the 100(1 − γ)% confidence interval for η based on GCI is where R η (ν) denotes the 100ν% percentile of R η .

Bootstrap Confidence Interval
Suppose that B bootstrap samples are available; then, the bootstrap estimator of η (denoted asη) can be defined aŝ whereα is calculated by applying Equation (15) andα y can be obtained by executing the following steps. Let y = (y 1 , y 2 , ..., y m ) be an original random sample of size m drawn from BS(α y , β y ). The MLEs of α y and β y (denoted asα y andβ y , respectively) are also obtained by maximizing the log-likelihood function by using the BFGS quasi-Newton nonlinear optimization algorithm. Suppose that y * = (y * 1 , y * 2 , ..., y * m ) are bootstrap samples of size m drawn from BS(α y ,β y ), then the bootstrap MLEs of α y given byα * y,1 ,α * y,2 , ...,α * y,B can be computed. Next, the bootstrap bias estimate based on B replications ofα y is obtained aŝ whereα * (·) y y,j . Therefore, the corrected estimate forα * y can be written as By substitutingα y in Equation (36), the 100(1 − γ)% confidence interval for η based on BCI becomes whereη(ν) denotes the 100ν% percentile ofη.

Bayesian Credible Interval
Based on the Bayesian method, the IG priors with known hyperparameters denoted as IG(β y |c 1 , d 1 ) and IG(α 2 y |c 2 , d 2 ) are considered as the priors for β y and α 2 y , respectively. Consequently, the marginal posterior distribution of β y becomes and the posterior conditional distribution of α 2 y given β y becomes The process for applying BayCI for the difference between the CVs can be summarized in the following steps.

1.
Set the values for a 1 , b 1 , a 2 , b 2 , c 1 , d 1 , c 2 , d 2 , r and r y , where r y > 0 is a constant.

Repeat
Step 3 N times.

Simulation Studies
A Monte Carlo simulation study was conducted by using R statistical software to evaluate the performances of the four methods used for constructing confidence intervals for the CV and the difference between the CVs of BS distributions under various combinations of parameters. We evaluated the performances of GCI, BCI, BayCI, and the HPD interval by measuring their coverage probabilities and average lengths based on 5000 independently generated replications, with 5000 pivotal quantities for GCI, B = 500 for BCI, and N = 5000 for BayCI and the HPD interval. Note that we set hyperparameters a 1 = b 1 = a 2 = b 2 = c 1 = d 1 = c 2 = d 2 = 10 −4 and r = r y = 2 for BayCI and the HPD interval [34]. For the nominal confidence level of 0.95, the best-performing method has a coverage probability close to or greater than 0.95, and the shortest average length. Since β is the scale parameter and the median of the BS distribution, β = β y = 1 without loss of generality was applied in this simulation study.
For the CV of a BS distribution, we used sample sizes n = 10, 20, 30, 50, or 100 and α = 0.1, 0.25, 0.5, 0.75, 1, or 2. The simulation results reported in Table 1 show that the coverage probabilities of GCI, BayCI, and the HPD interval were greater than or close to 0.95, even for a small sample size and/or a high value of α. Conversely, although BCI had the shortest average lengths, its coverage probabilities were the lowest and under 0.95 but improved when n was increased. When considering the average lengths of the other three methods, the HPD interval outperformed the others in all cases. In addition, the average lengths of the four methods tended to decrease and were similar when the sample size was increased. For the difference between the CVs of BS distributions, we used sample sizes (n, m) = (10, 10), (20,20), (30,30), (50, 50), (100, 100), (10,20), (30,20), (30, Tables 2 and 3, respectively. The coverage probabilities of GCI, BayCI, and the HPD interval were greater than or close to 0.95 for all cases irrespective of whether n = m ( Table 2) or n = m (Table 3). For both equal and unequal sample sizes, the coverage probabilities for BCI were the lowest, whereas its average lengths were similar to the others or the shortest in all cases. Meanwhile, the average lengths of the HPD interval were shorter than GCI and BayCI under all circumstances. Moreover, the performances of the four methods in terms of the average length improved and were close to each other when sample sizes (n, m) were increased.

An Empirical Application
The BS distribution has been successfully used to analyze air pollution concentration data, as its properties are similar to those of the lognormal distribution. For example, Leiva et al. [37] applied the BS distribution to examine sulfur dioxide concentration data and, later on, PM10 (particulate matter 10 µm) concentration data in Santiago, Chile [38]. In the second study, the authors proposed a criterion based on a control chart attribute for assessing the environmental risk.
In Chiang Mai, Thailand, air pollution from agricultural burning and forest fires from February to May has become a serious problem. We used datasets of PM 2.5 concentration data from Chiang Mai collected in March and April 2019 to illustrate the efficacies of the confidence intervals for the CV and the difference between the CVs of BS distributions derived using GCI, BCI, BayCI, and the HPD interval. The average daily PM 2.5 concentrations were measured at 9.00 AM in the area of Chang Phueak, Chiang Mai, Thailand. The datasets were obtained from the Pollution Control Department [39]. Since the data comprise positive values, they can be fitted to lognormal, exponential, gamma, Weibull, or BS distributions. The Akaike information criterion (AIC) and Bayesian information criterion (BIC) were applied to check the fitting of the data to these distributions. The results in Tables 4 and 5 show that the AIC and BIC of the BS distribution were the smallest, thereby ensuring its suitability for application to these datasets. To construct the BayCI and the HPD interval for the CV and difference between CVs of BS distributions using real data, we applied a 1 = b 1 = a 2 = b 2 = 10 −4 and r = 2 for both sampling areas. The sample mean, sample variance, and CV of the data are 96.0645, 2896.9290, and 0.5603 for the March dataset and 78.3000, 887.6655, and 0.3805 for the April dataset, respectively. Subsequently, the difference between the CVs was 0.1798. The 95% confidence interval based on GCI, BCI, BayCI, and the HPD interval for the CV and the difference between the CVs for the BS distributions are reported in Table 6. From the results, it can be seen that BCI had the shortest average length, while that of the HPD interval was shorter than GCI and BayCI, which is in agreement with the results from the simulation study. Once again, BCI attained the lowest coverage probability, so it is not recommended for constructing confidence intervals for the CV and the difference between the CVs of the BS distributions of the two real datasets. Meanwhile, the coverage probabilities of the HPD interval for the CV and the difference between the CVs of the BS distributions were greater than or close to 0.95. Thus, under these circumstances, the HPD interval provided the best-performing confidence intervals for the CV and the difference between the CVs of the BS distributions of these two datasets.

Conclusions
We proposed confidence intervals for the CV and the difference between the CVs of BS distributions using the concepts of GCI, BCI, BayCI, and the HPD interval. The coverage probabilities and average lengths of the four methods were evaluated through Monte Carlo simulations. The results show that, for all of the scenarios tested, the HPD interval had a reasonable coverage probability, while its average lengths were shorter than those of GCI and BayCI, and so it can be recommended for constructing the confidence intervals for the CV and the difference between the CVs of BS distributions. Although BCI had the shortest average lengths in all cases, its coverage probabilities were the lowest and under 0.95. Therefore, it cannot be recommended due to this shortcoming. Moreover, the average lengths of GCI, BayCI, and the HPD interval were similar when the sample sizes increased. Hence, GCI and BayCI can be considered as alternative methods to construct the confidence intervals for these two scenarios. In future research, we will investigate the confidence interval for the mean and the CV, as well as the difference between the means and the CVs of BS distributions with censored data.