Estimating Parameters of the Exponentially-Modified Logistic Distribution by the Maximum likelihood: Grey Wolf Optimization

1. Introduction

The two-parameter exponentially-modified logistic distribution is one of the new distributions that have been generalized recently. It is a distribution with more flexibility than other similar distributions for fitting data in many scientific fields, such as biological and psychological evolution, energy resource prediction, and technological and economic diffusion. Practically, it can be considered a better alternative than symmetrical distributions such as the logistic and normal distributions, as well as many other non-symmetrical statistical distributions in different application cases, especially when a little skewness exists. This can be viewed as this distribution’s main contribution. The generalization of this distribution has emerged from the importance of logistic and exponential distributions, which have both been combined to form a new, more reliable distribution called the two-parameter exponentially-modified logistic distribution. Reyes et al. first presented this distribution in 2018 [1]. The logistic distribution is recognized to be similar to the normal distribution; they are both members of the location-scale family, but the difference is that the logistic distribution has heavier tails [2]. The importance of logistic distribution is that it has the ability to be used in many scientific areas like physical science, finance, and many applications in reliability and survival analysis, besides the major utility of its distribution function in logistic regression, logit models, and neural networks [3]. More properties and details about logistic distribution are available in [4,5]. On the other hand, the exponential distribution, which is concerned with the measurement of time needed until the occurrence of a specific event [3], previously was the basis of reliability and life expectancy evaluation for many lifetime data distributions, but in further research in reliability theory, it was revealed that modeling by the exponential distribution is only useful for the first approximation and can’t be enough for a lot of problems in many cases. More details about the exponential distribution can be found in [6]. In the last two decades, in order to give better-fitting solutions and increase fit effectiveness for model functions that have no closed-form and require a numerical method in lifetime data analysis, many generalizations and modified extensions of the exponential distribution have been suggested to become more flexible and capable for modeling real-world data, especially when the characteristics of classical distributions are limited and, practically, they cannot provide a good fit in many situations [3,7,8,9,10]. Various exponentiated distributions have been generalized; for instance, in 1998, the exponentiated exponential (EE) distribution was first introduced by Gupta & Kundu [11] and is considered the first extension of the exponential distribution family. The exponentiated Weibull (EW) distribution was discovered in 2006 by Pal et al. [12], the exponentiated Gamma (EG) distribution was generalized in 2007 by Nadarajah & Gupta [13], another extension was proposed in 2011 by Nadarajah & Haghighi [14], the exponentiated log-logistic (ELL) distribution with two parameters was discovered in 2019 by Chaudhary [15], and many other well-known distributions recently have been extended and modified by the exponential distribution family. In general, there are many different statistical methodologies for estimating the parameters of any distribution, such as the maximum likelihood (ML) method, the Method of Moments (MOM), the Least Squares (LS) method, the Bayes method, and so forth. ML is the most widely used methodology among all statistical methods because of its high performance and well-known asymptotic properties for parameter estimators such as bias, consistency, efficiency, and so forth in comparison with any other method [16]. The basic principle of the ML methodology is to find the estimator values for the parameters of concern that maximize the likelihood function of the model, but in most cases, an explicit solution is rarely available because of the presence of nonlinear functions. Therefore, iterative algorithms to maximize the likelihood function are needed [17].

The primary goal of this research is to choose the best algorithm for calculating the maximum likelihood estimation of the location α and scale β parameters of the exponentially-modified logistic distribution and to show the applicability of this distribution in many areas. However, for the two-parameter exponentially-modified logistic distribution, explicit solutions to the likelihood equations were not found, and this is the main problem highlighted by this study. To solve this problem, ML estimates are found through the use of iterative numerical techniques based on traditional or non-traditional algorithms. In general, the Newton method is a common type of traditional iterative technique that is commonly used to solve the equation system generated by partial derivatives of the likelihood function to find estimated values of the parameter of interest for statistical distributions, but its major drawback is that it uses a gradient-based search algorithm to find the best parameter values based on the inverse of the hessian matrix, making it only applicable to functions that can be differentiated at least twice. At the same time, to stay in the same category and avoid such limitations, another type of traditional iterative techniques based on direct search algorithms without any need for the gradient information of the likelihood function can be used, such as the Nelder Mead (NM) algorithm. However, all traditional numerical algorithms start from a randomly selected initial guess and move towards the optimum solution iteratively, with no guarantee that the final solution is globally optimal. Besides that, the gradient-based methods can’t be used for discontinuous functions. Moreover, the final solution may get stuck at local optimum points, and the global optimum may never be reached [18]. To avoid such drawbacks, the use of non-traditional algorithms like meta-heuristic algorithms is preferable and recommended for solving complex problems, especially when traditional algorithms fail. At the same time, meta-heuristic algorithms guarantee global convergence with more simplicity, flexibility, and a derivation-free mechanism [19]. A variety of meta-heuristic algorithms have been applied in the ML method for estimating the parameters of many statistical distributions and models, such as in [20,21,22,23,24]. Among the heuristic algorithms, Particle Swarm Optimization (PSO), a well-known algorithm used in many fields [25], and new algorithms distinguished by their simplicity, flexibility, ease of implementation, and lower number of required parameters, such as the Grey Wolf Optimization (GWO) [26], the Whale Optimization Algorithm (WOA) [27], have proven to be dependable in solving real optimization problems where the objective function is nonlinear. In particular, many studies, including [28,29,30,31], considered the PSO algorithm for estimating the distribution parameters. In [32], GWO is applied for the estimation of three parameters of a new statistical distribution named the Marshall Olkin Topp Leon exponential distribution. Wang et al. [33] applied three meta-heuristic algorithms, including the GWO, PSO, and cuckoo search algorithms (CSA), and four numerical methods. The outcomes of the experimental work show that the GWO gives the most accurate and efficient results for estimating the parameters of the Weibull distribution, which is the best distribution in comparison with another distribution fitting wind energy data. Wadi [34] estimated the parameters of five statistical distributions, including Rayleigh, Weibull, Gamma, Burr Type XII, and generalized extreme value distributions, by using the ML method based on two heuristic algorithms: GWO and WOA. The results showed that the Gamma distribution based on GWO and WOA outperformed other distributions in modeling wind speed data and that GWO was more robust and faster than WOA. Furthermore, Wadi and Elmasry [35] added three more metaheuristic optimization algorithms to GWO and WOA, namely the Grasshopper Optimization Algorithm (GOA), Moth-Flame Optimization (MFO), and Salp Swarm Algorithm (SSA), which are used for estimating the parameters of the same distributions as Rayleigh and Weibull as well as other distributions such as inverse Gaussian, Burr Type XII, and Generalized Pareto to describe different wind speed data. According to the performance criteria, the Weibull distribution based on GWO, WOA, and MFO has the best goodness of fit. Al-Mhairat and Al-Quraan [36] estimated distribution parameters for Weibull, Rayleigh, and Gamma distributions by using the ML method based on three heuristic algorithms, which are PSO, GWO, and WOA, by implementing wind speed data. According to the performance indicators, such as root mean square error (RMSE) and coefficient of determination (R²), the gamma distribution based on the PSO algorithm achieved the best results, but if the mean absolute error (MAE) indicator is considered, then the method based on the GWO algorithm provides the best results. The sine-cosine algorithm (SCA) [37] is another simple, flexible, and easy-to-implement heuristic algorithm with a small number of required parameters. It is used for solving many optimization problems in various scientific areas, such as [38,39,40]. However, it has never previously been used to estimate the parameters of a statistical distribution. It can therefore be applied to the ML method to see whether or not implementing it to estimate location and scale parameters for the two-parameter exponentially-modified logistic distribution will produce better outcomes due to the features mentioned above. This paper’s contribution is to employ and examine these four meta-heuristic optimization algorithms by applying the ML method to obtain the estimator’s values for the location and scale parameters of the two-parameter exponentially-modified logistic distribution based on the GWO, WOA, SCA, and PSO algorithms, then comparing and evaluating their performances by conducting an extensive Monte-Carlo simulation study. To best of our knowledge, this is the first study to obtain the ML estimators for the location α and scale β parameters of the two-parameter exponentially-modified logistic distribution by using various meta-heuristic algorithms.

The following is how the rest of the article is organized: In Section 2, the two-parameter exponentially-modified logistic distribution and its properties are presented. In Section 3, the ML estimation methodology for GWO and the other numerical techniques based on other meta-heuristic algorithms used in this study are introduced. In Section 4, the efficiencies of the parameter estimators are compared via a comprehensive Monte-Carlo simulation study. Four applications of real datasets are implemented in Section 5. In the final section, the study ends with some conclusions.

2. Materials and Methods

2.1. Two-Parameter Exponentially-Modified Logistic Distribution

The generalization of this distribution is made by the combination of a logistic distribution with parameters for location α and scale β and the same scale parameter for the exponential distribution. As a result, the two-parameter exponentially modified logistic distribution is produced, with the left tail influenced by exponential distribution and the right tail distributed by logistic distribution. For the sake of simplicity, in the remaining portion of the study, this distribution will be denoted by "EMLOG" distribution.

If X is a random variable with a parameterized location α and scale β that follows an EMLOG distribution, X∼ EMLOG (α, β), then the probability density function (pdf) of X is:

$$\begin{gathered} f\left(x;\alpha ,\beta \right)=\frac{1}{\beta \left(e^{\left(\frac{x-\alpha }{\beta }\right)}+1\right)} \left[\left(e^{-\left(\frac{x-\alpha }{\beta }\right)}+1\right)log\left(e^{\left(\frac{x-\alpha }{\beta }\right)}+1\right)-1\right], \\ x \in R, \alpha \in R, \beta >0 \end{gathered}$$

(1)

X’s cumulative distribution function (cdf) seems to be as follows:

$$F\left(x;\alpha ,\beta \right)=1-e^{-\left(\frac{x-\alpha }{\beta }\right)}log\left(e^{\left(\frac{x-\alpha }{\beta }\right)}+1\right)$$

(2)

See the EMLOG distribution in Figure 1, where the plots of the EMLOG distribution are illustrated for certain values of α and β. The general formula for the EMLOG distribution’s k^th moment (μ_k) expression (for k = 1,2,3, ⋯) is:

$$E\left(X^k\right)=\sum _{i=0}^k\left(\begin{array}{c}k\\ j\end{array}\right)\left(2^i-2\right){\pi }^i\left|B_n\right|{\beta }^{k-i}\Gamma \left(k-i+1\right)$$

(3)

where B_n and Γ(.) refer to Bernoulli numbers and the Gamma function, respectively. The mean, variance, skewness (γ₁), and kurtosis (γ₂) values of the random variable X are computed using formula (3), see [1].

$$E\left(X\right)=\alpha +\beta$$

(4)

$$Var\left(X\right)=\left(1+\frac{{ \pi }^2}{3}\right){\beta }^2$$

(5)

$${\gamma }_1=\frac{2}{{\left(1+\frac{{ \pi }^2}{3}\right)}^{3/2}}=0.2251$$

(6)

$${\gamma }_2=\frac{\frac{7{\pi }^4}{15}+2{\pi }^2+9}{{\left(1+\frac{{ \pi }^2}{3}\right)}^2}=4.0318$$

(7)

2.2. Maximum Likelihood Estimation

The ML estimation for the parameters of interest are the values in the parameter space that maximize the likelihood function; in most cases, for calculation simplicity, the likelihood function’s logarithm is used. In this study, the log-likelihood (ln L) function is given below for estimating the unknown parameters α and β for the EMLOG distribution.

$$\begin{gathered} ln L\left(\alpha ,\beta \right)=-nlog\left(\beta \right)-\sum _{i=1}^{n}log\left(e^{ z_i}+1\right) \\ +\sum _{i=1}^{n}log\left[\left(e^{-z_i}+1\right)log\left(e^{ z_i}+1\right)-1\right] \end{gathered}$$

(8)

where z_i=(x_i-α)/β. In order to estimate the likelihood parameters for the ln L function for the EMLOG distribution, the partial derivatives with respect to the parameters of interest are taken and equated to zero. The likelihood equations are given as shown below.

$$\frac{\partial \mathit{ln}L\left(\alpha ,\beta \right)}{\partial \alpha }=\sum _{i=1}^n\frac{e^{z_i}}{e^{z_i}+1}+\sum _{i=1}^n\frac{e^{-z_i}log\left(e^{z_i}+1\right)-1}{\left(e^{-z_i}+1\right)log\left(e^{z_i}+1\right)-1}=0$$

(9)

and

$$\frac{\partial \mathit{ln}L\left(\alpha ,\beta \right)}{\partial \beta }=-n\beta +\sum _{i=1}^n\frac{\left(x_i-\alpha \right)e^{z_i}}{e^{z_i}+1} +\sum _{i=1}^n\frac{\left(x_i-\alpha \right)\left(e^{z_i}log\left(e^{z_i}+1\right)-1\right)}{\left(e^{-z_i}+1\right)log\left(e^{z_i}+1\right)-1}=0$$

(10)

As we can see from Equations (9) and (10), they have nonlinear functions, and an explicit solution for the likelihood equations cannot be obtained. Therefore, iterative algorithms are needed to solve these equations and obtain ML estimates for the location and scale. In this study, GWO, WOA, SCA, and PSO are some effective and powerful meta-heuristic algorithms considered as numerical techniques for estimating the likelihood estimators for the EMLOG distribution, and they are briefly introduced in the next few subsections.

2.2.1. Grey Wolf Optimization (GWO)

The GWO is a new intelligent method that is obtained from swarms based on meta-heuristic techniques that were modeled by the grey wolf leadership hierarchy in the process of trapping and hunting prey in nature. Mirjalili et al. made the first proposal in 2014 [26]. GWO has become a widely known critical device in swarm intelligence for optimization in almost all areas, such as engineering, physics, and many other applications in various scientific fields [41]. GWO is a straightforward population-based probabilistic algorithm motivated by the hunting and socialization of grey wolves. According to the swarm intelligence categorization, it is classified as the only algorithm for solving continuous real-life optimization problems that relies on a leadership hierarchy [42]. For the hunting process, the groups of grey wolves are categorized into four types to compose hierarchical commands. These types are called alpha (α), beta (β), delta (δ), and omega (ω), in that order. Figure 2 shows this hierarchy.

At the first level of the hierarchy, type alpha (α) represents the dominant gray wolf that makes decisions and gives orders to the other wolves in the pack. Type beta (β) represents the gray wolf that helps the alpha (α) type in making decisions and observing the movements of other wolves at the next level of the hierarchal chain. A group of alpha types will be replaced by beta types when a group of alphas dies or becomes older. Types delta (δ) and omega (ω) are the third and fourth types of gray wolves, respectively. For sure, the delta (δ) wolf type dominates the omega (ω) wolf type, and both of them represent the lowest level in the hierarchy; they are allowed to eat after the alpha (α) and beta (β) types have finished eating [43,44]. Figure 3 depicts a flow diagram of the GWO.

GWO algorithm steps: [41,44,45,46]

The GWO algorithm consists of three main phases, which are: 1) searching; 2) encircling; and 3) hunting the prey. They are summarized as follows:

For the first iteration, randomly assign the positions of the N number of wolves as well as the GWO parameters for the mathematical modeling as follows:

• Control parameter (a), which is an important parameter that declines linearly for each iteration in the range [0,2] used in this algorithm. This parameter can indeed be determined using the formula:

  $$a=2\ast \left(1-\frac{t}{T_{max}}\right)$$

  (11)

  where the iteration in progress (current) and the entire number of iterations are denoted by t and T_max, respectively.
• The coefficient vectors, A and C, can be found using the following formulas:

  $$A=2a \ast r_1-a$$

  (12)

  $$C=2 \ast r_2$$

  (13)

  where r₁ and r₂ are arbitrary vectors ranging from [0,1].

• Calculate the fitness value of each wolf type according to the fitness function, which is the same as the objective function represented by the ln L function in this study, and this fitness value refers to each wolf’s site in the pack. The highest value of the fitness function is considered the best position and assigned to the wolf of type alpha (α). The second and third highest fitness values are assigned to beta (β) and delta (δ) types of wolves, respectively. Steps 1 and 2 represent the searching phase of GWO.
• Update the position of each wolf in the pack surrounding the prey by calculating the distance between the current location (denoted by D) and the next location (denoted by) using the equations below.

  $$D=\left|C.X_p(t)-X(t)\right|$$

  (14)

  $$X(t+1)=X_p(t)-A\ast D$$

  (15)

  where $X(t)$ is the current position vector at iteration t and X_p(t) is the best solution’s position vector (optimal) when iterating to the t^th time. This step represents the phase of encircling behavior.
• Calculate the average value of the first three best solutions that refer to alpha (α), beta (β), and delta (δ) types of wolves because they have the best positions in the population. Besides that, they have the best knowledge of the prey’s potential location, which forces and obliges all the other wolves, including omega (ω), to change their current positions toward the best position, which has been determined by the following equations:

  $$X(t+1)=\frac{X_1(t)+X_2(t)+X_3(t)}{3}$$

  (16)

  where

  $$\left\{\begin{array}{c}{X}_1(t)=X_{\alpha }(t)-A_1\ast D_{\alpha }\\ X_2(t)=X_{\beta }(t)-A_1\ast D_{\beta }\\ X_3(t)=X_{\delta }(t)-A_1\ast D_{\delta }\end{array}\right.$$

  (17)

  and

  $$\left\{\begin{array}{c}{D}_{\alpha }=\left|C.X_{\alpha }(t)-X(t)\right|\\ D_{\beta }=\left|C.X_{\beta }(t)-X(t)\right|\\ D_{\delta }=\left|C.X_{\delta }(t)-X(t)\right|\end{array}\right.$$

  (18)

The fitness value of the new position is calculated, and according to that, the wolves’ alpha (α), beta (β), and delta (δ) types are updated. This step represents the hunting phase. The attack happens with respect to the changing value of the coefficient vector A, which depends on the parameter (a) and decreases from 2 to 0 for each iteration step while the algorithm is working. Therefore, A is a value generated randomly within the interval [-2a, 2a]. If |A| < 1, the wolves are ready to attack the prey. Otherwise, the wolves are forced to diverge to explore a better location.

4.

Finally, go back to step 3 to continue iteration until the convergence is satisfied by reaching the stopping criterion and the needed number of maximum iterations to overall acquire the most appropriate (optimal) solution. The solution values are called the GWO estimates of the parameters.

2.2.2. Whale Optimization Algorithm (WOA)

The WOA is also a unique swarm-based intelligent meta-heuristic methodology, recommended in 2016 by Mirjalili as well as Lewis for continuous optimization problems [27]. The inspiration for this algorithm has come from mimicking the hunting behavior of a specific type of whale (called a humpback) that applies a hunting strategy called the bubble-net feeding technique by creating bubbles along a circle around the prey, then slowly shrinking, encircling, and approaching the prey in a spiral shape through random search concerning each search agent’s location until finally, the hunt is complete [47,48,49].

WOA algorithm steps: [49]

The algorithm has three main phases that are: 1) encircling; 2) attacking by using the bubble net method, which includes shrinking encircling besides the mechanisms for spiral position updating; and 3) searching to catch prey. The steps are listed below.

• Initialize the position of the whale population randomly in the search space for the first iteration.
• Initiate the WOA parameters (a), A, and C, which are similar to GWO parameters previously calculated by Equations (11), (12), and (13), respectively, as well as other parameters such as parameter (b), which is a fixed value used to define the shape of the logarithmic spiral, and (l), a number drawn at random from the interval [–1,1]. Finally, the probability parameter (P) is set to 0.5 to give an equal chance of simulating both the shrinking surrounding and spiral approach movements of whales.
• Evaluate each whale’s fitness value in relation to the fitness function, which is considered the same as the objective function represented by the objective function ln L in this study. The best whale position in the initialized population is found and saved.
• If P < 0.5 and |A| < 1, the ongoing whale’s location is updated using the same Equations (14) and (15) as in GWO. Otherwise, if |A| > 1, one of the whales is chosen at random, and its position is updated using the following formulas:

  $$D=\left|C.X_{rand}(t)-X(t)\right|$$

  (19)

  $$X(t+1)=X_{rand}(t)-A\ast D$$

  (20)

  where X_rand is the position vector of any whale chosen at random from the current whale population.
• If P > 0.5, the current whale’s location is updated by the following formulas:

  $$D^{\prime }=\left|X_p(t)-X(t)\right|$$

  (21)

  $$X(t+1)=D^{’}{e}^{bl}cos(2\pi l)+X_p(t)$$

  (22)

  where D^’ refers to the path length between both the i^th whale and the best solution (prey) currently available.
• Verify that no whale’s updated position exceeds the search space, then go back to step 4 to continue iterating until the number of repetitions required for convergence is achieved. The solution values are called the WOA parameter estimates.

2.2.3. Sine Cosine Algorithm

The SCA is a population-based meta-heuristic technique proposed by Mirjalili in 2016, which is motivated by the mathematical trigonometric sine and cosine functions [37]. It’s been utilized to overcome a wide range of optimization issues in several areas by initializing, within the search space, a collection of a population of solutions that are iteratively assessed in relation to the objective function under the control of a set of developed optimization parameters. After that, the algorithm keeps the better solution and continuously updates it until convergence is satisfied by reaching the maximum number of iterations. This updated best position represents the best solution [50,51].

SCA algorithm steps: [52,53]

The main two phases of the SCA algorithm are 1) exploration (diversification), considered a global lookup search, and 2) exploitation (intensification), considered a local lookup search. The steps of these phases are summarized as follows:

• Initialize the position of N numbers of the population solutions randomly within the search space for the first iteration as well as the random parameters r₁, r₂, r₃, and r₄ of this algorithm, which are incorporated to strike a balance between exploration and exploitation capabilities and thus to avoid settling for local optimums. The parameter r₁ helps in determining whether an updated solution position or the movement direction of the next position is towards the best solution in the search space (r₁ < 1) or outwards from it (r₁ > 1). The r₁ parameter falls linearly from a constant (a) to 0, as seen in the equation:

  $$r_1=a-t\ast \left(\frac{a}{T_{max}}\right)$$

  (23)

The parameter r₂ is set within the interval, which helps in determining how large the extended movement of the solution towards or away from the intended target will be. The r₃ parameter is a random weight score to emphasize (r₃ > 1) or underemphasize (r₃ < 1) the significant effect of the intended target on distance calculation. The final random parameter, which is a random value defined in [0,1], can be considered a switch to choose between the trigonometric functions of sine and cosine elements.

2.

Evaluate the fitness value of each solution using using the fitness effect represented by the objective function in this study. Each fitness value refers to the position of each solution. The best (highest) value in the population is found and saved.

3.

Update the main parameters, which are r₁ by using Equation (23), and r₂, r₃, and r₄ randomly.

4.

Update the positions of all solution agents by utilizing the given equation:

$$X_i^{i+1}=\left\{\begin{array}{c}{X}_i^t+r_1\times sin(r_2)\times \left|r_3{P}_i^t-X_i^t\right|, { r}_4<0.5 \\ X_i^t+r_1\times cos(r_2)\times \left|r_3{P}_i^t-X_i^t\right|, { r}_4\ge 0.5\end{array}\right.$$

(24)

where $X_i^t$ denotes the position of the current solution in the i^th dimension at the t^th iteration and $P_i^t$ denotes the position of the target destination point in the dimension.

5.

Loop back to step 2 to continue iterating until the maximum number of iterations is reached. The solution values are called the SCA parameter estimates.

2.2.3. Particle Swarm Optimization

The PSO is considered one of the best-known population-based meta-heuristic algorithms dependent on swarm intelligence. It was proposed by Eberhart and Kennedy in 1995 [25]. PSO is a simulation of the continuous movements of particles in a swarm in a specific search area that mimics the movement behavior of bird flocks in nature using certain formulas until finally reaching the optimal solution [54]. It can be used to solve various constrained or unconstrained optimization problems, multi-objective optimization, non-linear programming, probabilistic programming, and combinatorial optimization issues [55].

PSO algorithm steps: [56,57]

• Initialize randomly the position and velocity of N number of population solutions (particles) for the first iteration as well as the algorithm parameters, which are c₁, c₂ representing acceleration coefficients, r₁, r₂ representing random numbers uniformly distributed among 0 and 1, and ω indicating the inertia weight parameter.
• Evaluate the fitness value of each solution (particle) by the fitness function ln L in this study. Each fitness value refers to the position of each solution. The best (highest) value of each particle in the population is found, compared with its previous historical movement, and then saved as a personal best solution (pbest) value. At the same time, the best value for fitness of each and every particle is found, compared with the previous historical global best, and saved as a (gbest) value.
• Update each solution’s position and velocity using the following equations:

  $$V_i^{t+1}=\omega V_i^t+c_1{r}_1\left( {pbest}_i^t-X_i^t\right)+{ c}_2{r}_2\left({gbest}^t-X_i^t\right)$$

  (25)

  $$x_i^{t+1}=x_i^t+V_i^{t+1}$$

  (26)

  where $V_i^t$ represents particle i’s velocity at iteration t, $X_i^t$ is the location of particle i during iteration t, ${pbest}_i^t$ is the best position of a particle at iteration, and ${gbes}^t$ is the most optimal (best) location of the group at iteration t.
• Loop back to step 2 again until the convergence is satisfied. The solution values are called the PSO parameter estimates.

3. Results

In this section, an extensive Monte-Carlo simulation study is carried out to compare the efficiencies of ML estimators of the model parameters for varying sample sizes, utilizing new meta-heuristic algorithms such as GWO, WOA, SCA, and the well-known PSO algorithm. All computations for the simulation study are made by Matlab R2021a software. Simulation algorithms for GWO, WOA, and SCA are coded according to Section 3, and the "particleswarm" function in the global optimization toolbox uses default initial conditions to obtain PSO estimates. Each Monte Carlo simulation run is replicated 1,000 times. The location α and scale β parameters are considered to be (α = 0, 1, 2, 3) and (β= 1, 2), respectively, for different values of sample size (n), which is taken as n = 30, 50, 100, 150, and 200. The search space (SS) for both α and β parameters is selected to be [–20,20]. Thus, 8×5×1000 = 40000 different samples are generated. The resulting estimates for location and scale parameters in the simulations are denoted by $\widehat{\alpha }$ and $\widehat{\beta }$, respectively. To analyze and evaluate the estimators’ performance, the simulated mean, bias, variance, mean square error (MSE), and deficiency (Def) values given by the Equations (27)–(31) below are used.

$$Mean \left(\widehat{\theta }\right)=\frac{\sum _i^n{\widehat{\theta }}_i}{n}$$

(27)

$$Bias\left(\widehat{\theta }\right)=E\left(\widehat{\theta }\right)-\theta$$

(28)

$$Var\left(\widehat{\theta }\right)=\frac{1}{n-1}\sum _{i=1}^n{\left(\widehat{{\theta }_i}-Mean \widehat{\theta }\right)}^2$$

(29)

$$MSE\left(\widehat{\theta }\right)=Var\left(\widehat{\theta }\right)+{\left(Bias\left(\widehat{\theta }\right)\right)}^2$$

(30)

$$Def\left(\widehat{\alpha } , \widehat{\beta }\right)=MSE\left(\widehat{\alpha }\right)+MSE\left(\widehat{\beta }\right)$$

(31)

where, θ = (α, β) ∈ R×R⁺. The resulting simulated values of mean, bias, MSE, and Def for $\widehat{\alpha }$ and $\widehat{\beta }$ are given in Table 1, Table 2, Table 3 and Table 4. The simulated values show that the GWO provides the best results in comparison with the other algorithms. Based on the simulated bias results, we conclude that the GWO algorithm’s estimator values for both parameters α and β provide the smallest bias values for almost all sample sizes when the true value of β = 1. Otherwise, when β = 2, the smallest bias values in most cases belong to the PSO algorithm. Concerning WOA, it shows the lowest bias values (which are the same as GWO values) in some cases and higher values in others, which indicates that it’s not a stable method for the EMLOG model. Anyway, in all cases, SCA demonstrates the worst performance with the largest bias values. Concerning MSE values, it’s very clear that MSE values for GWO estimators outperform other methods for all n values. Also, it can be noticed from the Table 1, Table 2, Table 3 and Table 4 that the MSE values of α and β for WOA estimators give almost the same values as GWO estimators in particular cases, but they do not show the best performance for all samples. These cases are classified into two categories: Firstly, when β = 1, as the following: for 1) α = 0, and n = 30, 100 2) α = 1, and n = 150. 2) α = 2, n = 30, 100, 150, 200. 3) α = 3 and n = 50. Secondly, when β = 2, as follows: 1) for α = 1, and n = 50, 100, 200. 2) α = 2 with n = 30, 50, 100, and 150. 3) For all n values, and α must be 0 either 3. In these two categories, the WOA is considered the same as the GWO for computing and estimating both the location α and scale β parameters that provide very high performance, which is why it is so close to the GWO with respect to MSE values. Except for these two categories, the PSO is considered the second-best method after the GWO with respect to MSE values for estimating only the location parameter α. Also, it’s clear that for all n values, SCA gives the lowest performance with the largest MSE values for both α and β parameters. In terms of the Def criterion, the strongest performance with the lowest deficiency for all n values is demonstrated by GWO. In light of that, it can be said that the ML estimator values of and β parameters using GWO give the most efficient values among the others. When the two categories mentioned before are considered, we see that the WOA is as efficient as the GWO method concerning deficiency values. Otherwise, out of these two categories, the PSO will be the second-best method after the GWO for estimating both α and β parameters. The weakest performance with the highest deficiency values is provided by SCA in all cases. Based on these findings, we can conclude that the GWO yields the best results when calculating ML estimators for unknown EMLOG distribution parameters (α, β). Therefore, GWO is more recommended and preferable than any other meta-heuristic algorithm used in this study. This is consistent with the findings in the literature [33,34,35,36], which allows us to draw the conclusion that the GWO algorithm is a powerful meta-heuristic algorithm even when used with various performance criteria or statistical distribution types. The simulation results also indicate that the SCA algorithm is unsuitable for estimating the location and scale parameters of this distribution due to the high values of its deficiency criterion.

4. Applications

Four real datasets are modeled using the EMLOG distribution in this section to demonstrate the utility of this model in practice for the applications. The first two datasets, which are used in various engineering areas, show a highly symmetrical shape, and the modeling performance of the EMLOG distribution for these two datasets is compared with the other well-known and commonly used distributions in the statistical literature, such as Gamma, log-normal, log-logistic, and others. In the last two datasets, which have a skewed shape and are used in medical and banking fields, the EMLOG distribution’s performance is compared to that of the normal and logistic distributions, which are symmetrical and belong to the location-scale family, to examine the flexibility and better goodness of fit of the EMLOG distribution among these two distributions. The modeling performance of the EMLOG distribution has been compared using well-known different criteria, including log-likelihood values, the Akaike Information Criterion (AIC), the corrected AIC (AICc), the consistent AIC (CAIC), the Bayesian Information Criterion (BIC), and the Hannan-Quinn Information Criterion (HQIC). For extra information on these criteria and their implementation, see [58,59,60]. The mathematical expression of these criteria is given by

$$AIC=2P-2 ln L$$

(32)

$$AI{C}_c=AIC+\frac{2P\left(P+1\right)}{n-P-1}$$

(33)

$$CAIC=P [ln(n)+1]-2 ln L$$

(34)

$$BIC=P ln(n)-2 ln L$$

(35)

$$HQIC=2 P ln [ln(n)]-2 ln L$$

(36)

where ln L is the maximized likelihood function, the number of observations is n, and P is just the total number of model parameters. The model of probability is considered the best-fit model when it has lower values for these mentioned criteria in comparison with other probability distributions [61].

4.1. Dataset 1: Tensile Tensile Strength of 69 Carbon Fibers

This dataset contains the tensile strength (in GPa) of 69 carbon fibers evaluated under stress at 20 mm gauge lengths. This dataset was originally used for the first time in 1982 by Bader and Priest [62]. The data are given as: 1.312, 1.314, 1.479, 1.552, 1.700, 1.803, 1.861, 1.865, 1.944, 1.958, 1.966, 1.997, 2.006, 2.021, 2.027, 2.055, 2.063, 2.098, 2.140, 2.179, 2.224, 2.240, 2.253, 2.270, 2.272, 2.274, 2.301, 2.301, 2.359, 2.382, 2.382, 2.426, 2.434, 2.435, 2.478, 2.490, 2.511, 2.514, 2.535, 2.554, 2.566, 2.570, 2.586, 2.629, 2.633, 2.642, 2.648, 2.684, 2.697, 2.726, 2.770, 2.773, 2.800, 2.809, 2.818, 2.821, 2.848, 2.880, 2.954, 3.012, 3.067, 3.084, 3.090, 3.096, 3.128, 3.233, 3.433, 3.585, 3.858. Descriptive statistics, including the values of sample size (n), minimum (min), first quartile (1^st Qu.), mean, mode, median, third quartile (3^rd Qu.), maximum (max), variance (S²), skewness (γ₁), and kurtosis (γ₂) coefficients, respectively, are given by Table 5, which demonstrates a high degree of symmetry for the data.

Various well-known distributions such as Gamma, log-normal, Log-Logistic, Rayleigh, Weibull, and exponential are used for modeling this data; see [63]. A comparison between EMLOG distribution modeling performance and the mentioned commonly used distributions according to the lnL, AIC, AICc, CAIC, BIC, and HQIC criteria is given in Table 6.

The results in Table 6 show that the EMLOG distribution performance gives a better fit than its rivals in terms of the considered criteria.

4.2. Dataset 2: Strengths of Glass Fibers.

The Strength of Glass Fibers dataset was introduced for the first time by Smith and Naylor [64]. It’s made up of 63 observations about the strengths of glass fibers and is widely used in the statistical literature. The data are as follows: 0.55, 0.93, 1.25, 1.36, 1.49, 1.52, 1.58, 1.61, 1.64, 1.68, 1.73, 1.81, 2.00, 0.74, 1.04, 1.27, 1.39, 1.49, 1.53, 1.59, 1.61, 1.66, 1.68, 1.76, 1.82, 2.01, 0.77, 1.11, 1.28, 1.42, 1.5, 1.54, 1.6, 1.62, 1.66, 1.69, 1.76, 1.84, 2.24, 0.81, 1.13, 1.29, 1.48, 1.5, 1.55, 1.61, 1.62, 1.66, 1.7, 1.77, 1.84, 0.84, 1.24, 1.3, 1.48, 1.51, 1.55, 1.61, 1.63, 1.67, 1.7, 1.78, 1.89. Table 7 shows the descriptive statistics for this dataset, which present a high symmetrical degree for the data. The resulting values of fitting this data to the EMLOG distribution model in comparison with Gamma, lognormal, log-logistic, Rayleigh, and Exponential distributions in terms of the criteria that were chosen are obtained in Table 8.

The results show that the EMLOG model outperforms any competitor distribution in terms of modeling performance.

4.3. Dataset 3: Bladder Cancer Patients.

This dataset is a biologically uncensored univariate dataset that reflects the remission periods (in months) of 128 bladder cancer patients who were randomly selected by Lee and Wang in 2003 [65]. This dataset is frequently cited in the literature; for more information, see [66,67,68,69,70,71,72]. The data are given as: 0.08, 2.09, 3.48, 4.87, 6.94, 8.66, 13.11, 23.63, 0.20, 2.23, 3.52, 4.98, 6.97, 9.02, 13.29, 0.40, 2.26, 3.57, 5.06, 7.09, 9.22, 13.80, 25.74, 0.50, 2.46, 3.64, 5.09, 7.26, 9.47, 14.24, 25.82, 0.51, 2.54, 3.70, 5.17, 7.28, 9.74, 14.76, 26.31, 0.81, 2.62, 3.82, 5.32, 7.32, 10.06, 14.77, 32.15, 2.64, 3.88, 5.32, 7.39, 10.34, 14.83, 34.26, 0.90, 2.69, 4.18, 5.34, 7.59, 10.66, 15.96, 36.66, 1.05, 2.69, 4.23, 5.41, 7.62, 10.75, 16.62, 43.01, 1.19, 2.75, 4.26, 5.41, 7.63, 17.12, 46.12, 1.26, 2.83, 4.33, 5.49, 7.66, 11.25, 17.14, 79.05, 1.35, 2.87, 5.62, 7.87, 11.64, 17.36, 1.40, 3.02, 4.34, 5.71, 7.93, 11.79, 18.10, 1.46, 4.40, 5.85, 8.26, 11.98, 19.13, 1.76, 3.25, 4.50, 6.25, 8.37, 12.02, 2.02, 3.31, 4.51, 6.54, 8.53, 12.03, 20.28, 2.02, 3.36, 6.76, 12.07, 21.73, 2.07, 3.36, 6.93, 8.65, 12.63, 22.69. The positive skewness of this dataset in terms of fundamental descriptive statistics values is given in Table 9 below.

The EMLOG distribution’s modeling performance is compared to that of the normal and logistic distributions. The results given in Table 10 show that the EMLOG distribution has a better fit in comparison with the normal and logistic distributions in terms of the considered criteria.

4.4. Dataset 4: Waiting Times (in Minutes) of 100 Bank Customers.

This dataset contains 100 observations that indicate the waiting periods (in minutes) before service for 100 bank customers, as examined and described by Ghitany et al. [73]. The data are given as:

0.8, 0.8, 1.3, 1.5, 1.8, 1.9, 1.9, 2.1, 2.6, 2.7, 2.9, 3.1, 3.2, 3.3, 3.5, 3.6, 4.0, 4.1, 4.2, 4.2, 4.3, 4.3, 4.4, 4.4, 4.6, 4.7, 4.7, 4.8, 4.9, 4.9, 5.0, 5.3, 5.5, 5.7, 5.7, 6.1, 6.2, 6.2, 6.2, 6.3, 6.7, 6.9, 7.1, 7.1, 7.1, 7.1, 7.4, 7.6, 7.7, 8.0, 8.2, 8.6, 8.6, 8.6, 8.8, 8.8, 8.9, 8.9, 9.5, 9.6, 9.7, 9.8, 10.7, 10.9, 11.0, 11.0, 11.1, 11.2, 11.2, 11.5, 11.9, 12.4, 12.5, 12.9, 13.0, 13.1, 13.3, 13.6, 13.7, 13.9, 14.1, 15.4, 15.4, 17.3, 17.3, 18.1, 18.2, 18.4, 18.9, 19.0, 19.9, 20.6, 21.3, 21.4, 21.9, 23.0, 27.0, 31.6, 33.1, 38.5.

Descriptive statistics for this dataset are given in Table 11. From Table 12, we can see that the EMLOG distribution seems to provide a better fit than normal and logistic distributions according to the considered criteria.

It can be noticed that in the first and second datasets, which have a symmetrical shape, the EMLOG model outperforms unsymmetrical models such as the Gamma, lognormal, log-logistic, Rayleigh, Weibull, and exponential distributions, whereas in the third and fourth datasets, which have a skewed shape, the EMLOG model provides a better fit than symmetrical models like the normal and logistic distributions. All this can lead us to say that the EMLOG model can give a better fit than many popular distributions like Gamma, Weibull and others when the data is symmetrical. On the other hand, the EMLOG model can give a better fit than normal and logistic distributions in many cases when the data is skewed. This shows the extent of the distinction of this model and its suitability for many practical cases and applications.

5. Discussion and Conclusions

The EMLOG distribution is a new distribution obtained by combining a logistic distribution with an exponential distribution. This distribution can be used as an alternative to logistic and normal distributions with a superior fit in many situations, and it can outperform many prominent distributions such as Gamma, Weibull, and others. Because EMLOG distribution is more flexible than its alternative distributions, it has been used in recent years in sectors such as technology, energy, marketing, biology, psychology, and so on. Therefore, parameter estimation for this distribution is very important. In this study, ML estimates of the EMLOG distribution’s location α and scale β parameters are investigated, which cannot be obtained explicitly because of the complication of finding a solution for their nonlinear likelihood equations. In such cases, numerical methods are needed. The ML estimation method based on meta-heuristic algorithms such as GWO, WOA, SCA, and PSO is considered a wise alternative and more convenient than any other traditional methods. A Monte Carlo simulation study was conducted to compare the performance of GWO estimator values in comparison to other meta-heuristic techniques with respect to bias, MSE, and Def criteria. According to the simulation results, the ML estimates of GWO show the best performance in comparison with the WOA, SCA, and PSO methods.