Research Article | Open Access
Mohan D. Pant, Todd C. Headrick, "A Method for Simulating Burr Type III and Type XII Distributions through -Moments and -Correlations", International Scholarly Research Notices, vol. 2013, Article ID 191604, 14 pages, 2013. https://doi.org/10.1155/2013/191604
A Method for Simulating Burr Type III and Type XII Distributions through -Moments and -Correlations
This paper derives the Burr Type III and Type XII family of distributions in the contexts of univariate -moments and the -correlations. Included is the development of a procedure for specifying nonnormal distributions with controlled degrees of -skew, -kurtosis, and -correlations. The procedure can be applied in a variety of settings such as statistical modeling (e.g., forestry, fracture roughness, life testing, operational risk, etc.) and Monte Carlo or simulation studies. Numerical examples are provided to demonstrate that -moment-based Burr distributions are superior to their conventional moment-based analogs in terms of estimation and distribution fitting. Evaluation of the proposed procedure also demonstrates that the estimates of -skew, -kurtosis, and -correlation are substantially superior to their conventional product moment-based counterparts of skew, kurtosis, and Pearson correlations in terms of relative bias and relative efficiency—most notably when heavy-tailed distributions are of concern.
Burr  introduced twelve cumulative distribution functions (cdfs) with the primary purpose of fitting distributions to real-world data. Two popular Burr cdfs are referred to as the Burr Type III and Type XII distributions. The specific forms of these two cdfs are given as [1, Equations (11), (20)] where , , and are real-valued parameters that determine the mean, variance, skew, and kurtosis of a distribution. The parameter is negative for the Type III and positive for the Type XII distributions, whereas the parameter is positive for both the distributions .
The Burr Type III and Type XII distributions attract special attention because they include several families of nonnormal distributions (e.g., the Gamma distribution) with varying degrees of skew and kurtosis [2–5]. Further, these distributions have been used in a variety of applied mathematics contexts. Some examples include modeling events associated with forestry [6, 7], fracture roughness [8, 9], life testing [10, 11], operational risk , option market price distributions , meteorology , modeling crop prices , software reliability growth , reliability analysis , and in the context of Monte Carlo or simulation studies (e.g., ).
The quantile functions associated with (1) are expressed as [2, Equations (5), (6)] where with pdf and cdf as 1 and , respectively. The shape of a Burr distribution associated with (2) or (3) is contingent on the values of the shape parameters ( and ), which can be determined by simultaneously solving equations (16) and (17) from [2, p. 2211] for given values of skew and kurtosis.
In order for (2) or (3) to produce a valid Burr pdf, the quantile function is required to be a strictly increasing monotone function . This requirement implies that an inverse function exists. As such, the cdf associated with (2) or (3) can be expressed as , and subsequently differentiating this cdf with respect to will yield the parametric form of the pdf for as . We would also note that the simple closed form expressions for the pdfs associated with (1) can be given as 
Some of the problems associated with conventional moment-based estimates are that they can be (a) substantially biased, (b) highly dispersed, or (c) influenced by outliers [18, 19] and thus may not be good representatives of the true parameters. To demonstrate, Figure 1 gives the graphs of the pdf and associated with Burr Type XII distribution with skew and kurtosis . These values of skew and kurtosis have been used in a number of studies [18, 20–22]. Table 1 gives the parameters and sample estimates of skew and kurtosis for the distribution in Figure 1. Inspection of Table 1 indicates that the bootstrap estimates ( and ) of skew and kurtosis ( and ) are substantially attenuated below their corresponding parameter values with greater bias and variance as the order of the estimate increases. Specifically, for sample size of , the values of the estimates are only 63.53%, and 22.81% of their corresponding parameters, respectively. The estimates ( and ) of skew and kurtosis ( and ) in Table 1 were calculated based on Fisher’s -statistics formulae (see, e.g., [23, pages 299-300]), currently used by most commercial software packages such as SAS, SPSS, and Minitab, and so forth, for computing the values of skew and kurtosis (where for the standard normal distribution).
Another unfavorable quality of conventional moment-based estimators of skew and kurtosis is that their values are algebraically bounded by the sample size such that and . This constraint implies that if a researcher wants to simulate non-normal data with kurtosis as in Table 1, and drawing a sample of size from this population, the largest possible value of the computed estimate of kurtosis is only 15, which is only 71.43% of the parameter.
The method of -moments introduced by Hosking  is an attractive alternative to conventional moments and can be used for describing theoretical probability distributions, fitting distributions to real-world data, estimating parameters, and testing of hypotheses [18, 19, 24, 25]. In these contexts, we note that the -moment based estimators of -skew, -kurtosis, and -correlation have been introduced to address the limitations associated with conventional moment-based estimators [18, 19, 24–31]. Some qualities of -moments that make them superior to conventional moments are that they (a) exist for any distribution with finite mean, and their estimates are (b) nearly unbiased for any sample size and less affected from sampling variability, (c) more robust in the presence of outliers in the sample data, and (d) are not algebraically bounded by the sample size [18, 19, 24, 25]. For example, the estimates ( and ) of -skew and -kurtosis ( and ) in Table 2 are relatively closer to their respective parameter values with much smaller variance compared to their conventional moment-based counterparts as in Table 1. Inspection of Table 2 shows that for the sample size of , the values of the estimates are on average 97.71% and 97.28% of their corresponding parameters.
In view of the above, the main purpose of this study is to characterize the Burr Type III and Type XII distributions through the method of -moments in order to obviate the problems associated with conventional moment-based estimators. Another purpose of this study is to develop an -correlation-based methodology to simulate correlated Burr Type III and Type XII distributions. Specifically, in Section 2, a brief introduction to univariate -moments is provided. The systems of equations associated with the Type III and Type XII distributions are subsequently derived for determining the shape parameters ( and ) for user specified values of -skew and -kurtosis . In Section 3, a comparison between conventional and -moment-based Burr Type III and Type XII distributions is presented in the contexts of estimation and distribution fitting. Numerical examples based on Monte Carlo simulation are also provided to confirm the methodology and demonstrate the advantages that -moments have over conventional moments. In Section 4, an introduction to the coefficient of -correlation is provided, and the methodology for solving for intermediate correlations for specified -correlation structure is subsequently presented. In Section 5, the steps for implementing the proposed -moment procedure are described for simulating non-normal Burr Type III and Type XII distributions with controlled skew (-skew), kurtosis (-kurtosis), and Pearson correlations (-correlations). Numerical examples and the results of simulation are also provided to confirm the derivations and compare the new procedure with the conventional moment-based procedure. In Section 6, the results of the simulation are discussed.
2.1. Theoretical and Empirical Definitions of -Moments
-moments can be expressed as certain linear combinations of probability weighted moments (PWMs). Let be identically and independently distributed random variables each with pdf , cdf , and the quantile function , then the PWMs are defined as [18, Equation (6)] where . The first four -moments associated with can be expressed in simplified forms as [25, pages 20–22] where the coefficients associated with in (7) are obtained from shifted orthogonal Legendre polynomials and are computed as in [25, pages 20–22] or in [18, pages 4–5].
The notations and denote the location and scale parameters. Specifically, in the literature of -moments, is referred to as the -location parameter which is equal to the arithmetic mean, and (>0) is referred to as the -scale parameter and is one-half of Gini’s coefficient of mean difference [23, pages 47-48]. Dimensionless -moment ratios are defined as the ratios of higher-order -moments (i.e., and ) to . Thus, and are, respectively, the indices of -skew and -kurtosis. In general, the indices of -skew and -kurtosis are bounded in the interval , and as in conventional moment theory, a symmetric distribution has -skew equal to zero . The boundary region for -skew and -kurtosis for a continuous distribution is given by the inequality (see, )
Empirical -moments for a sample (of size ) of real-world data are expressed as linear combinations of the unbiased estimators of the PWMs based on sample order statistics . Specifically, the unbiased estimators of the PWMs are given as [19, pages 113-114] where and is the sample mean. The first four sample -moments are obtained by substituting from (9) instead of in (7). The sample -moment ratios (i.e., -skew and -kurtosis) are denoted by and , where and .
2.2. -Moments for Burr Type III Distributions
Substituting from (2), , and in (6), the th PWM for the Burr Type III distributions is given by Equation (10) can also be expressed as Let . Then, . Substituting in (11), the th PWM can be expressed as The integral in (12) is a beta function, , where and such that .
2.3. -Moments for the Burr Type XII Distributions
Substituting from (3), , and in (6), the th PWM for the Burr Type XII distributions is given by After some manipulations, (16) can be expressed as Let . Then, . Substituting in (17), the th PWM can be rewritten as
Given specified values of -skew and -kurtosis the systems of equations (15) and (21) can be simultaneously solved for real values of and . The solved values of and can be substituted in (2) and (3), respectively, for generating the Burr Type III and Type XII distributions. Further, the solved values of and can be substituted in (13)(14) and (19)(20), respectively, for computing the values of -mean and -scale associated with the Type III and Type XII distributions. In the next section, two examples are provided to demonstrate the aforementioned methodology and the advantages that -moments have over conventional moments in the contexts of estimation and distribution fitting.
3. Comparison of -Moments with Conventional Moments
An example to demonstrate the advantages of -moment-based estimation over conventional moment-based estimation is provided in Figure 2 and Tables 4–7. Given in Figure 2 are the pdfs of the distribution (3, 10), Chi-square (), extreme value (0, 1), and logistic (0, 1) distributions superimposed, respectively, by the Burr Type XII, Type III, Type XII, and Type III pdfs (dashed curves) in both (a) conventional moment- and (b) -moment-based systems. The conventional moment-based parameters of skew and kurtosis associated with these four distributions, given in Table 4, were computed by using equations (11)−(13) from [2, page 2211]. The values of shape parameters ( and ) given in Table 4 were determined by solving equations (16) and (17) from [2, page 2211]. The values of and were used in (4) and (5) to superimpose the conventional moment-based Burr Type XII, Type III, Type XII, and Type III distributions as shown in Figure 2(a).
The -moment-based parameters of -skew and -kurtosis associated with the four distributions in Figure 2, given in Table 5, were obtained in three steps as follows: (a) compute the values of PWMs () using (6), (b) substitute these PWMs into (7) to obtain the values of the first four -moments, and (c) compute the values of and using and . The values of shape parameters ( and ) given in Table 5 were determined by solving the systems of equations (15) and (21), respectively. These values of and were used in (4) and (5) to superimpose the -moment-based Burr Type XII, Type III, Type XII, and Type III distributions as shown in Figure 2(b).
To superimpose the Type III or Type XII distribution, the quantile function in (2) or (3) was transformed into (a) , and (b) , respectively, where , and , are the values of mean, standard deviation, whereas , and , are the values of -mean, -scale obtained from the original distribution and the respective Burr Type III or Type XII approximation, respectively.
The advantages of -moment-based estimators over those based on conventional moments can also be demonstrated in the context of Burr Type III and Type XII distributions by considering the Monte Carlo simulation results associated with the indices for the percentage of relative bias (RB%) and standard error (St. error) reported in Tables 6 and 7.
Specifically, a Fortran  algorithm was written to simulate 25,000 independent samples of sizes and , and the conventional moment-based estimates ( and ) of skew and kurtosis ( and ) and the -moment based estimates ( and ) of -skew and -kurtosis ( and ) were computed for each of the samples based on the parameters and the values of and listed in Tables 4 and 5. The estimates ( and ) of and were computed based on Fisher’s -statistics formulae [23, pages 47-48], whereas the estimates ( and ) of and were computed using (7) and (9). Bias-corrected accelerated bootstrapped average estimates (Estimate), associated 95% confidence intervals (95% Bootstrap C.I.), and standard errors (St. error) were obtained for each type of estimates using 10,000 resamples via the commercial software package Spotfire S+ . Further, if a parameter was outside its associated 95% bootstrap C.I., then the percentage of relative bias (RB%) was computed for the estimate as
The results in Tables 6 and 7 illustrate that the -moment-based estimators are superior to their conventional moment-based counterparts in terms of both smaller relative bias and error. These advantages are most pronounced in the context of smaller sample sizes and higher order moments. For example for the Distribution 1, given a sample of size , the conventional moment-based estimates ( and ) generated in the simulation were, on average, 41.46% and 5.62% of their corresponding parameters ( and ). On the other hand, for the same Distribution 1, the -moment-based estimates ( and ) generated in the simulation study were, on average, 93.88% and 92.03% of their corresponding parameters ( and ). Thus, the relative biases of estimators based on -moments are essentially negligible compared to those associated with the estimators based on conventional moments. Also, it can be verified that the standard errors associated with the estimates and are relatively much smaller and more stable than the standard errors associated with the estimates and .
Inspection of the graphs in Figures 2(a) and 2(b) and the Monte Carlo simulation results in Tables 6 and 7 indicate that the -moment-based Type III and Type XII pdfs provide a more accurate approximation of the four distributions than those based on conventional moment theory.
3.2. Distribution Fitting
Figure 3 shows the conventional moment- and the -moment-based Burr Type XII pdfs superimposed on the histogram of ankle circumference data obtained from 252 adult males (http://lib.stat.cmu.edu/datasets/bodyfat) as cited in Headrick [35, page 48].
The conventional moment-based estimates ( and ) of skew and kurtosis ( and ) and the -moment-based estimates ( and ) of -skew and -kurtosis ( and ) were computed for the sample of size participants. The estimates of and were computed based on Fisher’s -statistics formulae [23, pages 47-48], whereas the estimates of and were computed using (7) and (9), respectively. These sample estimates were then used to solve for the values of shape parameters () using (a) equations (16)−(17) from [2, page 2211] and (b) (21). The solved values of were subsequently used in (5) to superimpose the parametric plots of the Burr Type XII pdfs shown in Figure 3.
Inspection of the two panels in Figure 3 illustrates that the -moment-based Type XII pdf provides a better fit to the sample data. The Chi-square goodness of fit statistics along with their corresponding values given in Table 3 provide evidence that the conventional moment-based Type XII distribution does not provide a good fit to the actual data, whereas the -moment-based Type XII distribution fits very well. The degrees of freedom for the Chi-square goodness of fit tests were computed as (class intervals) − 4 (estimates of the parameters) − 1 (sample size).
4. -Correlations for the Burr Type III and Type XII Distributions
Let and be two random variables with cdfs and , respectively. The second -moments of and can be defined as  The second -comoments of toward and toward are given as The -correlations of toward and toward are subsequently defined as The -correlation given in (25) (or, (26)) is bounded in the interval . A value of implies that and have a strictly and monotonically increasing (or decreasing) relationship. See Serfling and Xiao  for further details on the topics related to the -correlation.
The extension of the Burr Type III and Type XII distributions to multivariate level can be obtained by specifying quantile functions as given in (2) and (3) with a specified -correlation structure. Specifically, let denote standard normal variables with cdfs and the joint pdf associated with and given by the following expressions: where in (29) is the intermediate correlation between and . Using the cdfs in (27) and (28) as zero-one uniform deviates, that is, , , the quantile function defined in either (2) or (3) can be expressed as a function of , or (e.g., or ). Thus, the -correlation of toward can be determined using (25) with the denominator standardized to for the standard normal distribution as The variable in (30) is the standardized quantile function of (2) or (3) such that it has an -mean (or arithmetic mean) of zero and -scale equal to that of the standard normal distribution. That is, the quantile function is standardized by a linear transformation as where is the mean from (13) or (19) and is a constant that scales in (14) or (20) and in the denominator of (25) to . In particular, for the Type III and Type XII distributions can be expressed as The next step is to use (30) to solve for the values of the ICs such that the specified Type III and Type XII distributions have their specified -correlation structure.
Analogously, the -correlation of toward is given as Note that in general, the -correlation of toward in (30) is not equal to the -correlation of toward in (34). These -correlations are equal only when the values of shape parameters and associated with and are equal (i.e., when the two distributions are the same). Provided in Algorithm 1 is a source code written in Mathematica [36, 37], which shows an example for computing for the -correlation procedure. The steps for simulating correlated Burr Type III and Type XII distributions with specified values of -skew , -kurtosis , and with specified -correlation structure are given in Section 5.
5. The Procedure for Monte Carlo Simulation with an Example
The procedure for simulating Burr Type III and Type XII distributions with specified -moments and -correlations can be summarized in the following six steps.(1)Specify the -moments for transformations of the form in (2) and (3), that is, , and obtain the solutions for the shape parameters and by simultaneously solving the systems of equations (15) and (21) for the specified values of -skew () and -kurtosis () for each distribution. Specify a matrix of -correlations () for toward , where .(2)Compute the values of intermediate (Pearson) correlations (ICs), , by substituting the value of specified -correlation () and the solved values of and from Step (1) into the left- and the right-hand sides of (30), respectively, and then numerically integrating (30) to solve for . See Algorithm 1 for an example. Repeat this step separately for all pairwise combinations of ICs.(3)Assemble the ICs computed in Step (2) into a matrix and then decompose this matrix using Cholesky factorization. Note that this step requires the matrix to be positive definite. (4)Use elements of the matrix resulting from Cholesky factorization of Step (3) to generate standard normal variables () correlated at the IC levels as follows: where are independent standard normal random variables, and where is the element in the th row and th column of the matrix resulting from Cholesky factorization of Step .(5)Substitute from Step into the following Taylor series-based expansion for computing the cdf, , of standard normal distribution  where is the pdf of standard normal distribution and the absolute error associated with (36) is less than . (6)Substitute the uniform variables, , generated in Step into the equations of the form