What Determines EP Curve Shape?

Propose use kurtosis divided by skewness squared as shape factor, and use the global or conditional minimum/maximum of this shape factor for selecting and differentiating distribution families. Semi-empirical formulas for that lower/upper bound are calculated for various distribution families, with the aid of Computer Algebra System, for fitting hard to match distributions. Previous studies show high CV distribution is hard to fit and simulate, this study extends that conclusion to cases with low CV but still hard to match EP curves, characterized by having shape factors close to 1. The maximal likelihood approach of distribution fit can tell us which distribution family is better suited for an empirical distribution, but the shape factor range information can tell us why a distribution cannot fit well, or is not suitable at all. So the shape factor, in a sense, determines the EP curve shape.


Introduction
In reinsurance industry, losses for a contract are simulated and represented by the losses cumulative distribution function (CDF), survival, or quantile functions. The plots of these functions are called the EP curves with the following terminology [1]: for a given annual or aggregated loss, the probability of seeing annual loss exceeding that loss is the exceeding-probability (EP) or aggregate-exceedingprobability (AEP). The average of all annual losses exceeding that given loss is the AEP tail value at risk, called the AEP TVaR, or simply TVaR. The EP curve is represented by a table consisting of pairs of probability and loss. It is desirable to fit a parametric distribution to this table for a more succinct representation and more reasonable interpolations for values not in the table. Then which distribution family to use and what characteristics of the data are needed or determine the distribution are the questions to answer.
The (scaled) Beta distribution is widely used in reinsurance for fit loss or loss ratio, perhaps because the Beta distribution has only two parameters and very simple formulas for mean and standard deviation using these parameters, whose inverse function also has simple formulas, so that the two statistics of mean and standard deviation can be used to easily determine the parameters.
For about 85% of the perils, this approach works well, in the sense that the TVaR of the fitted distribution for quantile of interest, such as the 0.96, 0.99 or 0.996 quantile TVaR which is needed for pricing and risk monitoring, is close to a few percent of the original data TVaR. The remaining 15% perils, such as the North American Tornado Hail (NATH), Australia Wind Storm (AUWS), Hawaii Wind Storm (HIWS), and Mexico Earthquake (MXEQ), can have more than 10% deviations.
The maximum likelihood estimation method is a way to find alternative fitting distributions [2,3]. Instead of finding approximations of the smoothed empirical distribution, we optimize an objective function whose optimum solution gives us the candidate distribution form. Suppose the annual losses x i occurred n i times in our observation; to find a probability function that gives probabilities p i for these losses, we just maximize the objective function Q i p i n i . It is easily seen that for the optimum solution we have p i p j ¼ n i n j : the relative occurring frequency is maintained in the probability function. In the objective function, if we replace the p i by a power function of p i , the conclusion still holds, but not if we use a logarithm or exponential function.
While the maximum likelihood approach works well for many perils and identifies a few best fitted distribution families (Mathematica has more than 200 distribution families that can be used for extensive searches), it did not work for the NATH  This is not contradictory to the maximum likelihood principal, since in any implementation, only known forms of the probability density function (PDF) and as-small-as-possible numbers of parameters can be used. To overcome this limitation, we need to look into the particularity of those distributions and come up with or select more suitable function forms for the PDF or CDF. In [4] it is found that a high coefficient of variation (CV) distribution is hard to fit or simulate. But the NATH has a small CV of 1.28. The skewness and kurtosis alone also not differentiate them from other distributions.
Trial and error found the empirical rule that these hard distributions have small values of kurtosis divided by skewness squared, Table 1. This finding prompted us for the study of the property of kurtosis/skewness^2 (K/S^2), henceforth will be called the shape factor (SF).
Numerical optimization or solution will be our primary tool for this SF study. Analytical deduction, symbolic algebra, and symbolic limit from computer algebra system (CAS) Mathematica will be another major tool, as well as Mathematica's plot functions. Those plots can help reveal the patterns or tendencies of functions. The found pattern can in turn aid in taking special directional/constraint limit or substitutions in CAS to get the analytical formula for SF bound when it is possible.
The overall lower bound we find of SF is presented in Section 2, through the triple analytical, graphical, and numerical methods. Followed by in-detail studies of SF of various selected distribution families, which are either widely used in practice, such as the Beta distribution in Section 3 and the generalized Gamma distribution in Section 6, or is most simple to simulate, such as the Kumaraswamy distribution in Section 4. The most inclusive distribution, BetaPrime distribution, is in Section 5, for which we do not get an analytical formula, so the empirical formula for SF lower bound is provided. Some distributions that have wide matching capabilities, but for the NATH may have fitted distribution facing numerical difficulties, such as the Fleishman distribution, whose fit has non-monotonically increasing polynomial form and hence is hard to solve for inverse CDF, are only briefly mentioned in Section 7. The top distribution found through maximum likelihood fit, the generalized hyperbolic distribution (GH), even with the most complex PDF, has unexpectedly simple and beautiful analytical formulas for SF lower bound; the results are in the final Section 8. All our studies will focus on SF bound deductions and applications.
This property of the shape factor identified our exceptional perils as possessing very steep and asymmetric PDF whose SF are small.

Are there better definitions of shape factor?
To measure the steepness or the convexity, we can get similar inequality to Eq. (2) by Hölder inequality for absolute moment: From absolute central moment define: For nonnegative random variables such as the reinsurance contract loss distribution, use the following inequality for moment: From another application of Hölder inequality, we get yet other measures of convexity from absolute central moment: ½ 2 is when the numerical calculation error with extreme parameters arrive at negative kurtosis, then the calculated SF are meaningless (An example of BesselK function inaccuracy brings about negative kurtosis for generalized hyperbolic distribution can be found in [8]).
Even though both SF and SF1 are invariant under linear transformation of the distribution, and both measure the convexity, SF ≥ SF1 can additionally measure the asymmetry, combining these two into one quantity. Since most distributions in reinsurance are not symmetric, SF is preferred over SF1. That only SF measured both asymmetry and convexity, while the others cannot, can be seen from Figure 1, for the case of exponential distribution family with PDF e Àx n nx À1þn , x ∈ 0; ∞ ð Þ, n > 0, which is WeibullDistribution n; 1 ½ or GammaDistribution 1; 1; n; 0 ½ , where only SF has a nontrivial interior global minimum.
An intuitive reason for why using shape factor SF in favor of skewness and kurtosis alone is provided by studying the simple example power distribution family with PDF nþ1 n x 1 n , x ∈ 0; 1 ½ , n , À 1∥n > 0 (or BetaDistribution 1=n þ 1; 1 ½ ). This distribution family has the largest value of skewness and kurtosis, and at the same time the smallest shape factor SF when n turns to À1, where the PDF is the steepest, but the skewness and kurtosis take the indistinguishable value of infinity. In comparison, the shape factor SF takes the finite and distribution family specific value of 1.125. The shape factor SF thus makes meaning out of the meaningless infinities.

Alternative way of defining shape factor for symmetric distribution
For symmetric distribution, CM 3 ½ ¼ 0, our SF will be indiscriminately infinity. We can now employ SF1 in place of SF. Other measures from ACM such as SF2 and SF3 may also be candidates. From the SF3 plot Figure 2 of NormalDistribution μ; σ ½ we see that min 0 , r , 1 SF3 r ½ ¼ 0:919824. The lower the value of SF3 2 ½ , the higher the min 0 , r , 1 SF3 r ½ . We can use either SF3 2 ½ or min 0 , r , 1 SF3 r ½ as a shape factor for symmetric distribution to describe the convexity of the ACM curve. The second measure  has the merit of independence to the power order r, by additional efforts of numerical minimization. For our power distribution family, the maximum of the minimum is: max n > 0 min 0 , r , 1 SF3 r ½ ¼ 0:942085, higher than the Normal distribution family.
When all SF, SF1, and SF2 are available, however, we will prefer SF to SF1 and SF2 since its dependency on parameters show simpler patterns than the other two; this can be shown from their contour plots for Beta distribution Figures 3-5, where SF contours are almost lines.

Lower bound of SF for well-known distributions
Using numerical optimization [9,10], for most of the top-fitted distributions from the maximum log likelihood approach, we get the minimum SF values, with distribution definition in [11] whose naming and parameterization for probability distributions will be used throughout this chapter, in Table 2.
From this table, we know that most of the distributions are not able to describe NATH since NATH has SF 1.834. More involved numerical integration and optimization also eliminated the Beckmann Distribution [12], with admissible SF range  3.63-8.16, being the top four-parameter-distribution in another distribution fit case study that has SF 4.58.
The Alpha-Skew-Normal Distribution from [13] has minimum SF4.95061 when α is 2.07764, from its proposition 2.3, is thus also not eligible for NATH.
The global lower bound of SF for parametric distribution can be used to filter out those distributions whose values are larger than the losses data SF, so that we can focus on distributions that do not violate the bound. In the following sections we will study typical distribution SF bound, beginning with the Beta distribution.

Beta distribution
Regardless of the fact that multitude distribution types have been used for the frequency and severity distribution of individual contract losses, the aggregated portfolio losses for the majority of perils can be fitted by a compound Poisson distribution with Beta distribution as the severity, somehow an attest of its prevalence. Beta distribution has min SF ¼ 1:0, so we need an in-detail study of why it cannot fit NATH.
When matching a BetaDistribution α; β ½ for skewness 5.99378 and kurtosis 65.8902, we must have β < 0. When matching a Beta distribution for CV(=std/ mean, the standard deviation divided by the mean) 1.2829 and either skewness

Distribution
Min SF Location of the Min Table 2.
Lower bound of SF for some well-known distributions.
5.99378 or kurtosis 65.8902, we must have either both α < 0 and β < 0 or at least one of α or β less than 0. Since CV, skewness, and kurtosis are scale invariant, so no scaled Beta distribution can at the same time match any two of the three statistics CV, skewness, and kurtosis.

Minimum shape factor for given CV
Using Mathematica, we can solve the parameter α and β by cv and std for BetaDistribution α; β ½ : α ! cv À std À cv 2 std cv 3 , β ! cv 2 À 2cvstd À cv 3 std þ std 2 þ cv 2 std 2 cv 3 std : Since α > 0, we must have: We also know std must be between 0 and 0.5 for these solutions to exist. By computer-aided exploration through contour plot, we can find the location of the std where SF takes minimum for a given cv.
The overall observation is that when cv < 1, SF approaches infinity in the middle value of std, and decreases when deviating from it. When cv > 1, SF approaches its minimum in the middle value of std and increases when departing. Together with the fact that std has an allowable upper bound of cv/(1 + cv^2) and lower bound of 0, the minimum of SF must be attained either at the global extreme where the derivative of SF with respect to std is zero or at the two boundaries when cv > 1, and attained at the two boundaries when cv < 1.
Using Mathematica to take the derivative of the shape factor with respect to std to find the std where shape factor attained extreme values, and solving it for the intersection with std upper and lower bound, we know the minimal shape factor for Beta distribution for a given CV when CV is below 0.707107 or above 2.48239 (intersecting std upper bound) is attained at std upper bound cv 1þcv 2 with value: Þ 2 , when cv , 0:707107∥cv > 2:48239: When CV is between 0.707107 and 1.024766 (intersecting std lower bound) the minimal shape factor is attained at std lower bound 0 with value: min 0 , std , cv When CV is between 1.024766 and 2.48239, the minimum SF is attained at std that is the zero derivative points of the shape factor. The piecewise curve plot of the minimum SF for given CV is in Figure 6. The formula for the central piece, minshape, is given in Figure 7 which is too complex for manual derivation without the aid of computer algebra system.
From the curve we know when CV = 1.28, the minimal shape factor is 1.88, larger than 1.83 of NATH. In the best effort to match the input, we may elect to relax CV, for example, to 1.3, then the minimum shape factor is 1.85. With the constraint of a given CV, the minimum shape factor of the Beta Distribution may be significantly larger than its global minimum 1, so that it cannot attain to the wanted SF value.

Shape factor range for given skewness
By solving Beta distribution parameters α and β through skewness sk and kurtosis kt, and examining the contour plot of β, we can see the allowable region is bound by two parabolas, Figure 8.
For a fixed skewness, α is monotonic increasing with respect to kurtosis; on the other hand, β has a singular point in some kurtosis, below that kurtosis is positive and monotonic increasing(in the region where α is positive), Figure 9.
Solving for that singular point we get the permissible kurtosis upper bound 3 þ 3 2 sk 2 , and solve for β ¼ 0 get the permissible kurtosis lower bound 1 þ sk 2 .
Observe that the upper bound is when β turns to infinity, we can also get a simpler derivation of the upper bound by representing skewness and shape factor in α and β, letting β ! ∞, and then eliminating α to get shape factor as a function of skewness (Mathematica cannot solve equation for skewness which includes square root expression, we get around that by solving equation for the square of skewness, and then abandoning the negative solution introduced by this square).
A third way of more tedious calculation is through solving α by skewness and β, substituting the real solution into shape factor, and then take the limit for β ! ∞.
All three methods get the same upper bound of SF ¼ 3 2 þ 3 sk 2 :  So for Beta distribution, the allowable region of skewness and kurtosis is bound below by kurtosis = skewness^2 + 1 where β ! 0, and above by kurtosis = 3 + 1.5*skewness^2 where β ! ∞: For the given skewness of 5.99378 of NATH, the maximum allowable kurtosis is 56.88813, less than the wanted 65.8902. So NATH cannot be fitted by any affine transformation of Beta distribution, certifying NATH as a trying case for distribution fitting. We will use it to test many of the well-known distributions in later sections. We also see surprisingly that unlike many of the other distribution families whose shape factors are too high, the Beta distributions have the shape factor range too low, or too close to 1. This suggests us to search for distributions with shape factors ranges in between.

Kumaraswamy distribution
Using the same approach as in the Beta distribution, we first study the skewness and kurtosis tendency of KumaraswamyDistribution α; β ½ [14], since the latter tested  to be a better choice in our experiment and is also the easiest for simulation, Figures 10-12; and then study the SF bound for given skewness. From these plots, we see an overall rough tendency of the skewness, kurtosis and shape factor. For a given α, the shape factor converges to a finite limit when β ! ∞. For a given skewness or a given kurtosis, there exists a maximum allowable α that is arrived when β ! ∞. In the parameters space of (α,β), for a given α, the kurtosis is increasing with respect to β in the top left portion where the skewness is positive, and decreasing in the right bottom portion where the skewness is negative. And in the parameters space of (α,β), for a given α, the shape factor is decreasing with respect to β in the top left portion where the skewness is positive, and increasing in the right bottom portion where the skewness is negative. But we will see later that the tendencies are more delicate than the monotonicity shown through visual observation.
Combining the tendency of shape factor and the contour plot for given skewness, kurtosis, and shape factor as in Figure 13, we may guess that for a given positive skewness, when α turn to its upper limit and β turn to infinity, the shape  factor will converge to its minimum. We use Mathematica to calculate the asymptotic expansion of the Gamma function and the quotient of Gamma function at infinity for orders up to 4 or 2, take the symbolic limit for β ! ∞, to get these boundary values, Figures 14 and 15.
We thus have a simple formula for boundary value of Kumaraswamy distribution shape factor:  Derivation of Kumaraswamy distribution kurtosis upper bound for given α and shape factor boundary value for given α when β ! ∞.

13
What Determines EP Curve Shape? DOI: http://dx.doi.org/10.5772/intechopen.82832 Its plot Figure 16 has two branches, the dividing point is α ! 3:602349425719043 where the skewness is zero, and below it is mainly the positive skewness region while above it is the negative skewness region.
The minimum value at the left branch of Figure 16 is 1.91227 and arrived at α = 0.641149. When α > 1000 the numerical value for that boundary can be negative and is thus unreliable. The value 1.91227 is not the global minimum of the shape factor: for α = 0.641149 the shape factor plot Figure 17 with respect to β decreases first, at the point 10.6095 arriving at the minimum value of 1.80935, and increasing after the point 10.6095.
In principle, the extreme value of the shape factor for a given skewness will arrive either at the upper boundary where β ! ∞ or at the lower boundary where  α ! 0, or at some middle point where the contour plot of the skewness and the contour plot of the kurtosis will be tangent to each other. The Mathematica contour plot does not work for a very small α, but by numerical minimization we know the global minimum of the Kumaraswamy distribution shape factor is 1.03709 when α ¼ 1:80143 * 10 À9 , β ¼ 0:247044. The conditional minimum of the shape factor when skewness = 5.99378 is about 1.04753 when α ¼ 10 À10:5 , β ¼ 0:149286 through list calculation; this is higher than 1 + 1/S^2 = 1 + 1/ 5.99378199789956^2 = 1.02784, the lower boundary of Beta distribution.
The Mathematica contour plot works for large α, and we see the shape factor is increasing along the contour of skewness, which attains its maximum when β ! ∞. For example, for NATH skewness 5.99378199789956, the maximum shape factor is 1.97131, arriving at α = 0.5239510562868946. The maximum shape factor of Kumaraswamy distribution for given skewness is in Figure 18, which is algebraically represented by the parametric curve of Eq. (6) and Eq. (8).
So the permissible shape factor range of the Kumaraswamy distribution still spans the lower end of the whole allowable range of (1,∞), but higher than that of the Beta distribution. Affine transformed Kumaraswamy distribution can fit all the first four moments of NATH, with the fitted distribution TVaR close to NATH TVaR in the error range of 5-6%, while the best effort affine transformed Beta distribution is in the error range of 9-10%.
To further improve the fit, we need additional freedom in parameters, such as the GB1 distribution [15], since KumaraswamyDistribution α; β ½ ≈GeneralizedBetaDistributionI 1; β; α; 1 ½ , and the maximum shape factor plot in Figure 18 is lower than that of LogNormalDistribution, the upper bound of GB1. The following section will study a sibling distribution to GB1, fitted as good as GB1, but is more widely known.

BetaPrime distribution
Beta distribution and Kumaraswamy distribution are a few exceptions which have analytical formulas for the shape factor bounds; for other distributions to be studied, numerical optimization and empirical plot or formula will be the only feasible approach.
The minimum shape factor of Beta Distribution is 1, but that of the transformed is 1.5: Empirically, the larger the third parameter α, the smaller the minimum shape factor. The smallest shape factor we get of the BetaPrimeDistribution is 1.125, when α = 446.49537:
With the transformation of p-> 10^w, α-> 10^-z, q-> 4*10^z + y, we can study the GB2 shape factor change tendency with respect to α, Figure 19, and shape factor change tendency with respect to p, Figure 20.
The GB2 shape factor is mainly determined by α and p, only slightly changing with respect to q when q is smaller than 5. The change with respect to α and p is similar, having two peaks, or three peaks if we regard the two sides of the infinity as two branches since that border is not easily crossable for searching or optimization algorithms.
GB2 shape factor's dependency with p and α, or w and z through transformation p = 10^w, α = 10^-z, is mostly unaffected by q except for right-most values of z. They are μ-shaped (Figure 21), this is different from Hyperbolic Distribution (in Section 8), whose shape factor dependency with λ is V-shaped. We guess V-shaped curves have unique global minimums, but μ-shaped curves will show bifurcation behavior: the converged solution in optimization will be very different when the initial point or interval is slightly different.
The knowledge that the shape factor curve attained extreme values in À3.3,-1.25 and 1 with respect to z, and attained extreme values in À2.65, À1.11 and 1 with respect to w, can be used to set the initial interval, the paramount factor determining the quality of the numerical optimization solution, for solving the GB2 fitting problem.

Minimum shape factor for GB2
The skewness and kurtosis matching problem for GB2 is very sensitive to the initial parameter ranges given. A study of the minimum shape factor of GB2 with  GB2 distribution skewness kurtosis and shape factor vs. α or z, vs. p or w plots for fixed y = q-4/α. respect to each parameter will give us permissible ranges for those parameters. Direct work with shape factor encounters problems from Mathematica's numerical optimization function NMinimize, minimizing the log shape factor instead can overcome this difficulty. The plot is in Figure 22.
In the range (0.0001, 5.0) of p, the numerical minimum shape factor plot of GB2 is a very smooth curve. The fitted formula of GB2 min SF for given p by Mathematica's machine learning function FindFormula is Eq. (9). min K S 2 ¼ 1:1593871374775397 þ 1:4702458297305288 * 0:5148499158800361 1 p 0:3215433282777008 (9) As a test, for NATH the log shape factor is Log 1:83408 ½ ¼0:60654412, the solution of Eq. (9) for p with NATH SF is p ¼ 0:608342; the minimum log shape factor of GB2 for this p is 0.60603997, only 0.08% smaller than input.
From the contour plot Figure 20 we know for given α, the shape factor of GB2 has two singular points with p or 10 w . The minimization for given α needs to carry out in each of the three regions cut by these two singular points. The plot is in Figure 23. With a new parameterization, p ¼ λ α , q ¼ 4þν α , the minimization of shape Figure 22.
The numerical minimum GB2 shape factor for given p in horizontal axis.

Figure 23.
The numerical minimum GB2 shape factor for given α or given pα in horizontal axis.
factor for GB2, for given λ ¼ pα, is easier to perform. The plot is included in Figure 23 as well. Figures 22 and 23 show that the permissible parameters for NATH are p , 0:63, α > 0:5, pα , 0:5: This is confirmed by GB2 fit practice. The best fit by GB2 for NATH is at w ¼ À0:329075005, p ¼ 0:468732, with about 5% error from input TVaR. The discontinuity of fitted GB2 TVaR with respect to parameter change is also observed, this w value is such a critical point.

Generalized gamma distribution
The generalized gamma distribution in Mathematica is the Amoroso distribution [16], with the parameter correspondence: α $ α, β $ θ, γ $ β, μ $ a: For generalized gamma distribution GammaDistribution α; β; γ; μ ½ , the shape factor depends only on α and γ. It seems the smaller the α, and the bigger the γ, the smaller the K S 2 . When α ¼ 3:318512677036329 Â 10 À12 , γ ¼ 8811:572418686921, K S 2 ¼ 1:125, close to the global minimum 1 of K/S^2. So there arises the question: the generalized gamma and GB2 can match smaller shape factors than Hyperbolic Distribution (Section 8), why they cannot fit as good as the latter for NATH with shape factor 1.83409?
One explanation is that the numerical solution for GB2 or generalized Gamma distribution is trapped in the shape factor curve right branch by the combined constraints of skewness and kurtosis, which is not the branch that can attain 1.125, unlike the generalized hyperbolic distribution whose shape factor has a global minimum in λ = 0.
The lower bound of shape factor from unbounded distributions seems, in general, to be higher than bound distributions'. Outside of the latter's upper bound and near the former's lower bound, for a SF value slightly larger than 1.5, in practice, most parametric distributions have difficult matching both the kurtosis and skewness: the comparatively best one is selected for study in the next section.

Hyperbolic distribution
Taking a sequence of numerical minimization of the shape factor, for various values of fixed λ, we get the empirical minimum shape factor curve for generalized hyperbolic distribution (GH),HyperbolicDistribution λ; α; β; δ; μ ½ , in Figure 24. We observed that when λ > À0.6, the minimum shape factor is attained when α^2-β^2-> 0 and β-> 0, that is, it is attained by a skew hyperbolic t distribution [17][18][19]. When looking at the plot of shape factor with respect to λ, we feel that it must have some simple formula. So we utilize Mathematica symbolic calculation to expand the shape factor with asymptotic expansion for BesselK λ; a ½ , or K λ a ð Þ in [20], with respect to α^2-β^2 and then take the symbolic limit, Figure 25.  Derivation of the GH shape factor limit when λ > À2.
With this knowledge, the NATH permissible λ interval is [À0.8439,0.4454], and the left end point still have 0.99 TVaR larger than the input TVaR, but now only by 4.05%, better than the 5% error of GB2.

Conclusion and discussions
We proposed using the ratio of kurtosis by squared skewness as the best candidate for shape factor that can characterize the distribution asymmetry, as well as the PDF steepness. The closer this factor to 1, the more asymmetry and the steeper the PDF. The asymptotic approximation and symbolic limit is used to calculate the boundary of this factor for various distributions: the Beta, the Kumaraswamy, and the Hyperbolic Distributions, for example. This range information of the shape factor, with the surprisingly simple formulas in the three above examples (Eqs. 5-8, 10, 11), can be used to select or eliminate candidate distributions for fitting. The plot of the shape factor together with plot for skewness and kurtosis can aid in setting the initial value or parameter intervals when fitting distribution to data by numerical optimization, which usually would not work well without this information.
The idea of the shape factor and the careful study of each distribution for this shape factor is the preliminary for the numerical optimization that finally finds the best fit. The information provided by shape factor plot is rough but the numerical optimization's dependence on initial value or intervals is delicate, exemplified by GB2 case. The optimization function NMinimize and FindMinimum in Mathematica sometimes can only find a local optimum at best. As shown in [21,22], the DyHF and the CMODE algorithms are the two best no-adjustment-needed global optimization algorithms. Now that the C 2 oDE algorithm is better than these two [23], it would be desirable to see how it works on the GB2 fit problem. With a foolproof universally applicable global optimization algorithm, the ado with shape factor and their boundaries will no longer be needed, or be used merely as some validations; but before that time, the hard earned knowledge about shape factor through CAS is still indispensable. This is a good topic for subsequent research.
Our shape factor idea is only a small step ahead of the skewness-kurtosis plot of Pearson [6] and McDonald et al [15,[24][25][26]. Or we just made the idea implicitly in their plot explicit. But with this clearly defined form, anyone can readily start calculating it for any interested candidate distribution.
Our formula Eq. (5) is not new, since Beta distribution has the same range of skewness, kurtosis, and shape factor as the scaled Beta distribution, the B1 distribution in [15]. Our presentation is an example of how our method can be used to easily arrive at those formulas. Theoretically equivalent expressions are not equivalent in application. With data distributions usually not having small skewness, Eq. (5) says that the Beta distribution has a shape factor roughly in the range of (1, 1.5), this not only reveals an intrinsic property of Beta distribution, but is also more easily applicable in practice than the skewness-kurtosis plot.
The residual error of all the distributions tested so far indicates that the power function or simple exponential function PDF is not enough to provide the additional freedom of shifting for the EP curve on the condition of matched first four moments. Other forms such as mixtures, combinations, or transformations of distributions may need to be considered. A previous study indicated the following transformations are good candidates [4,[27][28][29][30][31][32]: EWGU, KGG, EG, EWED, LIG, THT. Further research will be done along these lines.