A Type I Generalized Logistic Distribution: Solving Its Estimation Problems with a Bayesian Approach and Numerical Applications Based on Simulated and Engineering Data

Abstract

The family of logistic type distributions has been widely studied and applied in the literature. However, certain estimation problems exist in some members of this family. Particularly, the three-parameter type I generalized logistic distribution presents these problems, where the parameter space must be restricted for the existence of their maximum likelihood estimators. In this paper, motivated by the complexities that arise in the inference under the likelihood approach utilizing this distribution, we propose a Bayesian approach to solve these problems. A simulation study is carried out to assess the performance of some posterior distributional characteristics, such as the mean, using Monte Carlo Markov chain methods. To illustrate the potentiality of the Bayesian estimation in the three-parameter type I generalized logistic distribution, we apply the proposed method to real-world data related to the copper metallurgical engineering area.

1. Introduction

There is a wide variety of probabilistic models that have been proposed to analysis data generated from continuous random variables [1,2]. One of these models is the logistic distribution [3]. Methods for estimating the parameters of some of these distributions have been well studied for different dimensionality of the parameter vector [4].

The logistic distribution has been largely studied and generalized from different perspectives; see, for example, the works proposed in [5,6,7,8]. One of these models is the type I generalized logistic (IGL) distribution, which can be considered as a family of proportional reversed hazard functions with the logistic distribution being its base. The IGL distribution has been used to model data with an unimodal probability density function (PDF). For a comprehensive account of the theory and applications of logistic distributions, we refer the readers to [2,3,9].

A three-parameter IGL distribution was proposed in [10] and methods of estimation for its parameters were also presented in this work. The authors demonstrated that its log-likelihood function is strictly decreasing in one of its parameters, and maximized only when these parameters tend to infinity. Then, its maximum likelihood (ML) estimators do not exist and a hybrid method must be employed to optimize the log-likelihood function by plugging it at the posterior mean, but at the same time solving the problem of non-existence.

A technique for estimating the location and scale parameters of the IGL distribution was proposed in [11] from a U-statistic constructed by using the best linear functions of order statistics as kernels. An efficiency comparison of the proposed estimators with respect to the ML estimator was also made by these authors. In a comparative study, the authors considered that there is an advantage when using a U-statistic for estimating the parameters from a practical point of view.

Some empirical evidence, based on stock markets and data sets of financial expected values, suggests that the true distribution of such expected values can be positively skewed in times of high inflation and more peaked than the normal distribution. The logistic distribution is often used in econometrics to model the consumer inflation rates due to its heavy tails compared to the normal distribution [12]. Its generalizations, which additionally allow for asymmetries, have been applied to model extreme risks in the context of the German stock market [13], and has been proved to be a simpler model than its harmonized index of consumer prices in inflation rates [14].

A method for determining explicit estimators of location and scale parameters was proposed in [15] by approximating the likelihood function using the Tiku method [16]. An application in genomic experiments utilizing the IGL distribution to determine whether a gene is differentially expressed under different conditions was presented in [17]. The authors simplified the computational complexity carrying out likelihood ratio (LR) tests for several thousands of genes. Moreover, an approximate LR test was also proposed to generalize the two-class LR method into multi-class microarray data. The IGL PDF can be symmetrical, left-skewed, right-skewed, and reversed-J shaped, which provides us with high flexibility in the modeling of the data.

The literature on the study of the performance of the joint estimators of the three-parameter IGL distribution is limited. Particularly for the ML estimators because, as mentioned, an estimation problem exists and the corresponding parameter space must be restricted for their existence. Despite this problem, an R package named glogis [18] was published without further discussion about this problem [19].

We can mention some other commonly used classical (or frequentist) estimation methods in addition to the ML approach. For example, we have the ordinary and weighted least squares methods [20], which are defined in terms of order statistics, making these methods robust and helpful for estimating the parameters of a skew distribution. Also, other estimation approaches related to percentiles, minimum distances, and bootstrapping can be considered [21]. Furthermore, the probability-weighted moment (PWM) method may be mentioned, which has been investigated by many researchers. This method was originally proposed in [22], is commonly employed in theoretical and empirical studies, and has determined frameworks that are more stable against outliers. The PWM method may be generally used to estimate parameters of a distribution whose inverse form cannot be formulated explicitly. The L-moment method was introduced in [23], which is linear combination of the PWM method. The L-moment method [24] is more easily related to the distribution shape and spread than the PWM method. However, when estimating parameters, there is no reason to distinguish between the L-moment and PWM methods, because they lead to identical parameter estimates, at least for their usual settings. An R package named lmom [25] has been implemented for the L-moment method.

When a distribution can be written as a mixture model [26,27], we may use the expectation-maximization algorithm proposed in [28] to efficiently estimate its parameters. Nevertheless, this algorithm may present some disadvantages at the step of maximization due to the dependence on initial values for multimodal likelihood functions [29]. To solve this and other issues discussed in the following sections, an alternative method based on a Bayesian approach was proposed in [30], which is known as the stochastic expectation maximization algorithm. This algorithm includes a step of stochastic simulation based on the Gibbs sampling and the Metropolis–Hastings technique from the posterior distribution of the parameters, that is, it it takes prior information and utilizes it to deal with such disadvantages.

An approach based on importance sampling was used in [31] to estimate the shape and scale parameters of the generalized exponential distribution, in addition to the Gibbs and Metropolis samplers to predict the behavior of the distribution. To the best of our knowledge, this approach has been not employed for the IGL distribution.

The main objective of this investigation is to use a Bayesian approach to obtain reliable estimators of the three-parameter IGL distribution. Specifically, we apply a Monte Carlo Markov chain (MCMC) algorithm to guarantee the coverage of the stationary distribution support. Our secondary objectives are: (i) to conduct a simulation study to assess the performance of some posterior distributional characteristics, such as the mean, utilizing Monte Carlo Markov chain methods; and (ii) to apply the proposed method to real-world data related to engineering area.

The outline of the paper is as follows. Section 2 presents the three-parameter IGL distribution. In Section 3, we introduce a classical approach to estimate the parameters of this distribution, whereas in Section 4, we provide the numerical approximation of the posterior distributions to make Bayesian inference by using MCMC methods. In Section 5, Monte Carlo simulation results are reported and discussed. Section 6 applies the obtained results to a real-world data set from the area of copper metallurgical engineering to illustrate the potentiality of the Bayesian analysis in the three-parameter IGL distribution. Finally, in Section 7, we sketch the conclusions of this study.

2. The Type I Generalized Logistic Distribution and Related Distributions

Let the random variable W follow an IGL distribution with location ($\mu$), scale ($\sigma$), and shape (b) parameters. Then, the PDF and cumulative distribution function of W are, respectively, given by

$$f_W(w;\mu ,b,\sigma )=\frac{b}{\sigma }\frac{{\mathrm{e}}^{-(w-\mu )/\sigma }}{{\left(1+{\mathrm{e}}^{-(w-\mu )/\sigma }\right)}^{b+1}},F_W(w;\mu ,b,\sigma )={\left(1+{\mathrm{e}}^{-(w-\mu )/\sigma }\right)}^{-b},$$

(1)

with $w\in \mathbb{R},\mu \in \mathbb{R},b>0$, $\sigma >0$, and “e” denoting the exponential function. The family of IGL distributions consists of members with their positively skewed PDF and coefficient of kurtosis greater than that of the logistic distribution. If W is an IGL distributed random variable, then $-W$ has a type II generalized logistic (IIGL) distribution [32]. Additional generalizations of the logistic distribution are discussed in [1]. For $b=1$, the distribution is symmetric, for $b<1$, it is skewed to the left (negative asymmetry), and for $b>1$, it is skewed to the right (positive asymmetry).

The moment generating function of the IGL distribution is stated as

$${\mathcal{M}}_W\left(r\right)={\mathrm{e}}^{\mu r}\frac{\Gamma (1-\sigma r)\Gamma (b+\sigma r)}{\Gamma \left(b\right)},-b<r<1,$$

(2)

where $\Gamma$ is the usual gamma function defined as $\Gamma \left(z\right)={\int }_0^{\infty }{x}^{z-1}{\mathrm{e}}^{-x}\mathrm{d}x$, for $z>0$. From (2), it follows that

$$\mathrm{E}\left(W\right)=\mu +\sigma (\psi \left(b\right)-\psi \left(1\right)),\mathrm{Var}\left(W\right)={\sigma }^2({\psi }^{{}^{\prime }}\left(1\right)+{\psi }^{{}^{\prime }}\left(b\right)),$$

(3)

with $\psi \left(z\right)={\Gamma }^{{}^{\prime }}\left(z\right)/\Gamma \left(z\right)$ and ${\psi }^{{}^{\prime }}\left(z\right)=\mathrm{d}\psi \left(z\right)/\mathrm{d}z$ being the digamma and trigamma functions, respectively. From the formulas established in (3), the coefficient of skewness for W, corresponding to the third standardized moment, is formulated as

$$\mathrm{Skew}\left(W\right)=\frac{{\psi }^{{}^{″}}\left(b\right)-{\psi }^{{}^{″}}\left(1\right)}{{({\psi }^{\prime }\left(b\right)+{\psi }^{\prime }\left(1\right))}^{3/2}}=:{\mathrm{Skew}}_W\left(b\right).$$

(4)

Note that the expression given in (4) does not depend on the parameter $\sigma$.

For the sake of comparison of how flexible a generalized logistic distribution can be, when incorporating a fourth parameter, in relation to the IGL case, we introduce the type IV generalized logistic (IVGL) distribution with parameters $\mu$, $\sigma$, $\alpha$ and $\beta$. If $M\sim \mathrm{IVGL}(\alpha ,\beta )$, then we get

$$\mathrm{Skew}\left(M\right)=\frac{{\psi }^{{}^{″}}\left(\alpha \right)-{\psi }^{{}^{″}}\left(\beta \right)}{{({\psi }^{\prime }\left(\alpha \right)+{\psi }^{\prime }\left(\beta \right))}^{3/2}}=:{\mathrm{Skew}}_M(\alpha ,\beta ).$$

(5)

In Figure 1a, we show the effect of the shape parameter b for the IGL distribution, when the other parameters $\mu$ and $\sigma$ take different values. It is clear from this figure that the distribution is negatively skewed for $b<1$ and is positively skewed for $b>1$. The IGL distribution is flexible in the sense that the rate of change of the coefficient of skewness is greater for negative values than for the positive values in relation to the other generalized logistic distributions. In Figure 1b, we show the flexibility to capture asymmetry of the IVGL distribution, which is more pronounced for values of the parameter $\beta$, approximately when $\beta <0.4$. Thus, this additional parameter of the IVGL distribution is not relevant. For the IIGL distribution, the skewness is a decreasing function of b. Note that if the value of b tends to infinity, then the IIGL distribution has heavier tails. We observe that the coefficient of skewness is approximately constant for values of the shape parameter greater than 0.025. Therefore, the IGL distribution can be used in robustness studies of classical procedures, since extreme values are often observed in real-world data.

3. Inference: Classical Approaches

Next, we describe the moment method of estimation. The expressions stated in (3) are useful to find the moment estimators of $\mu$, $\sigma$ and b, ${\widehat{\mu }}_{\mathrm{M}}$, ${\widehat{\sigma }}_{\mathrm{M}}$ and ${\widehat{b}}_{\mathrm{M}}$ namely, respectively. The estimator (moment-based) of Skew$\left(W\right)$ for a sample $W_n={(W_1,\dots ,W_n)}^{\top }$ from the IGL$(\mu ,\sigma ,b)$ distribution is given by

$$\widehat{\mathrm{Skew}}\left(W_n\right)={\displaystyle \frac{{(n-3)}^{-1}{\sum }_{i=1}^n{(W_i-{\overline{W}}_n)}^3}{{\left({(n-1)}^{-1}{\sum }_{i=1}^n{(W_i-{\overline{W}}_n)}^2\right)}^{3/2}}}.$$

(6)

If $-2<\widehat{\mathrm{Skew}}\left(W_n\right)<1.1396$, with $\widehat{\mathrm{Skew}}\left(W_n\right)$ being defined in (6), then it is possible to obtain a unique solution ${\widehat{b}}_{\mathrm{M}}$ from a condition stated as

$${\displaystyle \frac{{(n-3)}^{-1}{\sum }_{i=1}^n{(W_i-{\overline{W}}_n)}^3}{{\left({(n-1)}^{-1}{\sum }_{i=1}^n{(W_i-{\overline{W}}_n)}^2\right)}^{3/2}}}={\displaystyle \frac{{\psi }^{{}^{″}}\left({\widehat{b}}_{\mathrm{M}}\right)-{\psi }^{{}^{″}}\left(1\right)}{{\left({\psi }^{{}^{\prime }}\left({\widehat{b}}_{\mathrm{M}}\right)+{\psi }^{{}^{\prime }}\left(1\right)\right)}^{3/2}}}.$$

(7)

Thus, over the condition given in (7), we obtain the moment estimators for $\mu$ and $\sigma$ by analytically solving the equations established as

$$\begin{array}{c}\hfill {\displaystyle \frac{{\sum }_{i=1}^n{(W_i-{\overline{W}}_n)}^2}{n-1}}={\widehat{\sigma }}_{\mathrm{M}}^2({\psi }^{{}^{\prime }}\left(1\right)+{\psi }^{{}^{\prime }}\left({\widehat{b}}_{\mathrm{M}}\right)),{\overline{W}}_n={\widehat{\mu }}_{\mathrm{M}}+{\widehat{\sigma }}_{\mathrm{M}}(\psi \left({\widehat{b}}_{\mathrm{M}}\right)-\psi \left(1\right)).\end{array}$$

Moreover, in [10], note that $({\widehat{\mu }}_{\mathrm{M}},{\widehat{\sigma }}_{\mathrm{M}},{\widehat{b}}_{\mathrm{M}})$ are consistent and asymptotically unbiased for ${(\mu ,\sigma ,b)}^{\top }$ when ${\widehat{b}}_{\mathrm{M}}$ exists. However, from a practical numerical point of view, we see in our simulation study that the degree of sample asymmetry and the sample size affect the existence of ${\widehat{b}}_{\mathrm{M}}$, without being sure of the restrictions for $\widehat{\mathrm{Skew}}\left(W_n\right)$, and so ${\widehat{b}}_{\mathrm{M}}$ does not satisfy the condition given in (7).

Next, we consider the ML method of estimation. Assuming $W_n$ is a sample of size n from the IGL$(\mu ,\sigma ,b)$ distribution, the log-likelihood function of $\theta ={(\mu ,\sigma ,b)}^{\top }$ is formulated as $\ell :=\ell (\mu ,\sigma ,b)$, where

$$\ell (\mu ,\sigma ,b)=nlog\left(b\right)-nlog\left(\sigma \right)-\frac{1}{\sigma }{\sum }_{i=1}^n(W_i-\mu )-(b+1){\sum }_{i=1}^{n}log\left(1+{\mathrm{e}}^{-(W_i-\mu )/\sigma }\right).$$

(8)

From (8), we obtain a closed-form expression for estimating b by the ML method, that is, we reach

$$\widehat{b}(\mu ,\sigma ):=\widehat{b}={\overline{T_n}}^{-1},$$

(9)

with ${\overline{T}}_n(\mu ,\sigma ,W_n):={\overline{T}}_n=(1/n){\sum }_{i=1}^{n}log(1+{\mathrm{e}}^{-(W_i-\mu )/\sigma })$. When plugging the estimator $\widehat{b}$ into the log-likelihood function, we get the profile log-likelihood function, denoted in general by ${\ell }^{*}(\theta ,g\left(\theta \right))$ (or concentrated log-likelihood function, whose terminology was borrowed from the ML theory, where the insertion of certain partial solutions $\varphi =g\left(\theta \right)$ leads to such a concentration). Hence, the corresponding profile log-likelihood function is defined as

$${\ell }^{*}(\mu ,\widehat{b},\sigma )=-n{\overline{T}}_n-nlog\left(n{\overline{T}}_n\right)+H(\mu ,\sigma ,W_n),$$

(10)

where H is a function that does not depend on $\widehat{b}$ given by

$$H(\mu ,\sigma ,W_n)=n(log\left(n\right)-log\left(\sigma \right)-1)-n{\overline{Y}}_n,$$

(11)

with ${\overline{Y}}_n=(1/n){\sum }_{i=1}^n(W_i-\mu )/\sigma$. Therefore, with a usual optimization algorithm to maximize the expression stated in (10), we find the ML estimate $\widehat{\sigma }$. However, the ML estimation in the IGL distribution has some problems. In [10], it was shown that there is a path in the parameter space along which the likelihood function becomes unbounded.

We present the Fisher information matrix (FIM) of $\theta$ in order to use it as a guide to define priori distributions for Bayesian inference on $\theta ={(\mu ,\sigma ,b)}^{\top }={({\theta }_1,{\theta }_2,{\theta }_3)}^{\top }$. The FIM of $\theta$, say $\mathit{K}\left(\theta \right)$, have elements given by

$$-{\kappa }_{rs}=-n\mathrm{E}\left(\frac{{\partial }^2\ell }{\partial {\theta }_r\partial {\theta }_s}\right),r,s\in (1,2,3),$$

with ℓ being the corresponding log-likelihood function. Thus, the FIM of $\theta$ is stated as $\mathit{K}\left(\theta \right)=n{\mathit{K}}_1\left(\theta \right)$, where ${\mathit{K}}_1\left(\theta \right)$ is the FIM when $n=1$, that is,

$${\mathit{K}}_1\left(\theta \right)=\mathrm{E}\left(\frac{\partial \ell }{\partial {\theta }_r}\frac{\partial \ell }{\partial {\theta }_s}\right)=-\mathrm{E}\left(\frac{{\partial }^2\ell }{\partial {\theta }_r\partial {\theta }_s}\right),r,s\in (1,2,3).$$

In particular, we have

$$\begin{array}{ccc}\hfill {\kappa }_{12}& =& -\frac{n}{\sigma (b+1)},{\kappa }_{13}=-\frac{nb}{{\sigma }^2(b+2)}(\psi (b+1)-\psi \left(2\right)),{\kappa }_{11}=-\frac{nb}{{\sigma }^2(b+2)},\hfill \\ \hfill {\kappa }_{22}& =& -\frac{n}{b^2},{\kappa }_{23}=-\frac{n}{\sigma (b+1)}(\psi \left(b\right)-\psi \left(2\right)),\hfill \\ \hfill {\kappa }_{33}& =& -\frac{n}{{\sigma }^2}\left(1+\frac{b}{b+2}\left({\psi }^{\prime }(b+1)+{\psi }^{\prime }\left(2\right)+{(\psi (b+1)-\psi \left(2\right))}^2\right)\right).\hfill \end{array}$$

(12)

Note that the FIM is independent of $\mu$ so that we can write ${\mathit{K}}_1\left({\theta }^{★}\right):={\mathit{K}}_1\left(\theta \right)$, with ${\theta }^{★}=(\sigma ,b)$. Moreover, we observe that, if $b=1$, difficulties appear when applying the Fisher asymptotic theory, having a singular FIM. Additional details can be found in [33].

4. Inference: A Bayesian Approach

One of the first steps to do Bayesian statistics is generating information in terms of a PDF. Obtaining a priori distribution for a single parameter can be simple, if some experience in similar situations has been established. In the Bayesian inference, the posterior PDF $p\left(\theta \right|W_n)$ is constructed from the likelihood function $p\left(W_n\right|\theta )$, and the prior PDF $p\left(\theta \right)$ for the parameters is obtained by means of what is known as the Bayes rule, that is,

$$p\left(\theta \right|W_n)=\frac{1}{p\left(W_n\right)}p\left(W_n\right|\theta )p\left(\theta \right)\propto p\left(W_n\right|\theta )p\left(\theta \right)=:p^{★}\left(\theta \right|W_n).$$

(13)

In this context, $p\left(W_n\right)$ is known as the marginal likelihood function and is given by

$$p\left(W_n\right)={\int }_{\theta }p\left(W_n\right|\theta )\mathrm{d}\theta ,$$

(14)

which is often extremely hard to calculate when the likelihood function has a complex structure, as is the case when considering the IGL distribution. Hence, the Bayesian estimator of any function g of $\theta$ under squared error loss function is defined as

$$\widehat{g}\left(W_n\right):=E\left(g\left(\theta \right)\right)=\frac{1}{p\left(W_n\right)}{\int }_{\theta }g\left(\theta \right)p\left(W_n\right|\theta )p\left(\theta \right)\mathrm{d}\theta$$

(15)

and under absolute error loss function, it is the median of $g\left(\theta \right)$, that is,

$$\tilde{g}\left(W_n\right):=Med\left(g\left(\theta \right)\right).$$

(16)

Next, we conduct the prior specification. For a parameter vector, establishing a priori joint distribution can be much more complicated. It is also difficult to obtain information on the dependence of the parameters and express it in terms of a joint distribution.

We can consider a priori distribution for the parameter vector $\theta$, say $\theta ={({\theta }_1,{\theta }_2,{\theta }_3)}^{\top }$, as $p\left(\theta \right)=p\left({\theta }_1\right|{\theta }_2,{\theta }_3)\times p\left({\theta }_2\right|{\theta }_3)\times p\left({\theta }_3\right)$. Thus, by considering that each parameter is independent of the others, the last expression remains as $p\left(\theta \right)=p\left({\theta }_1\right)\times p\left({\theta }_2\right)\times p\left({\theta }_3\right)$.

The posterior PDF of $\theta ={(\mu ,\sigma ,b)}^{\top }$ is, up to the constant factor $p\left(W_n\right)$, stated as

$$p^{★}\left(\theta \right|W_n)=\frac{b^n}{{\sigma }^n}\frac{{\mathrm{e}}^{-S_n/\sigma +\mu /\sigma }}{{\sum }_{i=1}^n{\left(1+{\mathrm{e}}^{-(W_i-\mu )/\sigma }\right)}^{b+1}}p\left(\theta \right),$$

(17)

where

$$p\left(W_n\right)={\int }_0^{+\infty }{\int }_0^{+\infty }{\int }_{-\infty }^{+\infty }\frac{b^n}{{\sigma }^n}\frac{{\mathrm{e}}^{-S_n/\sigma +\mu /\sigma }}{{\sum }_{i=1}^n{\left(1+{\mathrm{e}}^{-(W_i-\mu )/\sigma }\right)}^{b+1}}\mathrm{d}\mu \mathrm{d}\sigma \mathrm{d}b,$$

(18)

with $S_n:={\sum }_{i=1}^n{W}_i$ and $p\left(\theta \right)$ acquiring the structures specified above (appropriate for its application). Therefore, the posterior PDF of $\theta ={(\mu ,\sigma ,b)}^{\top }$ given the data is formulated as

$$p(\mu ,\sigma ,b|W_n)=\frac{1}{p\left(W_n\right)}\frac{b^n}{{\sigma }^n}\frac{{\mathrm{e}}^{-S_n/\sigma +n\mu /\sigma }}{{\sum }_{i=1}^n{\left(1+{\mathrm{e}}^{-(W_i-\mu )/\sigma }\right)}^{b+1}}p(\mu ,\sigma ,b).$$

(19)

Note that the integrals stated in the expressions defined in (15) and (18), as well as the solution of the equation established in (16), are not simple to obtain explicit closed forms. Moreover, in this paper, we can specify the prior PDF without supposing independence between $\mu$, $\sigma$ and b, via the Jeffreys prior. Note that the Jeffreys prior does not work well for models of multidimensional parameters [34]. Nonetheless, we use it in the next section to find the Bayesian estimates in a particular case for the sake of comparison.

Next, we proceed with the MCMC estimation. We focus on those traditional tuned search procedures for proposed and associated priors. We test the method proposed in [35], available for R, Python, MatLab, and C++. Using the R language for our simulated scenarios here, it was not possible to achieve reasonable results.

Considering the MCMC algorithm, that is, ${\theta }^{(k+1)}$ with distribution defined by $q(·\mid {\theta }^{\left(k\right)})$ being accepted or rejected according to the MH algorithm, we set

$$\alpha ({\theta }^{\left(k\right)};{\theta }^{(k+1)}):=min\left\{1,\frac{p^{★}\left({\theta }^{(k+1)}\right|W_n)q({\theta }^{\left(k\right)}\mid {\theta }^{(k+1)})}{p^{★}\left({\theta }^{\left(k\right)}\right|W_n)q({\theta }^{(k+1)}\mid {\theta }^{\left(k\right)})}\right\},$$

(20)

where $p^{★}$ is defined in (17) and q is a specific conditional PDF also called candidate kernel. Considering the likelihood functions previously defined, and the general prior specification given, the posterior distribution has a complex form. Then, the Bayesian estimation can be implemented by employing MCMC methods, which make it simpler to obtain efficient sampling from the marginal posterior distributions. A crucial step in designing an effective sampling regime, implemented in the acceptance/rejection step stated in (20), is the choice of a kernel q.

5. Simulation Studies

We recall that, for large samples, the biases of the ML and Bayesian estimators are negligible. However, for small and moderate sample sizes, the second-order biases of the ML estimators may be large. We use Monte Carlo simulation to evaluate the finite sample performance of the ML and Bayes estimators for the parameter $\theta ={(\mu ,\sigma ,b)}^{\top }$ of the IGL distribution, with the PDF given in (1).

In order to analyze the point estimation results from ML and Bayesian methods, we compute, for each sample size, $n\in \{15,30,50,100\}$, each scenario specified, representing the degree skewness as: severe negative $\theta ={(0,2,0.05)}^{\top }$ –case 1–; moderate negative $\theta ={(0,4,0.01)}^{\top }$ –case 2–: zero $\theta ={(0,1,1)}^{\top }$ –case 3–; and moderate positive $\theta ={(0,6,10)}^{\top }$ –case 4–. For each of the two estimation procedures, it is reported the true bias (E$(\widehat{\theta })-\theta$); the relative bias ($(\mathrm{E}(\widehat{\theta })-\theta )/\theta$), by estimating E$(\widehat{\theta })$ with Monte Carlo experiments); and the root mean square error, that is, $\sqrt{\mathrm{MSE}}$, where MSE is the mean squared error estimated from the 5000 Monte Carlo replicates. The values of IGL distributed random variables were generated using the inversion method. More specifically, the generation of the random samples coming from $W\sim \mathrm{IGL}(\mu ,\sigma ,b)$ can be easily obtained utilizing the PDF defined in (1) through the inverse transformation, that is, employing $W=\mu -\sigma (log(U^{1/b}-1))$, where U is uniformly distributed on $(0,1)$. To find the ML estimate, the log-likelihood function defined in (8) is maximized using the Newton–Raphson algorithm with first analytic and second derivatives, both implemented in the glogisfit function of the glogis package [36].

To find the Bayesian estimate, we employ the MCMC algorithm by generating candidate ${\theta }^{(k+1)}$ with probability $\alpha ({\theta }^{\left(k\right)};{\theta }^{(k+1)})$ given in (20). The Geweke and multivariate Gelman–Rubin criteria were used to evaluate convergence. Moreover, we utilize the Ljung-Box test to detect the autocorrelations in each chain.

Next, we conduct the classical estimation approaches. In order to show numerically the difficulties they have with the ML and moment methods when estimating the vector of parameters $\theta$, even if convergence is achieved in some replicates, we implement the methods with the generated simulations.

Table 1. Percentage of updates by means of the moment estimation for the indicated values.

$\mathit{\mu }$	$\mathit{\sigma }$	b	$\mathit{n}=15$	$\mathit{n}=30$	$\mathit{n}=50$	$\mathit{n}=100$
0	2	0.05	93.5%	84.5%	83.0%	77.0%
0	4	0.1	94.5%	85.5%	84.5%	80.1%
0	1	1	98.0%	98.5%	99.5%	100.0%
0	6	10	88.5%	80.5%	82.0%	70.5%

reports the percentages of updates by means of the moment estimation. In addition, we apply the minimum chi-square estimation procedure [37] implemented in [38] for the mipfp package of R without achieving reliable results. We do not show these results in the paper due to restrictions of space.

Table 2. Estimates with the moment method for 1000 replicates of the listed size sample, case, parameter, and indicator.

$n=15$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−19.09638	13.31917	0.49371	−18.81575	13.75909	0.51013
	Bias	−19.09638	11.31917	0.44371	−18.81575	9.75909	0.41013
	Relative bias	-	5.65958	8.8742	-	2.43977	4.1013
	Standard deviation	14.25511	4.4811	0.21426	14.80999	4.53213	0.22428
	MSE	567.87987	148.20384	0.24278	573.36827	115.78	0.21851
$n=15$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.2572	0.98217	1.39323	11.42986	4.90052	2.78161
	Bias	−0.2572	−0.01783	0.39323	11.42986	-1.09948	−7.21839
	Relative bias	-	−0.01783	0.39323	-	−0.18325	−0.72184
	Standard deviation	1.0315	0.26367	1.42171	5.86088	1.39999	3.27409
	MSE	1.13014	0.06984	2.17588	164.99163	3.16883	62.82485
$n=30$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−14.43447	11.63792	0.38711	−14.29753	12.14067	0.4039
	Bias	−14.43447	9.63792	0.33711	−14.29753	8.14067	0.3039
	Relative bias	-	4.81896	6.7422	-	2.03517	3.039
	Standard deviation	9.2901	3.64442	0.14531	9.57106	3.64783	0.15087
	MSE	294.6599	106.17127	0.13476	296.0245	79.57715	0.11512
$n=30$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.1717	0.97755	1.468	11.33066	4.89198	3.03369
	Bias	−0.1717	−0.02245	0.468	11.33066	−1.10802	−6.96631
	Relative bias	-	−0.02245	0.468	-	−0.18467	−0.69663
	Standard deviation	1.04284	0.2548	3.32295	5.87308	1.14218	4.2669
	MSE	1.117	0.06543	11.26101	162.87692	2.53227	66.73593
$n=50$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−10.79721	10.10306	0.31954	−10.5922	10.65166	0.33576
	Bias	−10.79721	8.10306	0.26954	-10.5922	6.65166	0.23576
	Relative bias	-	4.05153	5.3908	-	1.66291	2.3576
	Standard deviation	8.38061	3.34231	0.13118	8.62776	3.34787	0.13503
	MSE	186.81438	76.83061	0.08986	186.63291	55.45283	0.07382
$n=50$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.12718	0.97281	1.30513	9.72513	5.09282	3.83445
	Bias	−0.12718	−0.02719	0.30513	9.72513	−0.90718	−6.16555
	Relative bias	-	−0.02719	0.30513	-	−0.1512	−0.61656
	Standard deviation	0.92368	0.21061	1.61055	5.55714	0.89339	5.36295
	MSE	0.86936	0.04509	2.68698	125.45997	1.62112	66.77523
$n=100$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−8.57027	9.05183	0.272	−8.21102	9.56039	0.28589
	Bias	−8.57027	7.05183	0.222	−8.21102	5.56039	0.18589
	Relative bias	-	3.52591	4.44	-	1.3901	1.8589
	Standard deviation	5.91348	2.60411	0.09776	6.22217	2.70276	0.10234
	MSE	108.41877	56.50969	0.05884	106.13622	38.22284	0.04503
$n=100$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.19489	1.00876	1.23021	7.13316	5.4028	5.68398
	Bias	−0.19489	0.00876	0.23021	7.13316	−0.5972	−4.31602
	Relative bias	-	0.00876	0.23021	-	−0.09953	−0.4316
	Standard deviation	0.67739	0.15986	0.68861	5.93529	0.75779	7.43532
	MSE	0.49684	0.02563	0.52718	86.10964	0.93089	73.91197

reports the results obtained for the estimates with the moment method and 1000 replicates for the size samples, cases, parameter, and indicators mentioned. These results are not satisfactory and, for the ML method, we do not report the results because they are totally unsatisfactory.

Now, the Bayesian approach is conducted. In the case where the parameter vector $\theta$ of the IGL distribution incorporates the three parameters, that is, $\mu$, $\sigma$ and b, we employ the traditional MCMC algorithm as an alternative to the classical estimation. Understanding the choice of prior distributions is undoubtedly the most controversial aspect of any Bayesian analysis [39,40], we proceed to specify two priori classes (informative and non-informative) for a parameter space of the IGL distribution.

We use the most common convergence criteria in the literature on the topic related to the Gelman–Rubin, Geweke, and Ljung-Box methods. The Gelman–Rubin statistic is built with at least two chains, where the variance is compared within and between them. If the variability is similar within the chains, they are considered to have converged. In terms of its practical use, the value of the statistic must be less that 1.2 [41]. The Geweke test is constructed with a chain, which is usually done by taking the initial 10% and 50% of the chain and making a comparison of means between both partitions, where significant differences indicate that the chain does not converge. The Ljung-Box test is applied to each chain and tries to evaluate the correlation between chain values at different lags, in particular it was evaluated at a maximum of 15 lags. If there is any significant observed correlation, then the chain is correlated and is not considered a random sample. We employ the smallest p-value of the 15 lags, if this is not significant. Then, the others correlations, for lags less than or equal to 14, are not significant either.

We consider a prior distribution given by $p\left(\theta \right)=p(\mu ,\sigma ,b)=p\left(\mu \right)p\left(\sigma \right)p\left(b\right),$ where $p\left(\mu \right)$, $p\left(\sigma \right)$ and $p\left(b\right)$ are specified as follows. For the case where we have severe negative skewness (case 1), we use $\mu \sim \mathrm{Normal}({\gamma }_1,{\delta }_1)$, $\sigma \sim \mathrm{Gamma}({\alpha }_1,{\beta }_1)$ and $b\sim \mathrm{Beta}({\nu }_1,{\psi }_1)$. When we have moderate negative skewness (case 2), we assume $\mu \sim \mathrm{Normal}({\gamma }_2,{\delta }_2)$, $\sigma \sim \mathrm{Lognormal}(\eta ,\kappa )$ and $b\sim \mathrm{Inverse}$-$\mathrm{gamma}(\rho ,\tau )$. When symmetry exists (case 3), we state $\mu \sim \mathrm{Normal}({\gamma }_3,{\delta }_3)$, $\sigma \sim \mathrm{Gamma}({\alpha }_2,{\beta }_2)$ and $b\sim \mathrm{Gamma}({\alpha }_3,{\beta }_3)$. For the case of positive skewness (case 4), we employ ${\theta }_1$, ${\theta }_2$ and ${\theta }_3$ as the dependent joint. We numerically test different priors and values of the hyperparameters to find the convergence of the chains and the sufficient number of iterations. Understanding that there are other methods to set hyperparameter values, for example, empirical Bayes and James–Stein estimators, the development and implementation of these other methods will be studied in future work.

The criterion for eliciting the kernels was mainly oriented to choosing those structures that allowed us to make the algorithm more efficient, that is, to decrease the rejection rate in the step that involves the Metropolis–Hastings ratio, as stated in (20). Consequently for the simulation study, we have that:

• Case 1: $\mathrm{E}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=0$, $\mathrm{Var}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={({\mu }^{★}+{\sigma }^{★}/2+b^{★})}^{1/2}$,
  $\mathrm{E}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\sigma }^{★}+b^{★}$, $\mathrm{Var}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+3{\sigma }^{★}$,
  $\mathrm{E}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}/5+b^{★}$, $\mathrm{Var}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\sigma }^{★}$.
• Case 2: $\mathrm{E}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=0$, $\mathrm{Var}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\left({\mu }^{★}\right)}^2+{\sigma }^{★}/4+b^{★}/10$,
  $\mathrm{E}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\sigma }^{★}$, $\mathrm{Var}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+{\sigma }^{★}+b^{★}$,
  $\mathrm{E}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=b^{★}$, $\mathrm{Var}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+{\sigma }^{★}$.
• Case 3: $\mathrm{E}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★}$, $\mathrm{Var}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={({\sigma }^{★}+b^{★})}^2$,
  $\mathrm{E}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\sigma }^{★}$, $\mathrm{Var}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\left({\mu }^{★}/2\right)}^2+b^{★}$.
• Case 4: $\mathrm{E}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=0$, $\mathrm{Var}\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+{\sigma }^{★}/4+b^{★}/10$,
  $\mathrm{E}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\left({\mu }^{★}\right)}^2+{\sigma }^{★}$, $\mathrm{Var}\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+{\sigma }^{★}/6+b^{★}/10$,
  $\mathrm{E}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=9+{\mu }^{★2}+b^{★}/10$, $\mathrm{Var}\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★2}+{\sigma }^{★}/2$,

where ★ denotes the true value of the MCMC.

Next, we consider dependent and non-informative priors. A class of advisable priori distributions are of the non-informative type, whose construction method is relatively simple, and without privileging specific behavior for any element of the parameter vector. This can be achieved via a Jeffreys prior defined as

$$p(\theta )\propto det{(\mathit{K}({\theta }^{★}))}^{1/2},$$

where $\mathit{K}({\theta }^{★})$ is defined with the elements given in (12), and ${\theta }^{★}={(\sigma ,b)}^{\top }$. We observe that the behavior of the Jeffreys prior depends only on ${\theta }^{★}={(\sigma ,b)}^{\top }$. In Figure 2, we show the behavior of $det{(\mathit{K}({\theta }^{★}))}^{1/2}$ on the subspace of ${\theta }^{★}$, say $(\sigma ,b)\in (0,7.0)\times (0,12.5)$, highlighting four region: (a) $(0,0.5)\times (0,0.5)$; (b) $(0,0.5)\times (4,6)$; (c) $(4,6)\times (0,0.5)$; and (d) $(5,7)\times (9,11)$. Thus, in general, note that the PDF of a Jeffreys prior takes extreme values when the scale parameter ($\sigma$) is close to the origin and in greater degree when the shape (b) is close to the origin (see Figure 3a,b, respectively). Whereas if we move away from it (see Figure 3c,d), the Jeffreys prior PDF (in this figure denoted as z) takes smaller values. This is the reason why it is recommended to use it in cases where empirical asymmetry is positive, that is, $\theta ={(0,6,10)}^{\top }$, in order to avoid those parts of the domain where the PDF has very strong changes and thus damages the efficiency of the algorithm. The above does not disregard the use of the independent and non-informative priors. To achieve a reasonable rejection rate, these results indicate to us that a Jeffrey priori on the region (c) should be considered.

Table 3. Bayesian estimate (posterior mean) for 1000 replicates of the listed size sample, case, parameter, and indicator.

$n=15$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.019	3.102	0.081	0.061	3.711	0.093
	Bias	−0.019	1.102	0.031	0.061	−0.289	−0.007
	Relative bias	-	0.551	0.62	-	−0.072	−0.07
	Standard deviation	0.08	0.632	0.024	0.406	0.157	0.024
	MSE	0.007	1.614	0.002	0.169	0.108	0.001
$n=15$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	0.086	0.962	1.363	0.084	5.918	10.517
	Bias	0.086	−0.038	0.363	0.084	−0.082	0.517
	Relative bias	-	−0.038	0.363	-	−0.014	0.052
	Standard deviation	0.85	0.331	0.683	0.172	0.746	0.742
	MSE	0.729	0.111	0.598	0.037	0.563	0.818
$n=30$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.025	2.737	0.07	0.062	3.705	0.093
	Bias	−0.025	0.737	0.02	0.062	−0.295	−0.007
	Relative bias	-	0.368	0.4	-	−0.074	−0.07
	Standard deviation	0.122	0.719	0.02	0.558	0.221	0.018
	MSE	0.015	1.061	0.001	0.315	0.136	0
$n=30$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.097	0.996	1.419	0.088	5.956	10.429
	Bias	−0.097	−0.004	0.419	0.088	−0.044	0.429
	Relative bias	-	−0.004	0.419	-	−0.007	0.043
	Standard deviation	0.836	0.267	0.761	0.21	0.581	0.949
	MSE	0.709	0.072	0.755	0.052	0.339	1.085
$n=50$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.007	2.502	0.064	0.086	3.71	0.093
	Bias	−0.007	0.502	0.014	0.086	−0.29	−0.007
	Relative bias	-	0.251	0.28	-	−0.072	−0.07
	Standard deviation	0.151	0.666	0.018	0.634	0.241	0.014
	MSE	0.023	0.696	0.001	0.409	0.142	0
$n=50$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.186	1.017	1.404	0.078	5.961	10.379
	Bias	−0.186	0.017	0.404	0.078	−0.039	0.379
	Relative bias	-	0.017	0.404	-	−0.006	0.038
	Standard deviation	0.754	0.196	0.754	0.238	0.47	1.063
	MSE	0.604	0.039	0.731	0.063	0.223	1.273
$n=100$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.012	2.383	0.061	0.05	3.779	0.096
	Bias	−0.012	0.383	0.011	0.05	−0.221	−0.004
	Relative bias	-	0.192	0.22	-	−0.055	−0.040
	Standard deviation	0.194	0.572	0.016	0.707	0.272	0.011
	MSE	0.038	0.474	0	0.502	0.123	0
$n=100$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.221	1.027	1.314	0.067	5.96	10.389
	Bias	−0.221	0.027	0.314	0.067	−0.04	0.389
	Relative bias	-	0.027	0.314	-	−0.007	0.039
	Standard deviation	0.561	0.151	0.531	0.187	0.357	0.879
	MSE	0.364	0.023	0.38	0.039	0.129	0.924

and

Table 4. Bayesian estimate (posterior median) for 1000 replicates of the listed size sample, case, parameter, and indicator.

$n=15$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.026	2.747	0.069	0.065	3.638	0.087
	Bias	−0.026	0.747	0.019	0.065	−0.362	−0.013
	Relative bias	-	0.374	0.38	-	−0.09	−0.13
	Standard deviation	0.095	0.624	0.023	0.398	0.158	0.022
	MSE	0.01	0.947	0.001	0.163	0.156	0.001
$n=15$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	0.22	0.907	0.935	0.081	5.84	10.267
	Bias	0.22	−0.093	−0.065	0.081	−0.16	0.267
	Relative bias	-	−0.093	−0.065	-	−0.027	0.027
	Standard deviation	0.924	0.346	0.513	0.191	0.724	0.728
	MSE	0.903	0.128	0.267	0.043	0.549	0.602
$n=30$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.034	2.472	0.062	0.064	3.632	0.089
	Bias	−0.034	0.472	0.012	0.064	−0.368	−0.011
	Relative bias	-	0.236	0.24	-	−0.092	−0.11
	Standard deviation	0.141	0.717	0.02	0.543	0.22	0.017
	MSE	0.021	0.737	0.001	0.299	0.184	0
$n=30$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	0.017	0.966	1.079	0.092	5.903	10.214
	Bias	0.017	−0.034	0.079	0.092	−0.097	0.214
	Relative bias	-	−0.034	0.079	-	−0.016	0.021
	Standard deviation	0.845	0.27	0.571	0.245	0.568	0.935
	MSE	0.714	0.074	0.332	0.068	0.332	0.92
$n=50$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.014	2.304	0.058	0.09	3.644	0.09
	Bias	−0.014	0.304	0.008	0.09	−0.356	−0.01
	Relative bias	-	0.152	0.16	-	−0.089	−0.1
	Standard deviation	0.18	0.672	0.018	0.624	0.239	0.014
	MSE	0.033	0.545	0	0.398	0.184	0
$n=50$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.089	0.996	1.142	0.085	5.922	10.185
	Bias	−0.089	−0.004	0.142	0.085	−0.078	0.185
	Relative bias	-	−0.004	0.142	-	−0.013	0.018
	Standard deviation	0.751	0.197	0.57	0.282	0.46	1.045
	MSE	0.572	0.039	0.346	0.087	0.217	1.126
$n=100$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Mean	−0.022	2.258	0.057	0.064	3.72	0.093
	Bias	−0.022	0.258	0.007	0.064	−0.28	−0.007
	Relative bias	-	0.129	0.14	-	−0.07	−0.07
	Standard deviation	0.231	0.582	0.016	0.693	0.272	0.011
	MSE	0.054	0.406	0	0.485	0.152	0
$n=100$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Mean	−0.154	1.015	1.16	0.076	5.938	10.218
	Bias	−0.154	0.015	0.16	0.076	−0.062	0.218
	Relative bias	-	0.015	0.16	-	−0.01	0.022
	Standard deviation	0.554	0.152	0.43	0.223	0.349	0.856
	MSE	0.331	0.023	0.211	0.055	0.125	0.781

report the results obtained for the posterior mean and median estimates, respectively, with the Bayesian method and 1000 replicates for the size samples, cases, parameter, and indicators mentioned. These results are quite satisfactory and helpful. We use the convergence Gelman–Rubin, Geweke, and Ljung-Box criteria, with their results being reported in

Table 5. Results of the diagnostics for 1000 replicates of the listed size sample, case, parameter, and test.

$n=15$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Gelman-Rubin	-	1.2081	-	-	1.0644	-
	Geweke	0.9516	0.648	0.6535	0.6399	0.3411	0.1022
	Ljung-Box	0.112	0.2358	0.2497	0.0729	0.2563	0.1289
$n=15$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Gelman-Rubin	-	1.1755	-	-	1.0706	-
	Geweke	0.2091	0.438	0.3774	0.7777	0.5109	0.7443
	Ljung-Box	0.525	0.4431	0.2327	0.6502	0.1333	0.1462
$n=30$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Gelman-Rubin	-	1.017	-	-	1.0218	-
	Geweke	0.2203	0.0943	0.1702	0.9785	0.9788	0.427
	Ljung-Box	0.6654	0.0271	0.0022	0.737	0.4986	0.7545
$n=30$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Gelman-Rubin	-	1.0758	-	-	1.0414	-
	Geweke	0.5641	0.2097	0.4419	0.722	0.5419	0.7722
	Ljung-Box	0.2461	0.8968	0.193	0.0779	0.0902	0.3097
$n=50$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Gelman-Rubin	-	1.0209	-	-	1.0305	-
	Geweke	0.8796	0.542	0.6853	0.5699	0.2915	0.5914
	Ljung-Box	0.0822	0.4869	0.2854	0.2587	0.7124	0.6139
$n=50$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Gelman-Rubin	-	1.1348	-	-	1.1341	-
	Geweke	0.5635	0.4084	0.5817	0.9604	0.7169	0.6362
	Ljung-Box	0.1675	0.6055	0.0738	0.4626	0.2712	0.0961
$n=100$	Cases 1–2	$\mu =0$	$\sigma =2$	$b=0.05$	$\mu =0$	$\sigma =4$	$b=0.1$
	Gelman-Rubin	-	1.1659	-	-	1.0252	-
	Geweke	0.6261	0.5756	0.8681	0.8935	0.6845	0.9877
	Ljung-Box	0.0754	0.1719	0.0878	0.0942	0.2024	0.5346
$n=100$	Cases 3–4	$\mu =0$	$\sigma =1$	$b=1$	$\mu =0$	$\sigma =6$	$b=10$
	Gelman-Rubin	-	1.0038	-	-	1.1379	-
	Geweke	0.4396	0.1397	0.4124	0.0997	0.8876	0.5265
	Ljung-Box	0.7793	0.1469	0.3915	0.1579	0.2634	0.1191

. From these results, we observe that all criteria employed here tell us that all chains converge, except, although debatable, in the case where the distribution has a severe negative skewness and with a small sample size. In a future study, we hope to provide a detailed discussion on the selection of priori distributions. We recall that the main objective of our study is to show that the Bayesian estimation is an alternative to the estimation problem of the three-parameter IGL distribution.

6. Empirical Application

Next, we present an application of real-world data corresponding to the solvent extraction (SX) process of copper. This process considers several stages of extraction and re-extraction with their respective operating variables. These variables, as in any process, often generate difficulties to the operation as the lack of efficiency on the copper recovery itself. Generally, difficulties occur in the following stage—the electrowinning (EW). Thus, it results from the dragging of unwanted impurities into the electrolyte. One of the purposes in the metallurgical area is to analyze the process in the stages of solvent extraction and electro-obtaining, identifying the most relevant variables, according to their: (i) impact on the process; (ii) empirical and design operational conditions; and (ii) parameters of control. The pregnant leach solution (or impure aqueous solution, known as PLS in short), from the leaching process, feeds the SX area which, by an organic solvent, is transferred to a pure and concentrated copper sulfate solution, called rich electrolyte, which is sent to the electro-obtaining area. The total PLS flow that feeds the plant is one of the main factors that determine the efficiency of copper extraction [42].

The PLS daily flow throughout the extraction process in the area of SX sulfides (W) was registered between years 2015 and 2019, with a sample size of $n=728$. A descriptive summary, a histogram and the empirical cumulate distribution of the data set are shown in

Table 6. Descriptive summary of PLS daily flow data.

Mean	Median	Variance	Standard Deviation	Coefficient of Skewness
10,739.57	10,999.74	318,485.03	564.34	−3.50

and Figure 4, respectively. We can observe from Table 6 and Figure 4 that the data distribution has a moderate positive skewness. This is a very complex situation to estimate the parameters of the IGL distribution by the moment and ML methods, as it was shown in Section 2 and simulation study of Section 4. We use the Bayesian estimation method.

To apply Bayesian estimation, we consider independent and informative prior distributions for the parameter $\theta ={(\mu ,\sigma ,b)}^{\top }$ and $p(\mu ,\sigma ,b)=p\left(\mu \right)\times p\left(\sigma \right)\times p\left(b\right)$, where $\mu \sim \mathrm{Normal}(\gamma ,\delta ),$ $\sigma \sim \mathrm{Gamma}(\rho ,\tau )$ and $b\sim \mathrm{Gamma}(\rho ,\tau ).$ The hyperparameters are set by utilizing the criteria mentioned. For the kernels, we have:

E$\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★}$, Var$\left({\mu }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={({\sigma }^{★}/100+b^{★}/100)}^2$,

E$\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\sigma }^{★}$, Var$\left({\sigma }_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★}/4000+b^{★}/120$,

E$\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})=b^{★}$, Var$\left(b_t\right|{\mu }^{★},{\sigma }^{★},b^{★})={\mu }^{★}/\mathrm{10,000}+{\sigma }^{★}/1000$.

In the diagnostic analysis, we have that the multivariate Gelman–Rubin factor is 1.01. The p-values from the Geweke test are in

Table 7. p-value of the Geweke test for the indicated parameter with PLS daily flow data.

Parameter	$\mathit{\mu }$	$\mathit{\sigma }$	b
p-value	0.60	0.34	0.04

, from where we can observe that, only for the chains associated with the shape parameter b, this was significant at 5%. Moreover, in Figure 4, the histogram with the fitted IGL PDF (a) and the empirical cumulate distribution with the fitted IGL cumulative distribution function (b) are shown, from where we can see that a good fit was achieved in comparison to the other estimation methods considered in this article. The Bayesian estimate of $\theta ={(\mu ,\sigma ,b)}^{\top }$ is reported in

Table 8. Classical and Bayesian estimates of $\theta ={(\mu ,\sigma ,b)}^{\top }$ with PLS daily flow data.

	Parameter
Indicator	$\mathit{\mu }$	$\mathit{\sigma }$	b
Estimate	604.75	1184.63	3866.72
Posterior mean	11,000.27	177.55	0.49
Posterior median	11,000.24	177.55	0.49
Variance	1.18	4.15	0.00
Standar deviation	1.09	2.04	0.02

, and the autocorrelation functions (ACFs) of the chains are shown in Figure 5. Thus, by having the daily flow distribution of the PLS characterized via the estimated values of the parameters, and the solvent extraction process, the entire LX-SX-EW circuit is better controlled since the behavior of the PLS flow is better characterized. In this way, the control becomes more efficient.

7. Conclusions

In relation to the classical methods used for estimating the parameters of the type I generalized logistic distribution, that is, the maximum likelihood and moment estimators, we can conclude from a numerical point of view that: (a) in the case of the ML estimator, as expected, there are convergence problems, which can be seen clearly for the distribution with negative asymmetry; (b) there is no pattern in the ranges of the parameters or in the sample sizes where the estimators have good behavior; and (c) for moment estimation, the convergence problems are accentuated, whose evidences come to complement a more analytical point of view, as described in [10,43].

The proposed analysis with a real-world data set based on the type I generalized logistic distribution by using the Bayesian approach notoriously corrects the problems detected with the maximum likelihood and moment estimators. Another aspect in favor of our proposal is to be able to find congruent estimates independent of the range of values that parameters and sample size could have. In general, the Bayesian approach turned out to be robust in the estimation of parameters because the average of the estimates resembles the values of the parameters with which the different samples were simulated and with quite small variance. Thus, this work showed the impossibility of finding consistent solutions to a problem that was reflected in the classical methods of estimation. Therefore, in the light foregoing, the future research will be analytically addressing some solutions to improve or prevent the non-existence of the maximum likelihood estimator of the third parameter of the type I generalized logistic distribution. Exploration of different computational algorithms in the estimation of parameters for a class of distributions could be carried out in the line of the work proposed in [44].

As mentioned, some other commonly used classical estimation approaches are related to ordinary/weighted least squares [20], percentiles, minimum distances, and bootstrap methods [21]. The comparison of additional estimation methods to the considered ones in the present study is beyond our objectives. However, we are planning such a comparison for future work. Once we solve the problem of estimation in the type I generalized logistic distribution, different aspects related to sampling, inference, and modeling can be studied. Then, another important challenge to be implemented is related to calculating the sample size when estimating the mean or other parameters of the type I generalized logistic distribution [45], as well as formulating more complex modeling structures, such as regression, temporal, spatial, functional, partial least squares, and errors-in-variables settings [46,47,48,49,50]. We hope to report findings related to these topics in new investigations.