# Bayesian and Classical Inference for the Generalized Log-Logistic Distribution with Applications to Survival Data

Authors: Abdisalam Hassan Muse; Samuel Mwalili; Oscar Ngesa; Saad J. Almalki; Gamal A. Abd-Elmougod
Journal: Computational Intelligence and Neuroscience
 (2021)

---

## Abstract

The generalized log-logistic distribution is especially useful for modelling survival data with variable hazard rate shapes because it extends the log-logistic distribution by adding an extra parameter to the classical distribution, resulting in greater flexibility in analyzing and modelling various data types. We derive the fundamental mathematical and statistical properties of the proposed distribution in this paper. Many well-known lifetime special submodels are included in the proposed distribution, including the Weibull, log-logistic, exponential, and Burr XII distributions. The maximum likelihood method was used to estimate the unknown parameters of the proposed distribution, and a Monte Carlo simulation study was run to assess the estimators’ performance. This distribution is significant because it can model both monotone and nonmonotone hazard rate functions, which are quite common in survival and reliability data analysis. Furthermore, the proposed distribution’s flexibility and usefulness are demonstrated in a real-world data set and compared to its submodels, the Weibull, log-logistic, and Burr XII distributions, as well as other three-parameter parametric survival distributions, such as the exponentiated Weibull distribution, the three-parameter log-normal distribution, the three-parameter (or the shifted) log-logistic distribution, the three-parameter gamma distribution, and an exponentiated Weibull distribution. The proposed distribution is plausible, according to the goodness-of-fit, log-likelihood, and information criterion values. Finally, for the data set, Bayesian inference and Gibb’s sampling performance are used to compute the approximate Bayes estimates as well as the highest posterior density credible intervals, and the convergence diagnostic techniques based on Markov chain Monte Carlo techniques were used.


---

## Body

## 1. Introduction

Applied statisticians use many probability distributions for reliability and survival studies. The distributions could be applied in different fields such as medicine, engineering, economy, industrial and physical fields, and so many other fields. Exponential distributions, generalized exponential distributions, gamma distributions, generalized gamma distributions, extreme value distributions, Weibull distributions, log-logistic distributions, log-normal distributions, Burr XII distributions, and generalized Weibull distributions are among the most frequently used distributions in survival and reliability analysis.Typically, researchers in reliability and survival analysis are concerned with the development of new probability models. Log-logistic (LL) distribution is one of the parametric distributions that can be used as a life-testing distribution because of the simplicity of its cumulative distribution and survival function which can both be stated in closed form and because it belongs to the Scale-Shape family [1]. LL is one of the right-skewed, heavy-tailed functions that can be used as an alternative to a log-normal distribution. It resembles the log-normal distribution in shape but has heavier tails. Log-logistic distribution is particularly applicable to model nonmonotone (i.e., unimodal) hazard functions.It is well understood that the log-logistic model is not appropriate for modelling where the failure rate is monotonic when analyzing time-to-event data with parametric models. It is suitable to use an extension of the model which has a monotone hazard function. Departures from the monotonicity of distribution are typically studied in terms of its shape or more specifically in terms of its skewness (also referred to as asymmetry) and kurtosis.In this study, we focus on a modification of the log-logistic model because it resembles the log-normal distribution in shape but is better suited for the application in the analysis of survival data when dealing with incomplete data, such as censored observations which are common in such data [2]. The presence of incomplete observations causes difficulties when using log-normal or inverse Gaussian models, since the survival functions in these cases are complicated. On the other hand, since the logarithms of small positive numbers are large negative numbers, the log-normal distribution may give undue weight to very short survival times [1]. For the reasons stated above, we will focus on the log-logistic model whose hazard rate exhibits the aforementioned behaviour.However, due to the log-logistic model’s symmetric property, it may be inadequate for cases where the hazard rate is heavily tailed or skewed, as well as for modelling censored survival data [3–5]. In this study, we studied a modification (or generalization) of the log-logistic parametric survival model and referring to this as the generalized log-logistic distribution given in [6]. The generalized log-logistic distribution reflects the structure of the heavy tails and the skewness and it significantly outperformed the log-logistic distribution in general.In the statistical literature, with the aim of increasing the versatility of the log-logistic distribution in modelling survival time data, different generalized forms of the distribution have recently been proposed, including a new extension of the LL distribution with applications to actuarial data sets [7], alpha power transformed LL distribution [8, 9], transmuted four-parameter generalized LL distribution [10, 11], a new three-parameter LL distribution [12], extended log-logistic distribution [13], exponentiated LL geometric distribution [14], the LL Weibull distribution [15], beta LL distribution [16], McDonald LL distribution [2], transmuted LL distribution [17], Marshal-Olkin LL distribution [18], the Zografos-Balakrishnan LL distribution [19], and exponentiated LL distribution [20]. More details about the modifications and recent generalizations of the log-logistic distribution can be found in [21].In addition, other authors have studied the Bayesian inference of the LL distribution and some of its generalizations. dos Santos et al. [22] developed a Bayesian analysis of the transmuted LL distribution. Yahaya and Dewu [23] studied the Bayesian estimation of the scale parameter for the LL distribution using Chi-square and Maxwell priors. Abbas and Tang [24] studied the objective Bayesian analysis of the LL distribution using the reference and Jeffreys prior. Al-Shomrani et al. [25] focused on the application of the Markov chain Monte Carlo (McMC) techniques for estimating the unknown parameters of the LL distribution. Guure et al. [26] explored the Bayesian inference of the LL distribution for the interval-censored data. Kang et al. [27] proposed the noninformative priors for the LL distribution. Chaudhary and Kumar [28] studied the Bayesian estimation of the three-parameter exponentiated LL distribution. Akhtar et al. [29] discussed the Bayesian analysis of the LL distribution using the Laplace approximation. Chaudhary [30] proposed the Bayesian analysis of the two-parameter exponentiated LL distribution.The log-logistic distribution has large-scale applications in analyzing time-to-event data. The model is closed under both proportionality (multiplication) of failure time and proportionality of odds, though it is not a proportional hazard (PH) model. However, regarding this issue, Khan and Khosa [6] presented generalized log-logistic distribution that belongs to the proportional hazard models. The proposed distribution has similar properties to the 2-parameter log-logistic distribution and approaches the Weibull distribution in limit. However, its statistical and mathematical properties, as well as inferential procedures, have not received attention so far. On the other hand, they discussed the classical inference of the proposed distribution under the PH regression framework. However, much work still has to be done. In this paper, we focused on the Bayesian and classical inference of the generalized log-logistic distribution as a generalized distribution, not as a regression model.Additionally, for the applied cases, especially in the survival modelling, the GLL model could be applicable in the following cases: (1) modelling the “asymmetric monotonically right-skewed” heavy tail data sets; (2) modelling the “bathtub-shaped hazard rate” data sets like data set I; (3) in “survival analysis,” the GLL distribution could be chosen for modelling proportional hazard frameworks; (4) in the medical field, the GLL distribution could be considered in modelling the “bladder cancer data sets” which have “reversed bathtub-shaped HRF” as illustrated in data set I; and (5) in the reliability and survival analysis, the proposed distribution can be an alternative to the Weibull distribution since it can be closed under both accelerated failure time (AFT) and PH models since the Weibull distribution fails to model unimodal data. For these based on ground reasons, we are motivated to study and introduce the GLL distribution.Thus, the main goal of this research article is to propose and study a generalized log-logistic distribution, which extends the exponential, Weibull, log-logistic, and Burr XII distributions, with the hope that the proposed distribution may have a better fit compared to these distributions and other 3-parametric distributions in certain practical situations. In addition, we would provide a comprehensive account of the mathematical and statistical properties of the proposed model. The proposed model’s formulae are simple and tractable, and, with the use of modern computer software and its numerical capabilities, the proposed model could be a great addition to the arsenal of applied mathematicians and statisticians in the areas like medicine, engineering, economics, social sciences, and biology, among others. Finally, we discussed the Bayesian model formulation of the proposed distribution.The rest of the paper is organized as follows. Section2 describes the distribution functions for the GLL distribution, its submodel distributions, and some of its basic properties. Some mathematical properties of the GLL distribution are derived in Section 3. Section 4 describes the maximum likelihood for the estimation parameters of GLL distribution. Section 5 discusses the findings of a simulation study that was conducted to estimate and compare the performances of the proposed estimators. Section 6 presents an analysis of a real-life data set. The Bayesian model formulation for the proposed distribution is discussed in Section 7. Section 8 presents the Bayesian analysis of a real-life data set using Markov chain Monte Carlo techniques. Finally, Section 9 summarizes the study with some concluding remarks.


## 2. The Generalized Log-Logistic Distribution

The generalized log-logistic distribution is a continuous probability distribution with positive supportℝ on a subset of 0,∞ with three parameters. It is a generalization of the two-parameter log-logistic distribution. The generalization of log-logistic distribution for censored survival data can be traced back to Singh et al. [3] who discussed a generalized log-logistic distribution and applied it to censored survival data and proposed a generalized log-logistic model and introduced the shape parameter and then they used it to fit a lung cancer data. Prentice [31] proposed a generalization for quantile response data and discussed several of its uses.Since many continuous probability distributions are commonly applied for parametric models in survival analysis like the exponential, Gompertz, Weibull, log-normal, log-logistic, and the gamma distribution, GLL is also applicable for survival data analysis. There are a number of probability functions that are related to continuous probability distributions; we will concentrate on functions that are related to the lifetime distributions as a random variable in this study.


### 2.1. Hazard (Failure) Rate Function

SThe hazard (failure) rate function plays an important role in survival analysis. It is the most popular function for analyzing and modelling lifetime data because of its intuitive interpretation of the amount of risk to fail associated with a unit timet, applicable for describing the lifetime distribution of engineered and other components. The hazard rate is more informative than all of the other functions in lifetime distributions. Because of this, the authors in [6] started their work by defining the hazard rate of the GLL distribution. Cox and Oakes [32] described the reason why the hazard rate is considered when we are dealing with the survival data. They gave a number of reasons including the fact that hazard rate-based models are often convenient when there is incomplete information (censoring) or there are several types of failure rates; also hazard rate is a special form of the intensity function, and last but not least the hazard rate function can be derived from all other functions that we use to describe lifetime distributions.The hazard rate function describes how the instantaneous failure rate changes over time. For the GLL distribution, the hazard rate function plots are given in(1)hx;θ=αkkxα−11+ηxα,x≥0,k,α,η>0,where k>0,β>0,η>0 are parameters and θ=k,α,η′.It can be easily seen from equation (1) that the hazard rate function is monotonically decreasing for α≤1 and unimodal when α≤1. That is, it initially increases to a maximum at t=α−1/λα1/α and then decreases to zero monotonically as t⟶∞. The HRF plots are shown in Figure 1.Figure 1
The hazard curve of the GLL distribution.


### 2.2. Submodels

The proposed distribution consists of a number of important submodels that are widely used in parametric survival modelling. These include the log-logistic distribution, the standard log-logistic distribution, the Burr XII distribution, the Weibull distribution, and the exponential distribution. The propositions below relate the GLL to the log-logistic, standard log-logistic, Burr XII, Weibull, and exponential distributions.


#### 2.2.1. Log-Logistic Distribution

Proposition 1.
LetX∼GLLα,k,η. If ηdepends onkviak=η, then the hazard rate function of (1) reduces to the hazard rate function of the log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(2)hx;θ=αkkxα−11+ηxα.
If we replaceη=k, it gives us(3)hx;θ=αkkxα−11+kxα=αkkxα−11+kxα,which is the hazard rate function form of a log-logistic distribution with the two unknown parameters k,α. When θ=k,α′, k=1/βis the rate parameter.
It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=1/kα−11/α).


#### 2.2.2. Standard Log-Logistic Distribution

Proposition 2.
LetX∼GLLα,k,η. If ηdepends onkviak=η=1, then the hazard rate function of (1) reduces to the hazard rate function of the standard log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(4)hx;θ=αkkxα−11+ηxα.
If we replaceη=k=1, it gives us(5)hx;θ=α⋅11⋅xα−11+1⋅xα=αxα−11+xα,which is the hazard rate function form of a standard log-logistic distribution with one unknown parameter α. Hence, the proof.
It should be noted thatx>0, is the distribution‘s support, and α is the distribution's shape parameter. It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for 0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=α−11/α).


#### 2.2.3. Burr XII Distribution

Proposition 3.
LetX∼GLLα,k,η. If ηdepends onkviaη=kλ−1/α,λ>0, then the hazard rate function of (1) reduces to the hazard rate function of the Burr XII distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(6)hx;θ=αkkxα−11+ηxα.
If we replaceη=kλ−1/α, it gives us(7)hx;θ=αkkxα−11+kλ−1/αxα=αkkxα−11+kλ−α/αxα=αkxα−11+xα,which is the hazard rate function form of a Burr XII distribution with two unknown parameters α,k. Hence, the proof.
The Burr XII hazard function is monotonically decreasing forα≤1 and upside-down bathtub shapes curve for α>1 (which means that it initially increases, attains a maximum at x=α−11/α, and then decreases to zero at (x⟶∞).


#### 2.2.4. Weibull Distribution

Proposition 4.
LetX∼GLLα,k,η. If ηα⟶0, then the hazard rate function of the GLL (1) approaches the hazard rate function of the Weibull distribution.Proof.
If we now letηα⟶0, then, from the hazard rate function of the GLL given by(8)hx;θ=αkkxα−11+ηxα,we have that(9)hx;θ=αkkxα−11+0,which by simplifying gives(10)ht;θ=αkkxα−1,which is a hazard function of a Weibull distribution with the unknown parameters α,k. This property of the GLL enables it to handle monotonically increasing hazard satisfactorily with α>1 and λ close to zero (very small).
It is clear from (10) that, for 0<α<1, the hazard rate function decreases, for α>1, the hazard rate function increases, and for α=1, the hazard rate function decreases.
The distribution reduces to exponential forα=1.


#### 2.2.5. Exponential Distribution

Proposition 5.
Similarly, if we now letα=1, then the hazard rate function of (10) reduces to the hazard rate function of the exponential distribution.Proof.
From (10), we have that the hazard rate function is(11)ht;θ=αkkxα−1,and if we replace α=1,(12)ht;θ=k⋅11⋅t1−1,which by simplifying gives(13)ht;θ=k,which is the hazard rate function of an exponential distribution. This property makes the exponential distribution be inadequate to describe survival data. Hence, the proof.
The summary of the submodels for the proposed distribution is summarized in Table1.Table 1
Summary of submodels from the GLL distribution.
DistributionsαηkLog-logistic distributionαη=kk=ηWeibull distributionηα⟶0ηα⟶0kExponential distributionα=1η⟶0kStandard log-logistic distributionαη=k=1k=η=1Burr XII distributionαη=kλ−1/α,λ>0η=kλ−1/α,λ>0


### 2.3. The Probability Density Function

The pdf of the GLL distribution with three unknown parameters can be obtained by applying the following equation and the pdf plots are shown in Figure2.(14)fx;θ=hx,θexp−∫0xhxdx.Figure 2
The pdf curve of the GLL distribution.Simplifying gives(15)fx;θ=αkkxα−11+ηxαkα/ηα+1,x≥0,k,α,η>0.


### 2.4. The Survival (or Reliability) Function

The survival (reliability) function of the GLL distribution that represents the probability that observation does not fail untilt is given below and its plots are shown in Figure 3.(16)Sx;θ=fx;θhx;θ.Figure 3
The survival curves of the GLL distribution.Simplifying gives(17)Sx;θ=1+ηxα−kα/ηα,x≥0,k,α,η>0.


### 2.5. Cumulative Distribution Function of the GLL Distribution

The cumulative distribution function (CDF), also known as the lifetime distribution function, of the GLL distribution is of the form below and the CDF plots are shown in Figure4.(18)Fx;θ=1+ηxαkα/ηα−11+ηxαkα/ηα,x≥0,k,α,η>0,where k>0,β>0,η>0 are parameters and θ=k,α,η′.Figure 4
The CDF curve of the GLL distribution.


### 2.6. The Reversed Hazard Rate Function

The reversed hazard rate (also known as the retro hazard) is defined as the ratio of pdf to the corresponding CDF. The retro hazard is written as follows:(19)rx;θ=fx;θFx;θ.Reversed hazard rate function plays an important role in the analysis of censored data and in the estimation of the survival function. The following equation gives us the basic relationship between hazard rate function and the reversed hazard rate function.(20)rx;θ=hx;θSx;θ1−Sx;θ.The applications of hazard rate function in survival analysis are well known. Recently, the reversed hazard rate function has gained popularity among applied statisticians; for more information, see [33, 34]. Block et al. [33] showed that the hazard rate function plays an essential role in the analysis of right-censored data, while the retro hazard plays an essential role in the analysis of left-censored data.The reversed hazard rate function of the GLL distribution takes the form(21)rx;θ=fx;θFx;θ=αkktα−1/1+λxαkα/λα+11+λxαkα/λα.Simplifying gives(22)rx;θ=αkkxα−11+λxαkα/λα+1−1+λxα,x≥0,k,α,η>0.The reversed hazard rate plots are shown in Figure5.Figure 5
The reversed hazard rate curves of the GLL distribution.


### 2.7. The Cumulative Hazard Function

The cumulative hazard function of the GLL distribution takes the form(23)Hx;θ=−logSx;θ=∫0xhx;θdx.Simplifying gives(24)Hx;θ=kαλαlog1+λxα,x≥0,k,α,η>0,where k>0,α>0,λ>0 are parameters and θ=k,α,η′.


### 2.8. The Hazard Rate Average (FRA) Function

The HRA function ofX is expressed as(25)HRAx;θ=Hx;θx=∫0xhx;θdxx,x>0,where Hx;θ is the cumulative hazard function. An analysis of HRAx;θ on t enables us to find increasing hazard rate average and decreasing hazard rate average.


## 2.1. Hazard (Failure) Rate Function

SThe hazard (failure) rate function plays an important role in survival analysis. It is the most popular function for analyzing and modelling lifetime data because of its intuitive interpretation of the amount of risk to fail associated with a unit timet, applicable for describing the lifetime distribution of engineered and other components. The hazard rate is more informative than all of the other functions in lifetime distributions. Because of this, the authors in [6] started their work by defining the hazard rate of the GLL distribution. Cox and Oakes [32] described the reason why the hazard rate is considered when we are dealing with the survival data. They gave a number of reasons including the fact that hazard rate-based models are often convenient when there is incomplete information (censoring) or there are several types of failure rates; also hazard rate is a special form of the intensity function, and last but not least the hazard rate function can be derived from all other functions that we use to describe lifetime distributions.The hazard rate function describes how the instantaneous failure rate changes over time. For the GLL distribution, the hazard rate function plots are given in(1)hx;θ=αkkxα−11+ηxα,x≥0,k,α,η>0,where k>0,β>0,η>0 are parameters and θ=k,α,η′.It can be easily seen from equation (1) that the hazard rate function is monotonically decreasing for α≤1 and unimodal when α≤1. That is, it initially increases to a maximum at t=α−1/λα1/α and then decreases to zero monotonically as t⟶∞. The HRF plots are shown in Figure 1.Figure 1
The hazard curve of the GLL distribution.


## 2.2. Submodels

The proposed distribution consists of a number of important submodels that are widely used in parametric survival modelling. These include the log-logistic distribution, the standard log-logistic distribution, the Burr XII distribution, the Weibull distribution, and the exponential distribution. The propositions below relate the GLL to the log-logistic, standard log-logistic, Burr XII, Weibull, and exponential distributions.


### 2.2.1. Log-Logistic Distribution

Proposition 1.
LetX∼GLLα,k,η. If ηdepends onkviak=η, then the hazard rate function of (1) reduces to the hazard rate function of the log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(2)hx;θ=αkkxα−11+ηxα.
If we replaceη=k, it gives us(3)hx;θ=αkkxα−11+kxα=αkkxα−11+kxα,which is the hazard rate function form of a log-logistic distribution with the two unknown parameters k,α. When θ=k,α′, k=1/βis the rate parameter.
It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=1/kα−11/α).


### 2.2.2. Standard Log-Logistic Distribution

Proposition 2.
LetX∼GLLα,k,η. If ηdepends onkviak=η=1, then the hazard rate function of (1) reduces to the hazard rate function of the standard log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(4)hx;θ=αkkxα−11+ηxα.
If we replaceη=k=1, it gives us(5)hx;θ=α⋅11⋅xα−11+1⋅xα=αxα−11+xα,which is the hazard rate function form of a standard log-logistic distribution with one unknown parameter α. Hence, the proof.
It should be noted thatx>0, is the distribution‘s support, and α is the distribution's shape parameter. It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for 0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=α−11/α).


### 2.2.3. Burr XII Distribution

Proposition 3.
LetX∼GLLα,k,η. If ηdepends onkviaη=kλ−1/α,λ>0, then the hazard rate function of (1) reduces to the hazard rate function of the Burr XII distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(6)hx;θ=αkkxα−11+ηxα.
If we replaceη=kλ−1/α, it gives us(7)hx;θ=αkkxα−11+kλ−1/αxα=αkkxα−11+kλ−α/αxα=αkxα−11+xα,which is the hazard rate function form of a Burr XII distribution with two unknown parameters α,k. Hence, the proof.
The Burr XII hazard function is monotonically decreasing forα≤1 and upside-down bathtub shapes curve for α>1 (which means that it initially increases, attains a maximum at x=α−11/α, and then decreases to zero at (x⟶∞).


### 2.2.4. Weibull Distribution

Proposition 4.
LetX∼GLLα,k,η. If ηα⟶0, then the hazard rate function of the GLL (1) approaches the hazard rate function of the Weibull distribution.Proof.
If we now letηα⟶0, then, from the hazard rate function of the GLL given by(8)hx;θ=αkkxα−11+ηxα,we have that(9)hx;θ=αkkxα−11+0,which by simplifying gives(10)ht;θ=αkkxα−1,which is a hazard function of a Weibull distribution with the unknown parameters α,k. This property of the GLL enables it to handle monotonically increasing hazard satisfactorily with α>1 and λ close to zero (very small).
It is clear from (10) that, for 0<α<1, the hazard rate function decreases, for α>1, the hazard rate function increases, and for α=1, the hazard rate function decreases.
The distribution reduces to exponential forα=1.


### 2.2.5. Exponential Distribution

Proposition 5.
Similarly, if we now letα=1, then the hazard rate function of (10) reduces to the hazard rate function of the exponential distribution.Proof.
From (10), we have that the hazard rate function is(11)ht;θ=αkkxα−1,and if we replace α=1,(12)ht;θ=k⋅11⋅t1−1,which by simplifying gives(13)ht;θ=k,which is the hazard rate function of an exponential distribution. This property makes the exponential distribution be inadequate to describe survival data. Hence, the proof.
The summary of the submodels for the proposed distribution is summarized in Table1.Table 1
Summary of submodels from the GLL distribution.
DistributionsαηkLog-logistic distributionαη=kk=ηWeibull distributionηα⟶0ηα⟶0kExponential distributionα=1η⟶0kStandard log-logistic distributionαη=k=1k=η=1Burr XII distributionαη=kλ−1/α,λ>0η=kλ−1/α,λ>0


## 2.2.1. Log-Logistic Distribution

Proposition 1.
LetX∼GLLα,k,η. If ηdepends onkviak=η, then the hazard rate function of (1) reduces to the hazard rate function of the log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(2)hx;θ=αkkxα−11+ηxα.
If we replaceη=k, it gives us(3)hx;θ=αkkxα−11+kxα=αkkxα−11+kxα,which is the hazard rate function form of a log-logistic distribution with the two unknown parameters k,α. When θ=k,α′, k=1/βis the rate parameter.
It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=1/kα−11/α).


## 2.2.2. Standard Log-Logistic Distribution

Proposition 2.
LetX∼GLLα,k,η. If ηdepends onkviak=η=1, then the hazard rate function of (1) reduces to the hazard rate function of the standard log-logistic distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(4)hx;θ=αkkxα−11+ηxα.
If we replaceη=k=1, it gives us(5)hx;θ=α⋅11⋅xα−11+1⋅xα=αxα−11+xα,which is the hazard rate function form of a standard log-logistic distribution with one unknown parameter α. Hence, the proof.
It should be noted thatx>0, is the distribution‘s support, and α is the distribution's shape parameter. It is easy to verify that the hazard rate function of the log-logistic distribution is monotonically decreasing for 0<α≤1 and unimodal for α>1 (decreases and then increases with the maximum at x=α−11/α).


## 2.2.3. Burr XII Distribution

Proposition 3.
LetX∼GLLα,k,η. If ηdepends onkviaη=kλ−1/α,λ>0, then the hazard rate function of (1) reduces to the hazard rate function of the Burr XII distribution.Proof.
The hazard rate function of the generalized log-logistic distribution is given by(6)hx;θ=αkkxα−11+ηxα.
If we replaceη=kλ−1/α, it gives us(7)hx;θ=αkkxα−11+kλ−1/αxα=αkkxα−11+kλ−α/αxα=αkxα−11+xα,which is the hazard rate function form of a Burr XII distribution with two unknown parameters α,k. Hence, the proof.
The Burr XII hazard function is monotonically decreasing forα≤1 and upside-down bathtub shapes curve for α>1 (which means that it initially increases, attains a maximum at x=α−11/α, and then decreases to zero at (x⟶∞).


## 2.2.4. Weibull Distribution

Proposition 4.
LetX∼GLLα,k,η. If ηα⟶0, then the hazard rate function of the GLL (1) approaches the hazard rate function of the Weibull distribution.Proof.
If we now letηα⟶0, then, from the hazard rate function of the GLL given by(8)hx;θ=αkkxα−11+ηxα,we have that(9)hx;θ=αkkxα−11+0,which by simplifying gives(10)ht;θ=αkkxα−1,which is a hazard function of a Weibull distribution with the unknown parameters α,k. This property of the GLL enables it to handle monotonically increasing hazard satisfactorily with α>1 and λ close to zero (very small).
It is clear from (10) that, for 0<α<1, the hazard rate function decreases, for α>1, the hazard rate function increases, and for α=1, the hazard rate function decreases.
The distribution reduces to exponential forα=1.


## 2.2.5. Exponential Distribution

Proposition 5.
Similarly, if we now letα=1, then the hazard rate function of (10) reduces to the hazard rate function of the exponential distribution.Proof.
From (10), we have that the hazard rate function is(11)ht;θ=αkkxα−1,and if we replace α=1,(12)ht;θ=k⋅11⋅t1−1,which by simplifying gives(13)ht;θ=k,which is the hazard rate function of an exponential distribution. This property makes the exponential distribution be inadequate to describe survival data. Hence, the proof.
The summary of the submodels for the proposed distribution is summarized in Table1.Table 1
Summary of submodels from the GLL distribution.
DistributionsαηkLog-logistic distributionαη=kk=ηWeibull distributionηα⟶0ηα⟶0kExponential distributionα=1η⟶0kStandard log-logistic distributionαη=k=1k=η=1Burr XII distributionαη=kλ−1/α,λ>0η=kλ−1/α,λ>0


## 2.3. The Probability Density Function

The pdf of the GLL distribution with three unknown parameters can be obtained by applying the following equation and the pdf plots are shown in Figure2.(14)fx;θ=hx,θexp−∫0xhxdx.Figure 2
The pdf curve of the GLL distribution.Simplifying gives(15)fx;θ=αkkxα−11+ηxαkα/ηα+1,x≥0,k,α,η>0.


## 2.4. The Survival (or Reliability) Function

The survival (reliability) function of the GLL distribution that represents the probability that observation does not fail untilt is given below and its plots are shown in Figure 3.(16)Sx;θ=fx;θhx;θ.Figure 3
The survival curves of the GLL distribution.Simplifying gives(17)Sx;θ=1+ηxα−kα/ηα,x≥0,k,α,η>0.


## 2.5. Cumulative Distribution Function of the GLL Distribution

The cumulative distribution function (CDF), also known as the lifetime distribution function, of the GLL distribution is of the form below and the CDF plots are shown in Figure4.(18)Fx;θ=1+ηxαkα/ηα−11+ηxαkα/ηα,x≥0,k,α,η>0,where k>0,β>0,η>0 are parameters and θ=k,α,η′.Figure 4
The CDF curve of the GLL distribution.


## 2.6. The Reversed Hazard Rate Function

The reversed hazard rate (also known as the retro hazard) is defined as the ratio of pdf to the corresponding CDF. The retro hazard is written as follows:(19)rx;θ=fx;θFx;θ.Reversed hazard rate function plays an important role in the analysis of censored data and in the estimation of the survival function. The following equation gives us the basic relationship between hazard rate function and the reversed hazard rate function.(20)rx;θ=hx;θSx;θ1−Sx;θ.The applications of hazard rate function in survival analysis are well known. Recently, the reversed hazard rate function has gained popularity among applied statisticians; for more information, see [33, 34]. Block et al. [33] showed that the hazard rate function plays an essential role in the analysis of right-censored data, while the retro hazard plays an essential role in the analysis of left-censored data.The reversed hazard rate function of the GLL distribution takes the form(21)rx;θ=fx;θFx;θ=αkktα−1/1+λxαkα/λα+11+λxαkα/λα.Simplifying gives(22)rx;θ=αkkxα−11+λxαkα/λα+1−1+λxα,x≥0,k,α,η>0.The reversed hazard rate plots are shown in Figure5.Figure 5
The reversed hazard rate curves of the GLL distribution.


## 2.7. The Cumulative Hazard Function

The cumulative hazard function of the GLL distribution takes the form(23)Hx;θ=−logSx;θ=∫0xhx;θdx.Simplifying gives(24)Hx;θ=kαλαlog1+λxα,x≥0,k,α,η>0,where k>0,α>0,λ>0 are parameters and θ=k,α,η′.


## 2.8. The Hazard Rate Average (FRA) Function

The HRA function ofX is expressed as(25)HRAx;θ=Hx;θx=∫0xhx;θdxx,x>0,where Hx;θ is the cumulative hazard function. An analysis of HRAx;θ on t enables us to find increasing hazard rate average and decreasing hazard rate average.


## 3. Some Mathematical Properties of the GLL Distribution

In this section, we present some mathematical properties of the GLL distribution. The functions that we discussed in Section2 are not the only ways that we can define the GLL distribution, but there are other mathematical functions that we can use to describe the lifetime distributions of a random variable X. These include quantile function and its related results, moments and its related properties, rth central moments, residual life and reversed residual life functions, and other mathematical properties.


### 3.1. The Quantile Function and Related Results

The quantile function (which is the inverse of the CDF) is crucial in statistical and quantitative data analysis. A probability distribution can be defined in terms of either the quantile function or the cumulative distribution function [35]. The quantiles of the proposed distribution with various parameter values are given in Table 2.Table 2
Quantiles of the proposed distribution for different parameter values.
Quantilesk,α,η(0.5, 0.5, 0.5)(5.0, 1.5, 1.5)(4.0, 4.0, 2.5)(3.0, 2.0, 3.0)(5.0, 3.0, 2.0)0.10.02470.04490.14270.11110.09450.20.12500.07450.17250.16670.12160.30.36730.10260.19450.21820.14240.40.88890.13140.21340.27220.16080.51.99990.16270.23120.33330.17830.64.49990.19850.24890.40820.19610.710.88890.24210.26810.50920.21550.832.00000.30060.29060.66670.23850.9162.00000.39720.32221.00000.2707Theorem 1.
IfT∼GLLk,α,η, then the quantile function, lower quartile, median, and the upper quartile of the GLL distribution, respectively, are given by(26)Xq=F−1q;k,α,η=1/1−pηα/kα−11/αη,(27)Xq1=4/3ηα/kα−11/αη,(28)Xq2=Median=2ηα/kα−11/αη,(29)Xq3=4ηα/kα−11/αη.Proof.
The quantile function of GLL distribution is derived by finding the value ofQ for which(30)1−1+ηxα−kα/ηα=p,Xq=F−1q;k,α,η=1+ηqα−kα/ηα=1−p=11+ηqαkα/ηα=1−p=1+ηqαkα/ηα=11−p=1+ηqα=11−pηα/kα=ηqα=11−pηα/kα−1=ηq=11−pηα/kα−11/α,∴q=1/1−pηα/kα−11/αη,where p∈0,1. k>0,α>0,η>0. Hence the proof.
Similarly, we can prove (27)–(29) by applying the following values: the lower quartile = 1/4, median = 2/4 = 1/2, and the upper quartile = 3/4.
Lower quartile is(31)Xq1=4/3ηα/kα−11/αη.
Median is(32)Xq2=median=2ηα/kα−11/αη.
Upper quartile is(33)Xq3=4ηα/kα−11/αη.


#### 3.1.1. Skewness and Kurtosis

The following relationship defines the mathematical form of the Galton Skewness and Moors Kurtosis of the GLL model with three parameters:(34)SK=Q3/4+Q1/4−2Q2/4Q3/4−Q1/4,KM=Q7/8+Q3/8−Q5/8−Q1/8Q6/8−Q2/8,where Q describes different quartile values.The above equations can be determined as functions of the GLL quantile function. The advantages of these measures are that they are less sensitive in the presence of outliers and that they exist even when the distribution is lacking moments.


### 3.2. The Random Deviate Generation Functions

Let U be a random variable with a uniform distribution0,1 and an inverse CDF, F.. Then any sample drawn from F−1u is assumed to have been drawn from F.. As a result, using GLL k,α,η, the random deviate can be generated as follows:(35)x=1/1−u−1λα/kα1/αλ,0<u<1,where u follows U0,1 distribution.


### 3.3. Therth Moments and Related Results

Numerous important characteristics and properties of a probability distribution such as mean, variance, kurtosis, and skewness can be obtained from its moments. Moments are extremely important and play a central role in statistical analysis, especially in applications. The important moment functions, such as the moments,rth moment, rth central moment, mean, variance, skewness, and kurtosis of the proposed distribution, are presented.Theorem 2.
IfT∼GLLk,α,η, then the rth power, negative moments, and logarithmic moments are given, respectively, by(36)ETr=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1,forαkαηα>r,(37)ET−r=λα+rkαΓkα/ηα+1Γkα/ηα−r/αΓr/α+1.Proof.
We have(38)ETr=∫0∞trft;k,α,ηdt=∫0∞trαkktα−11+ηtαkα/ηβ+1dt=αkΓkα/ηα+1∫0∞trktα−11+ηtαdt=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1,forαkαηα>r.
Similarly, we can prove (37).


#### 3.3.1. Mean and Variance

Corollary 1.
IfT∼GLLk,α,η, then the mean and variance are given, respectively, as follows.The mean of the GLL distribution is(39)μ=ET=kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+1.This is provided thatαkα/ηα>1.The Variance of the GLL distribution is(40)σ2=VT=ET2−ET2=kαηα+2Γkα/ηα−2/αΓ2/α+1Γkα/ηα+1−kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+12.This is provided thatαkα/ηα>2.


### 3.4. Therth Central Moments

Corollary 2.
IfT∼GLLk,α,η, then the cumulants of the first, second, and rth central moments, are given, respectively, by(41)c1=μ1′=ET=kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+1,c2=μ2′−μ12′=ET2−ET2=kαηα+2Γkα/ηα−2/αΓ2/α+1Γkα/ηα+1−kβηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+12,cr=μr′−∑n=1r−1r−1n−1cnμr−m′=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1−∑n=1r−1r−1n−1cnkαηα+r−nΓkα/ηα−r−n/αΓr−n/α+1Γkα/ηα+1.Hence, from Corollary 2, we can derive the skewness and kurtosis of the GLL distribution by computing, respectively:(42)Skewness=c3σ23/2,Kurtosis=c4σ22.


### 3.5. Residual and Reverse Residual Life

The residual life has broader applications in survival analysis and risk management. The residual lifetime of the GLL random variable is calculated as follows:(43)Rtx=Sx+tSt,Rtx=1+ηx+tα−kα/ηα1+ηtα−kα/ηα.In addition, the reverse residual life of the generalized log-logistic random variable can be calculated as follows:(44)R^tx=Sx−tSt,R^tx=1+ηx−tα−kα/ηα1+ηtα−kα/ηα.From Table3, the GLL distribution is clearly numerically versatile in terms of means and variance. Furthermore, the values of CS show that it can be right-skewed, nearly symmetrical, or slightly left-skewed. The CK values show that the GLL distribution can be mesokurtic, leptokurtic, or platykurtic. All of these characteristics demonstrate the GLL distribution flexibility, which remains appealing for modelling purposes.Table 3
1st five moments, standard deviation, skewness, and kurtosis of the GLL distribution for some parameter values.
Momentsk,α,η(0.5, 0.5, 0.5)(1.0, 1.5, 1.5)(1.5, 2.0, 2.5)(2.0, 5.0, 3.0)(1.0, 1.0, 2.0)(4.0, 4.5, 0.2)(5.0, 4.0, 0.5)μ1′0.10340.20650.24320.27950.15470.22810.1813μ2′0.05670.12920.14820.17410.08930.05540.0354μ3′0.03880.09250.10360.12040.06190.01410.0073μ4′0.02940.07150.07870.09000.04710.00370.0016μ5′0.02370.05810.06310.07110.03800.00100.0004SD0.21460.29430.29840.30980.25570.05750.0509CV2.07431.42501.22701.10811.65290.25210.2805CS2.37431.17840.91090.61001.6648−0.1784−0.0871CK7.88423.03182.52402.02384.66282.80812.7479The mean and variance plots for different values of alpha and kappa parameters are shown in Figure6, while the skewness and kurtosis plots are shown in Figure 7.Figure 6
The mean and variance plot for several combinations of alpha and kappa parameters.Figure 7
The skewness and kurtosis plot for several combinations of alpha and kappa.


## 3.1. The Quantile Function and Related Results

The quantile function (which is the inverse of the CDF) is crucial in statistical and quantitative data analysis. A probability distribution can be defined in terms of either the quantile function or the cumulative distribution function [35]. The quantiles of the proposed distribution with various parameter values are given in Table 2.Table 2
Quantiles of the proposed distribution for different parameter values.
Quantilesk,α,η(0.5, 0.5, 0.5)(5.0, 1.5, 1.5)(4.0, 4.0, 2.5)(3.0, 2.0, 3.0)(5.0, 3.0, 2.0)0.10.02470.04490.14270.11110.09450.20.12500.07450.17250.16670.12160.30.36730.10260.19450.21820.14240.40.88890.13140.21340.27220.16080.51.99990.16270.23120.33330.17830.64.49990.19850.24890.40820.19610.710.88890.24210.26810.50920.21550.832.00000.30060.29060.66670.23850.9162.00000.39720.32221.00000.2707Theorem 1.
IfT∼GLLk,α,η, then the quantile function, lower quartile, median, and the upper quartile of the GLL distribution, respectively, are given by(26)Xq=F−1q;k,α,η=1/1−pηα/kα−11/αη,(27)Xq1=4/3ηα/kα−11/αη,(28)Xq2=Median=2ηα/kα−11/αη,(29)Xq3=4ηα/kα−11/αη.Proof.
The quantile function of GLL distribution is derived by finding the value ofQ for which(30)1−1+ηxα−kα/ηα=p,Xq=F−1q;k,α,η=1+ηqα−kα/ηα=1−p=11+ηqαkα/ηα=1−p=1+ηqαkα/ηα=11−p=1+ηqα=11−pηα/kα=ηqα=11−pηα/kα−1=ηq=11−pηα/kα−11/α,∴q=1/1−pηα/kα−11/αη,where p∈0,1. k>0,α>0,η>0. Hence the proof.
Similarly, we can prove (27)–(29) by applying the following values: the lower quartile = 1/4, median = 2/4 = 1/2, and the upper quartile = 3/4.
Lower quartile is(31)Xq1=4/3ηα/kα−11/αη.
Median is(32)Xq2=median=2ηα/kα−11/αη.
Upper quartile is(33)Xq3=4ηα/kα−11/αη.


### 3.1.1. Skewness and Kurtosis

The following relationship defines the mathematical form of the Galton Skewness and Moors Kurtosis of the GLL model with three parameters:(34)SK=Q3/4+Q1/4−2Q2/4Q3/4−Q1/4,KM=Q7/8+Q3/8−Q5/8−Q1/8Q6/8−Q2/8,where Q describes different quartile values.The above equations can be determined as functions of the GLL quantile function. The advantages of these measures are that they are less sensitive in the presence of outliers and that they exist even when the distribution is lacking moments.


## 3.1.1. Skewness and Kurtosis

The following relationship defines the mathematical form of the Galton Skewness and Moors Kurtosis of the GLL model with three parameters:(34)SK=Q3/4+Q1/4−2Q2/4Q3/4−Q1/4,KM=Q7/8+Q3/8−Q5/8−Q1/8Q6/8−Q2/8,where Q describes different quartile values.The above equations can be determined as functions of the GLL quantile function. The advantages of these measures are that they are less sensitive in the presence of outliers and that they exist even when the distribution is lacking moments.


## 3.2. The Random Deviate Generation Functions

Let U be a random variable with a uniform distribution0,1 and an inverse CDF, F.. Then any sample drawn from F−1u is assumed to have been drawn from F.. As a result, using GLL k,α,η, the random deviate can be generated as follows:(35)x=1/1−u−1λα/kα1/αλ,0<u<1,where u follows U0,1 distribution.


## 3.3. Therth Moments and Related Results

Numerous important characteristics and properties of a probability distribution such as mean, variance, kurtosis, and skewness can be obtained from its moments. Moments are extremely important and play a central role in statistical analysis, especially in applications. The important moment functions, such as the moments,rth moment, rth central moment, mean, variance, skewness, and kurtosis of the proposed distribution, are presented.Theorem 2.
IfT∼GLLk,α,η, then the rth power, negative moments, and logarithmic moments are given, respectively, by(36)ETr=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1,forαkαηα>r,(37)ET−r=λα+rkαΓkα/ηα+1Γkα/ηα−r/αΓr/α+1.Proof.
We have(38)ETr=∫0∞trft;k,α,ηdt=∫0∞trαkktα−11+ηtαkα/ηβ+1dt=αkΓkα/ηα+1∫0∞trktα−11+ηtαdt=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1,forαkαηα>r.
Similarly, we can prove (37).


### 3.3.1. Mean and Variance

Corollary 1.
IfT∼GLLk,α,η, then the mean and variance are given, respectively, as follows.The mean of the GLL distribution is(39)μ=ET=kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+1.This is provided thatαkα/ηα>1.The Variance of the GLL distribution is(40)σ2=VT=ET2−ET2=kαηα+2Γkα/ηα−2/αΓ2/α+1Γkα/ηα+1−kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+12.This is provided thatαkα/ηα>2.


## 3.3.1. Mean and Variance

Corollary 1.
IfT∼GLLk,α,η, then the mean and variance are given, respectively, as follows.The mean of the GLL distribution is(39)μ=ET=kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+1.This is provided thatαkα/ηα>1.The Variance of the GLL distribution is(40)σ2=VT=ET2−ET2=kαηα+2Γkα/ηα−2/αΓ2/α+1Γkα/ηα+1−kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+12.This is provided thatαkα/ηα>2.


## 3.4. Therth Central Moments

Corollary 2.
IfT∼GLLk,α,η, then the cumulants of the first, second, and rth central moments, are given, respectively, by(41)c1=μ1′=ET=kαηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+1,c2=μ2′−μ12′=ET2−ET2=kαηα+2Γkα/ηα−2/αΓ2/α+1Γkα/ηα+1−kβηαΓkα/ηα−1/αΓ1/α+1Γkα/ηα+12,cr=μr′−∑n=1r−1r−1n−1cnμr−m′=kαηα+rΓkα/ηα−r/αΓr/α+1Γkα/ηα+1−∑n=1r−1r−1n−1cnkαηα+r−nΓkα/ηα−r−n/αΓr−n/α+1Γkα/ηα+1.Hence, from Corollary 2, we can derive the skewness and kurtosis of the GLL distribution by computing, respectively:(42)Skewness=c3σ23/2,Kurtosis=c4σ22.


## 3.5. Residual and Reverse Residual Life

The residual life has broader applications in survival analysis and risk management. The residual lifetime of the GLL random variable is calculated as follows:(43)Rtx=Sx+tSt,Rtx=1+ηx+tα−kα/ηα1+ηtα−kα/ηα.In addition, the reverse residual life of the generalized log-logistic random variable can be calculated as follows:(44)R^tx=Sx−tSt,R^tx=1+ηx−tα−kα/ηα1+ηtα−kα/ηα.From Table3, the GLL distribution is clearly numerically versatile in terms of means and variance. Furthermore, the values of CS show that it can be right-skewed, nearly symmetrical, or slightly left-skewed. The CK values show that the GLL distribution can be mesokurtic, leptokurtic, or platykurtic. All of these characteristics demonstrate the GLL distribution flexibility, which remains appealing for modelling purposes.Table 3
1st five moments, standard deviation, skewness, and kurtosis of the GLL distribution for some parameter values.
Momentsk,α,η(0.5, 0.5, 0.5)(1.0, 1.5, 1.5)(1.5, 2.0, 2.5)(2.0, 5.0, 3.0)(1.0, 1.0, 2.0)(4.0, 4.5, 0.2)(5.0, 4.0, 0.5)μ1′0.10340.20650.24320.27950.15470.22810.1813μ2′0.05670.12920.14820.17410.08930.05540.0354μ3′0.03880.09250.10360.12040.06190.01410.0073μ4′0.02940.07150.07870.09000.04710.00370.0016μ5′0.02370.05810.06310.07110.03800.00100.0004SD0.21460.29430.29840.30980.25570.05750.0509CV2.07431.42501.22701.10811.65290.25210.2805CS2.37431.17840.91090.61001.6648−0.1784−0.0871CK7.88423.03182.52402.02384.66282.80812.7479The mean and variance plots for different values of alpha and kappa parameters are shown in Figure6, while the skewness and kurtosis plots are shown in Figure 7.Figure 6
The mean and variance plot for several combinations of alpha and kappa parameters.Figure 7
The skewness and kurtosis plot for several combinations of alpha and kappa.


## 4. Maximum Likelihood Estimation (MLE)

In this section, the unknown parameters of the generalized log-logistic distribution based on a complete sample are estimated using the maximum likelihood method. LetX1,X2,…,Xn indicate a random sample of the complete GLL data, and then the sample’s likelihood function is given as(45)L=∏i=1nfxi,α,k,η,Lx;α,k,η=∏i=1nαkkxiα−11+ηxiαkα/ηα+1.The log-likelihood function may be expressed as(46)ℓ=nlogαk+α−1∑i=1nlogkxi−∑i=1nlog1+ηxiα−kη∑i=1nlog1+ηxiα.By taking the first derivatives of the log-likelihood function in equation (48) with respect to α,k,andη and fixing the outcome to zero, we have(47)∂ℓ∂α=nα+∑i=1nlogkxi−∑i=1nηxiαlogηxi1+ηxiα−kη∑i=1nηxiαlogηxi1+ηxiα,(48)∂ℓ∂k=nk+nkα−1−1η∑i=1nlog1+ηxi,(49)∂ℓ∂η=−∑i=1nηxiαlogηxi1+ηxiα−kη2∑i=1nηxiαlogηxi1+ηxiα.It is worth noting that the MLEsα^,k^andη^ of α,k,andη, respectively, can be obtained by equating the results to zero and numerically solving the system of nonlinear equations. Because the expected information matrix is complicated, the observed information matrix Jθ is used to construct confidence intervals for the model parameters. The observed information matrix is given by(50)Jθ=−∂2ℓ∂2α∂2ℓ∂α∂k∂2ℓ∂α∂η∂2ℓ∂2k∂2ℓ∂k∂∂∂2ℓ∂2η,where θ=α,k,η′. When the usual regularity conditions are met and the parameters are within the parameter space’s interior but not on the boundary, n≅θ−θ converges in distribution to N30,I−1θ, where Iθ is the expected information matrix. When Iθ is replaced by the observed information matrix evaluated at Jθ, the asymptotic behaviour remains valid. The asymptotic multivariate normal distribution N30,J−1θ can be used to generate 1001−τ% two-sided confidence intervals for the model parameters, where τ is the significant level.


## 5. Monte Carlo Simulation Study

In this section, we assess the performance of the MLEs estimators for a finite sample of sizen using a Monte Carlo simulation study. The simulation study based on the generalized log-logistic distribution is carried out to examine the average biases (ABs), the mean square errors (MSEs), the root mean square errors (RMSEs), and maximum likelihood estimates (MLEs) for the model parameters α,k, and η. The simulation experiment was carried out using a variety of simulations with varying sample sizes and parameter values. To generate random samples for the GLL, the quantile function is given in equation [26]. The simulation study was repeated 1500 times, each with sample sizes n=50,100,…,1500, and the following parameter scenarios in set I: α=0.9,k=0.5andη=2.5, and the following parameter scenarios in set II: α=0.8,k=0.4andη=2.0.The MLEs of the GLL model are determined via the nlminb () R-function with the argument method = “BFGS”; see supplementary materials (availablehere). For each piece of simulated data, say, (α^,k^,η^) for i=1,2,…,1000, the AB, RMSE, and MSE of the parameters were computed by(51)AB=1N∑i=1Nθ^−θ,MSE=1N∑i=1Nθ^−θ2,RMSE=1N∑i=1Nθ^−θ2,where θ=α,kandη.The MLE, AB, and RMSE values of the parametersα,kandη are displayed from various sample sizes. Based on these findings, we conclude that the MLEs perform quite well in estimating the model parameters and that the estimates are fairly stable and are nearer to the true values for these sample sizes. Table 4 and Figures 8–11 show that as the sample size increases, the MSE and RMSE decrease as expected. Furthermore, as the sample size increases, the AB decreases. In addition, the MLEs of the parameters of the model are very close to the true value. As a result, the maximum likelihood estimates and their asymptotic results can be applied to construct confidence intervals for the model parameters even for a small sample size.Table 4
Monte Carlo simulation results for the GLL distribution: MLE, AB, MSEs, and RMSEs.
ParametersnIIIMLEABRMSEMLEABRMSEα502.3201.4205.2732.3301.5307.1491001.2810.3812.3861.0970.2972.2633000.9950.0951.3690.9370.1371.5126000.9210.0210.52070.8400.0400.6199000.9080.0080.0670.8040.0040.05812000.9050.0050.0600.8030.0030.04915000.9040.0040.0540.8040.0040.045k501.2460.7462.3061.2170.8172.8021000.7920.2921.2350.6130.2131.0933000.5710.0710.6650.4630.0630.4676000.5110.0110.1350.4220.0220.1999000.5080.0080.0950.4050.0050.08112000.5070.0070.0820.4040.0040.06615000.5050.0050.0730.4050.0050.063η503.2800.7803.4043.0331.0333.9441002.7800.2801.8002.2810.2811.6143002.6120.1120.9042.0560.0560.8066002.5540.0540.5882.0460.0460.5439002.5420.0420.5002.0190.0190.44212002.5260.0260.4092.0140.0140.36015002.5200.0200.3702.0300.0300.340Figure 8
Plots for MLEs and biases of the GLL model for set I of the table.Figure 9
Plots for MSEs and RMSEs of the GLL distribution for the values of set I in the table.Figure 10
Plots for MLEs and biases of the GLL distribution for the values of set II in the table.Figure 11
Plots for MSEs and RMSEs of the GLL distribution for the values of set II in the table.


## 6. Data Analysis

In this section, the proposed distribution is fully applied to real-world data set which is taken from literature to demonstrate the ability of the new model. We compare the proposed distribution with the other three parametric survival distributions including gamma, log-normal, log-logistic, exponentiated Weibull, and the Weibull distribution. Also, we have compared the GLL distribution with some of its submodels with two-parameter distribution, namely, Weibull, log-logistic, and the Burr XII distributions.The density functions of the fitted models are as follows.(1)
Weibull distribution:(52)ft=αkktα−1exp−ktα.(2)
Log-logistic distribution:(53)ft=αkktα−11+ktα2.(3)
Burr XII distribution:(54)ft=αktα−11+tα−k−1.(4)
Exponentiated Weibull distribution:(55)ft=αkλktα−11−exp−ktαλ−1exp−ktα.(5)
Three-parameter log-logistic distribution (or shifted log-logistic distribution):(56)ft=α/βt−μ/βα−11+t−μ/kβα2.(6)
Three-parameter log-normal distribution:(57)ft=αβt−μβα−1exp−t−μkβα.(7)
Three-parameter Weibull distribution:(58)ft=exp−1/2logt−μ−α/β22πβx−μ.(8)
Three-parameter Gamma distribution:(59)ft=t−μα−1exp−t−μ/ββαΓα,wheret>μ.Certain analytical measures are taken into account in order to determine which distribution best fits the applied data. These analytical measures include four discrimination measures: AIC (Akaike Information Criterion), CAIC (Consistent Akaike Information Criterion), BIC (Bayesian Information Criterion), and HQIC (Hannan-Quin Information Criterion). In addition, there are two goodness-of-fit tests: Anderson–Darling (A∗) and Cramer-von Mises (W∗).The AIC is(60)AIC=2k−2l.The BIC is(61)BIC=klnn−2l.The CAIC is(62)CAIC=2nkn−k−1−2l.The HQIC is(63)HQIC=2klnlnn−2l,where l represents the log-likelihood function evaluated as the MLEs, n denotes the sample size, and k denotes the number of model parameters. The goodness-of-fit measures under consideration are as follows.The Anderson–Darling (A∗) test statistic is given by(64)A∗=−n−1n∑i=1n2l−1×lnGXi+ln1−GXn−i+1.The Cramer-von Mises (W∗) test statistic is given by(65)W∗=112n+∑i=1n2i−12n+GXi2,where xi is the ith observation in the sample and n is the sample size; xi is calculated when the data is sorted in ascending order.The best model is the one with the lowest AIC, BIC, CAIC, and HQIC, as well as theA∗, W∗, and K-S tests. Moreover, the best model is also chosen as the one having the highest value of the log-likelihood function, and p values for the K-S statistics are also used to compare the competitive models.


### 6.1. Likelihood Ratio Test for Submodels

The GLL distribution has five submodels, namely, log-logistic distribution, Weibull distribution, Burr XII distribution, exponential distribution, and the standard log-logistic distribution. Hence, we have employed the likelihood ratio criterion to test the following hypotheses:(1)
H0:ηα⟶0; that is, the sample is from Weibull distribution. H1:ηα⟶not0; that is, the sample is GLL(2)
H0:η=k; that is, the sample is from log-logistic distribution. H1:η≠k; that is, the sample is GLL(3)
H0:kλ−1/α,λ>0; that is, the sample is from Burr XII distribution. H1:kλ−1/α,λ≤0; that is, the sample is GLL(4)
H0:η=k=1; that is, the sample is from the standard log-logistic distribution. H1:η≠1,k≠1; that is, the sample is GLL(5)
H0:η=0&α=1; that is, the sample is from an exponential distribution. H1:η≠0&α≠1; that is, the sample is GLLThe likelihood ratio test (LRT) is given by(66)LR=−2lnLθ^∗;xLθ^;x,where θ^∗ represents the restricted Maximum likelihood estimates under the null hypothesis H0 and θ^ represents the unrestricted Maximum likelihood estimates under the alternative hypothesis H1. Under the null hypothesis, the LRT follows Chi-square distribution with degrees of freedom (df) dfalt−dfnull. If the p value is less than 0.05, the null hypothesis is rejected.


### 6.2. An Application to Bladder Cancer Data Set

The following real-world data set is used to demonstrate the proposed methodology. The data in Table5 below show the remission times (in months) of a sample of 128 bladder cancer patients. The data set is available in [36]. The descriptive statistics for the data set are shown in Table 6 and the likelihood ratio test statistics for the data set are given in Table 7.Table 5
The remission times (in months) of a sample of 128 bladder cancer patients.
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, 0.08, 2.09, 3.48, 4.87, 6.94, 8.66, 13.11, 23.63, 0.20, 2.23, 3.5, 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, 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.00, 3.36, 6.93, 8.65, 12.63, 22.69.Table 6
Descriptive statistics of data set I.
MeanMedianModeVarianceSkewnessKurtosisMinimumMaximum9.3656.3955110.4353.28615.4810.0879.05Table 7
Likelihood ratio test statistic for data set I.
DistributionHypothesisLRTp valuesW2H0:ηα⟶0vsH1:H0is false8.6760.003LL2H0:ηα=kvsH1:H0is false10.8190.001Burr XIIH0:kλ−1/α,λ>0vsH1:H0is false87.472<0.001ExH0:η=0&α=1vsH1:H0is false9.1820.010Standard LLH0:η=k=1vsH1:H0is false190.150<0.001For data set I, the asymptotic variance-covariance matrix for the estimated GLL parameters is given by(67)J−1=3.0929×10−41.7255×10−35.8513×10−41.7255×10−33.1612×10−25.9347×10−35.8513×10−45.9347×10−31.5958×10−3.The information criterion values in Table8 and the goodness-of-fit tests in Table 9 both demonstrate the superiority of the proposed model over the other competing models.Table 8
Information criterion for data set I.
DistributionAICBICCAICHQICGLL825.564834.120825.756829.040LN3826.723835.279826.916830.199LL2826.937835.641827.033829.254ExpW827.393835.949827.586830.869LL3827.458836.014827.651830.934G3831.955840.511832.148835.431W2832.163837.868832.259834.481W3832.665841.221832.858836.141Burr XII910.959916.663911.055913.276Table 9
MLE estimators of the model parameters, the log-likelihood, and goodness-of-fit statistics for data set I.
DistributionsEstimates (SEs)ℓW∗A∗K−S (p value)GLL (α, kη)α = 1.410 (0.174)−409.780.0190.1280.034(0.999)k = 0.134 (0.017)η = 0.077 (0.038)ExpW (α, kλ)α = 0.275 (0.146)−410.700.0450.2910.044(0.967)k = 0.676 (0.136)λ = 2.636 (1.161)LL3 (α,β,γ)α = 0.535 (0.061)−410.730.0190.1350.038(0.993)β = 1.863 (0.106)μ = −0.293 (0.358)LN3 (α,β,γ)α = 0.877 (0.090)−410.360.0170.1150.029(0.998)β = 1.925 (0.111)μ = −0.623 (0.372)G3 (α,β,γ)α = 1.098 (0.134)−412.980.1250.7780.067(0.618)β = 8.424 (1.238)μ = 0.075 (0.018)W3 (α,β,γ)α = 1.031 (0.072)−413.330.1340.8390.080(0.387)β = 9.743 (0.908)μ = 0.077 (0.013)W2 (α,β)α = 1.049 (0.068)−414.080.1310.7840.071(0.545)k = 9.576 (0.854)BXII (α,β)α = 2.342 (0.356)−453.480.7524.5640.251(<0.005)k = 0.233 (0.040)LL2 (α,β)α = 0.578 (0.043)−411.470.0430.3100.041(0.984)k = 1.805 (0.088)The estimated pdf and CDF of the proposed distribution corresponding to the real-world data set are shown in Figure12 and the Kaplan–Meier and PP plots for the proposed distribution are shown in Figure 13.Figure 12
Estimated pdf and CDF of the GLL distribution corresponding to data set I.Figure 13
PP and Kaplan–Meier plots of the GLL distribution corresponding to data set I.


#### 6.2.1. TTT Plot

The total time test (TTT) plot plays a central role in determining the best model to fit the given data in terms of the hazard rates. This plot depicts the various forms of the hazard rate. A straight line on the TTT plot indicates that the given data has a constant hazard rate. If the plot is convex, the hazard rate will be decreased; if it is concave, the hazard rates will be increased. The plot for the bathtub shape is first convex and then concave. Similarly, if the hazard rate has an inverted bathtub shape, it will increase first (or concave) and then decrease (or convex). The TTT plot is calculated by using the following formula:(68)Grn=∑i=1rxi:n+n−rxi:n∑i=1rxi:n,r=xi:n=1,2,…,n,where xi:n are the order statistics.The TTT and box plots of the data set are presented in Figure14. These plots indicate that the empirical hazard rate function of the 1st data set is bathtub shape, monotonically increasing.Figure 14
TTT and box plots of data set I.The estimated fitted pdfs and CDFs of data set I for the competitive models are shown in Figure15.Figure 15
Some estimated fitted densities and cumulative functions of data set I.
(a)(b)


## 6.1. Likelihood Ratio Test for Submodels

The GLL distribution has five submodels, namely, log-logistic distribution, Weibull distribution, Burr XII distribution, exponential distribution, and the standard log-logistic distribution. Hence, we have employed the likelihood ratio criterion to test the following hypotheses:(1)
H0:ηα⟶0; that is, the sample is from Weibull distribution. H1:ηα⟶not0; that is, the sample is GLL(2)
H0:η=k; that is, the sample is from log-logistic distribution. H1:η≠k; that is, the sample is GLL(3)
H0:kλ−1/α,λ>0; that is, the sample is from Burr XII distribution. H1:kλ−1/α,λ≤0; that is, the sample is GLL(4)
H0:η=k=1; that is, the sample is from the standard log-logistic distribution. H1:η≠1,k≠1; that is, the sample is GLL(5)
H0:η=0&α=1; that is, the sample is from an exponential distribution. H1:η≠0&α≠1; that is, the sample is GLLThe likelihood ratio test (LRT) is given by(66)LR=−2lnLθ^∗;xLθ^;x,where θ^∗ represents the restricted Maximum likelihood estimates under the null hypothesis H0 and θ^ represents the unrestricted Maximum likelihood estimates under the alternative hypothesis H1. Under the null hypothesis, the LRT follows Chi-square distribution with degrees of freedom (df) dfalt−dfnull. If the p value is less than 0.05, the null hypothesis is rejected.


## 6.2. An Application to Bladder Cancer Data Set

The following real-world data set is used to demonstrate the proposed methodology. The data in Table5 below show the remission times (in months) of a sample of 128 bladder cancer patients. The data set is available in [36]. The descriptive statistics for the data set are shown in Table 6 and the likelihood ratio test statistics for the data set are given in Table 7.Table 5
The remission times (in months) of a sample of 128 bladder cancer patients.
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, 0.08, 2.09, 3.48, 4.87, 6.94, 8.66, 13.11, 23.63, 0.20, 2.23, 3.5, 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, 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.00, 3.36, 6.93, 8.65, 12.63, 22.69.Table 6
Descriptive statistics of data set I.
MeanMedianModeVarianceSkewnessKurtosisMinimumMaximum9.3656.3955110.4353.28615.4810.0879.05Table 7
Likelihood ratio test statistic for data set I.
DistributionHypothesisLRTp valuesW2H0:ηα⟶0vsH1:H0is false8.6760.003LL2H0:ηα=kvsH1:H0is false10.8190.001Burr XIIH0:kλ−1/α,λ>0vsH1:H0is false87.472<0.001ExH0:η=0&α=1vsH1:H0is false9.1820.010Standard LLH0:η=k=1vsH1:H0is false190.150<0.001For data set I, the asymptotic variance-covariance matrix for the estimated GLL parameters is given by(67)J−1=3.0929×10−41.7255×10−35.8513×10−41.7255×10−33.1612×10−25.9347×10−35.8513×10−45.9347×10−31.5958×10−3.The information criterion values in Table8 and the goodness-of-fit tests in Table 9 both demonstrate the superiority of the proposed model over the other competing models.Table 8
Information criterion for data set I.
DistributionAICBICCAICHQICGLL825.564834.120825.756829.040LN3826.723835.279826.916830.199LL2826.937835.641827.033829.254ExpW827.393835.949827.586830.869LL3827.458836.014827.651830.934G3831.955840.511832.148835.431W2832.163837.868832.259834.481W3832.665841.221832.858836.141Burr XII910.959916.663911.055913.276Table 9
MLE estimators of the model parameters, the log-likelihood, and goodness-of-fit statistics for data set I.
DistributionsEstimates (SEs)ℓW∗A∗K−S (p value)GLL (α, kη)α = 1.410 (0.174)−409.780.0190.1280.034(0.999)k = 0.134 (0.017)η = 0.077 (0.038)ExpW (α, kλ)α = 0.275 (0.146)−410.700.0450.2910.044(0.967)k = 0.676 (0.136)λ = 2.636 (1.161)LL3 (α,β,γ)α = 0.535 (0.061)−410.730.0190.1350.038(0.993)β = 1.863 (0.106)μ = −0.293 (0.358)LN3 (α,β,γ)α = 0.877 (0.090)−410.360.0170.1150.029(0.998)β = 1.925 (0.111)μ = −0.623 (0.372)G3 (α,β,γ)α = 1.098 (0.134)−412.980.1250.7780.067(0.618)β = 8.424 (1.238)μ = 0.075 (0.018)W3 (α,β,γ)α = 1.031 (0.072)−413.330.1340.8390.080(0.387)β = 9.743 (0.908)μ = 0.077 (0.013)W2 (α,β)α = 1.049 (0.068)−414.080.1310.7840.071(0.545)k = 9.576 (0.854)BXII (α,β)α = 2.342 (0.356)−453.480.7524.5640.251(<0.005)k = 0.233 (0.040)LL2 (α,β)α = 0.578 (0.043)−411.470.0430.3100.041(0.984)k = 1.805 (0.088)The estimated pdf and CDF of the proposed distribution corresponding to the real-world data set are shown in Figure12 and the Kaplan–Meier and PP plots for the proposed distribution are shown in Figure 13.Figure 12
Estimated pdf and CDF of the GLL distribution corresponding to data set I.Figure 13
PP and Kaplan–Meier plots of the GLL distribution corresponding to data set I.


### 6.2.1. TTT Plot

The total time test (TTT) plot plays a central role in determining the best model to fit the given data in terms of the hazard rates. This plot depicts the various forms of the hazard rate. A straight line on the TTT plot indicates that the given data has a constant hazard rate. If the plot is convex, the hazard rate will be decreased; if it is concave, the hazard rates will be increased. The plot for the bathtub shape is first convex and then concave. Similarly, if the hazard rate has an inverted bathtub shape, it will increase first (or concave) and then decrease (or convex). The TTT plot is calculated by using the following formula:(68)Grn=∑i=1rxi:n+n−rxi:n∑i=1rxi:n,r=xi:n=1,2,…,n,where xi:n are the order statistics.The TTT and box plots of the data set are presented in Figure14. These plots indicate that the empirical hazard rate function of the 1st data set is bathtub shape, monotonically increasing.Figure 14
TTT and box plots of data set I.The estimated fitted pdfs and CDFs of data set I for the competitive models are shown in Figure15.Figure 15
Some estimated fitted densities and cumulative functions of data set I.
(a)(b)


## 6.2.1. TTT Plot

The total time test (TTT) plot plays a central role in determining the best model to fit the given data in terms of the hazard rates. This plot depicts the various forms of the hazard rate. A straight line on the TTT plot indicates that the given data has a constant hazard rate. If the plot is convex, the hazard rate will be decreased; if it is concave, the hazard rates will be increased. The plot for the bathtub shape is first convex and then concave. Similarly, if the hazard rate has an inverted bathtub shape, it will increase first (or concave) and then decrease (or convex). The TTT plot is calculated by using the following formula:(68)Grn=∑i=1rxi:n+n−rxi:n∑i=1rxi:n,r=xi:n=1,2,…,n,where xi:n are the order statistics.The TTT and box plots of the data set are presented in Figure14. These plots indicate that the empirical hazard rate function of the 1st data set is bathtub shape, monotonically increasing.Figure 14
TTT and box plots of data set I.The estimated fitted pdfs and CDFs of data set I for the competitive models are shown in Figure15.Figure 15
Some estimated fitted densities and cumulative functions of data set I.
(a)(b)


## 7. Bayesian Model Formulation

Given a set of datax=x1,x2,…,xn from GLL α,k,η, the likelihood function of the model is given by(69)Lα,k,η|x=αkn∏i=1nkxiα−1∏i=1n1+ηxiα−kα/ηα+1.The Bayesian model is built by specifying the prior distribution for the model parametersα,kandη and then multiplying with the likelihood function Lα,k,η|x for the given data x=x1,x2,…,xn to obtain the posterior distribution function using the Bayes theorem. The prior distribution of α,kandη is denoted as pα,k,η.The joint posterior is(70)pα,k,η|x∝Lα,k,η|xpα,k,η.


### 7.1. Prior Distribution

We assumed independent noninformative gamma priors for the parameters of the proposed model in this study due to the flexibility of gamma distributions in accommodating many possible shapes for the types of parameters involved in the proposed distribution. Furthermore, they enable efficient posterior calculations and the recovery of the noninformative distribution for each parameter. Many research papers in the literature consider taking these priors into account (see [28, 37–41]).For the model parameters, we assume independent gamma priors:α∼Ga1,b1,k∼Ga2,b2, and η∼Ga3,b3.(71)pα=b1a1Γa1αa1−1exp−b1α,α>0,a1>0,b1>0,pk=b2a2Γa2ka2−1exp−b2k,α>0,a2>0,b2>0,pη=b3a3Γa3ηa3−1exp−b3η,η>0,a3>0,b3>0.Hence, we have(72)pα,k,η=pαpkpη.


### 7.2. Posterior Distribution

The posterior expression can be obtained, up to proportionality, by multiplying the likelihood by the prior, and this can be written as(73)pα,k,η|x∝αa1+n−1ka2+n−1ηa3+n−1e−b1α+b2k+b3ηL1,where(74)L1=αkn∏i=1nkxiα−1∏i=1n1+ηxiα−kα/ηα+1.The posterior is complicated, and there are no closed-form inferences. As a result, we, propose using McMC techniques to simulate samples from the posterior, allowing for simple sample-based inferences.


### 7.3. Gibbs Sampler: Algorithm

Markov chains require a stationary distribution in order to perform Markov chain Monte Carlo calculations. These chains can be built in a variety of ways. Over the last decade, the following Monte Carlo sampling techniques for assessing high-dimensional posterior integrals have already been developed. Others are Metropolis-Hastings’s sampling, Monte Carlo importance sampling, Gibb’s sampling, and others. The most popular McMC sampling algorithm in the Bayesian survival inference computation literature is Gibbs’ sampling, which is primarily a special case of Metropolis-Hastings’s sampling. Gibb’s sampling is preferred in high-dimensional numerical computation.By using Gibbs’s sampling, we only need to know the full conditional distribution. To carry out Gibbs’s sampling, the basic scheme is as follows:Step 1: compute the posterior distribution, up to proportionality, and specify the full conditionals, using equation (71), of the model parameters α,ηandk as follows.(i)
Full conditional ofα given η,kandx:(75)pα|η,k,x∝αa1+n−1e−b1αL1.(ii)
Full conditional ofk given α,ηandx:(76)pk|α,η,x∝ka2+n−1e−b2kL1.(iii)
Full conditional ofη given α,kandx:(77)pη|α,k,x∝ηa3+n−1e−b3ηL1.Step 2: select an initial valueθ0=α0,k0,η0 to start the chain.Step 3: suppose that, at theith step, θ=α,η,k takes the value θi=αi,ki,ηi; then, from full conditionals, generate(78)αi+1,frompα|ki,ηi,x,ki+1,frompk|αi+1,ηi,x,ηi+1,frompη|αi+1,ki+1,x.Step 4: this completes a transition fromθi to θi+1.Step 5: repeat Step 3N times.


## 7.1. Prior Distribution

We assumed independent noninformative gamma priors for the parameters of the proposed model in this study due to the flexibility of gamma distributions in accommodating many possible shapes for the types of parameters involved in the proposed distribution. Furthermore, they enable efficient posterior calculations and the recovery of the noninformative distribution for each parameter. Many research papers in the literature consider taking these priors into account (see [28, 37–41]).For the model parameters, we assume independent gamma priors:α∼Ga1,b1,k∼Ga2,b2, and η∼Ga3,b3.(71)pα=b1a1Γa1αa1−1exp−b1α,α>0,a1>0,b1>0,pk=b2a2Γa2ka2−1exp−b2k,α>0,a2>0,b2>0,pη=b3a3Γa3ηa3−1exp−b3η,η>0,a3>0,b3>0.Hence, we have(72)pα,k,η=pαpkpη.


## 7.2. Posterior Distribution

The posterior expression can be obtained, up to proportionality, by multiplying the likelihood by the prior, and this can be written as(73)pα,k,η|x∝αa1+n−1ka2+n−1ηa3+n−1e−b1α+b2k+b3ηL1,where(74)L1=αkn∏i=1nkxiα−1∏i=1n1+ηxiα−kα/ηα+1.The posterior is complicated, and there are no closed-form inferences. As a result, we, propose using McMC techniques to simulate samples from the posterior, allowing for simple sample-based inferences.


## 7.3. Gibbs Sampler: Algorithm

Markov chains require a stationary distribution in order to perform Markov chain Monte Carlo calculations. These chains can be built in a variety of ways. Over the last decade, the following Monte Carlo sampling techniques for assessing high-dimensional posterior integrals have already been developed. Others are Metropolis-Hastings’s sampling, Monte Carlo importance sampling, Gibb’s sampling, and others. The most popular McMC sampling algorithm in the Bayesian survival inference computation literature is Gibbs’ sampling, which is primarily a special case of Metropolis-Hastings’s sampling. Gibb’s sampling is preferred in high-dimensional numerical computation.By using Gibbs’s sampling, we only need to know the full conditional distribution. To carry out Gibbs’s sampling, the basic scheme is as follows:Step 1: compute the posterior distribution, up to proportionality, and specify the full conditionals, using equation (71), of the model parameters α,ηandk as follows.(i)
Full conditional ofα given η,kandx:(75)pα|η,k,x∝αa1+n−1e−b1αL1.(ii)
Full conditional ofk given α,ηandx:(76)pk|α,η,x∝ka2+n−1e−b2kL1.(iii)
Full conditional ofη given α,kandx:(77)pη|α,k,x∝ηa3+n−1e−b3ηL1.Step 2: select an initial valueθ0=α0,k0,η0 to start the chain.Step 3: suppose that, at theith step, θ=α,η,k takes the value θi=αi,ki,ηi; then, from full conditionals, generate(78)αi+1,frompα|ki,ηi,x,ki+1,frompk|αi+1,ηi,x,ηi+1,frompη|αi+1,ki+1,x.Step 4: this completes a transition fromθi to θi+1.Step 5: repeat Step 3N times.


## 8. Bayesian Analysis

In this work, we assumed the independent gamma priors forα∼Ga1,b1,k∼Ga2,b2, and η∼Ga3,b3 with hyperparameter values a1=b1=a2=b2=a3=b3=1.0.


### 8.1. Convergence Diagnostics

The proposed model is built with the goal of calculating Bayesian estimates for GLL parameters using the McMC method. Due to the Ergodic property of the Markov chain, all inferences are based on the assumption that it will converge. Hence, the McMC convergence diagnostic is crucial. If the simulated sample gives an acceptable approximation for the posterior density, the inferences are correct. Several convergence diagnostic analyses are used to determine whether the chains have converged, including the following.


#### 8.1.1. Geweke’s Convergence Diagnostic

Geweke’s diagnostic, also called Geweke’sz-score diagnostic, focuses on comparing the first and last parts of a chain. It is, in fact, a frequentist comparison, of means, with 95 percent of the values falling between −2 and 2, as proposed by [42]. All three values of the three parameters for the three chains in Figure 16 are between −2 and 2.Figure 16
Geweke’s diagnostic plot for alpha, eta, and kappa parameters.


#### 8.1.2. Autocorrelation Diagnostics

The autocorrelation plot for the parameters is shown in Figure17.Figure 17
Autocorrelation plot for the alpha, eta, and kappa parameters.


#### 8.1.3. Heidelberger and Welch’s Convergence Diagnostic

Schruben [43] and Schruben et al. [44] proposed detecting nonstationarity in simulation output using a spectral analysis approach to estimate the sample mean variance. They applied the Cramer-von Mises statistic and Brownian bridge theory to test the null hypothesis of stationarity of the Markov chain.Heidelberger and Welch [45] applied the aforementioned test to introduce a comprehensive method for generating a confidence interval of a predetermined width for the mean of a parameter when the chain has an initial transient (a state when the algorithm has not reached stationarity yet). They computed a test statistic (based on the Cramer-von Mises test statistic) to reject or accept the null hypothesis that the Markov chain belongs to a stationary distribution. A single chain was subjected to diagnostic.


#### 8.1.4. Raftery and Lewis’s Diagnostic

Raftery and Lewis [46, 47] proposed “a method for a single chain that tests for chain convergence to the target distribution and estimates the run-lengths required to properly estimate quantiles of functions of the parameters.”In this study, we applied a quantile of interest (0.025), the desired level of accuracy of ±0.0005, and a probability of 0.95 to attain the indicated degree of accuracy.


#### 8.1.5. Brooks–Gelman–Rubin (BGR) Convergence Diagnostic

The fact that the lines for all of the parameters are close to 1 indicates convergence from BGR plots as shown in Figure18.Figure 18
BGR plots for alpha, eta, and kappa parameters.In this section, a summary of some common statistical convergence diagnostics tests is provided in Table10.Table 10
Summary of some statistical convergence diagnostic tests.
ParameterGeweke’s diagnosticRaftery and LewisHeidelberger-WelchPr > |z|Total no. of samp.p valueStationarity testHalfwidth testAlpha−1.199238230.072PassedPassedEta−0.571143380.690PassedPassedKappa0.414441060.980PassedPassed


#### 8.1.6. Ergodic Mean (Running Mean) Plot

The running mean, also known as the ergodic mean, is the average of all samples up to and including a specific iteration. It is used to observe the McMC chains’ convergence pattern. Figure19 shows a time-series graph of each parameter and it displays the running mean (or ergodic mean) plots for the three parameters of the GLL distribution. The running mean plots of alpha, eta, and kappa show that the chains converge to the values in Table 11 after N iterations.Figure 19
The ergodic mean plots for alpha, eta, and kappa.Table 11
Numerical summaries of posterior properties for the GLL model with gamma priors based on an McMC sample.
CharacteristicsChain 1Chain 2Chain 3αηkαηkαηkMean1.4440.0940.1441.4370.0930.1441.4410.0930.143SD0.1750.0410.0180.1720.0400.0180.1740.0400.018Naïve SE0.0020.0010.00020.0020.0010.00020.0020.0010.0002Time-series SE0.0030.00030.00030.0030.00030.00040.0030.0010.0003MC error0.0010.00040.00010.0010.00030.00020.0020.0020.0001Minimum0.9670.0030.0900.9330.0010.0900.9650.0010.0882.5th percentile1.1390.0270.1121.1340.0250.1121.1340.0240.112Q11.3190.0650.1311.3160.0640.1311.3160.0640.131Medium (Q2)1.4320.0900.1421.4250.0880.1431.4250.0900.142Q31.5500.1180.1551.5480.1170.1551.5530.1180.15597.5th percentile1.8250.1840.1821.8200.1800.1811.8200.1800.183Maximum2.2400.3070.2382.0870.2812.3922.2330.2900.233Mode1.4500.0900.1451.4500.0900.1451.3500.0900.145Variance0.0310.0020.00031.6360.0020.00030.0300.0020.0003Skewness0.4630.6550.4840.3720.5140.3750.3650.5240.423Kurtosis0.3720.8210.6120.1530.3440.3010.0440.4200.41295% credible interval(1.139, 1.825)(0.027, 0.184)(0.112, 0.182)(1.134, 1.820)(0.025, 0.180)(0.112, 0.181)(1.134, 1.820)(0.024, 0.180)(0.112, 0.183)95% HPD interval(1.113, 1.784)(0.021, 0.174)(0.112, 0.181)(1.104, 1.764)(0.021, 0.171)(0.110, 0.179)(1.107, 1.779)(0.019, 0.172)(0.108, 0.178)


### 8.2. Posterior Analysis

In this section, we present numerical and visual summaries of the posterior distribution for each of the three chains. The joint posterior distribution for the proposed model was estimated using the JAGS software [48]. For each proposed model, we ran three parallel chains with 50,000 iterations and a burn-in of 5,000. Chains were thinned by storing every fifth iteration to reduce autocorrelation in the sample. The use of various convergence diagnostic tools ensured convergence to the joint posterior.


#### 8.2.1. Numerical Summary

We have considered different quantities of interest and their numeric data based on an McMC sample of posterior properties for generalized log-logistic distribution. The McMC simulation results include the results of of the posterior mean, posterior standard deviation, naïve standard error, time-series standard error, Markov chain error, the posterior five-point summary statistics (minimum, lower quartile (Q1), median (Q2), upper quartile (Q3), and maximum), the posterior skewness, posterior kurtosis, 2.5th percentile, 97.5th percentile, and the credible interval followed by the highest probability density (HPD).The naïve standard error is defined as a measure of simulation error in the mean rather than posterior uncertainty.(79)naive SE=posterior SDn.The time-series SE adjusts the “naïve” SE for autocorrelation.


#### 8.2.2. Visual Summary

In this subsection, we have considered different graphs for a visual summary of the posterior properties; those include the box plot, density strip plots, histogram, and trace plots for the parameters. These graphs and plots provide a nearly complete picture of the parameters’ posterior uncertainty [49]. We applied the posterior sample αj,kjandηj, j=1,…,15000, to draw these graphs.(1) Box Plots. The boxes in Figure 20 represent interquartile ranges, and the line in the middle of each box is the median; the arms of each box extend to encompass the central 95 percent of the distribution, and their ends thus correspond to the 2.5 percent and 97.5 percent quartiles, respectively.Figure 20
The box plots for the alpha, eta, and kappa parameters.(2) Density and Histogram Plots. Histogram can provide information about the behaviour in the tails, skewness, data outliers, and the presence of multimodal behaviour. The graphs in Figure 21 can provide us with a nearly complete picture of the posterior uncertainty about the GLL parameters, while the graphs in Figure 22 show a comparison of the full density and partial density of the parameters.Figure 21
Kernel density estimate and the histogram plots for alpha, eta, and kappa parameters.Figure 22
Density plots for the parameters comparing the whole chains with their last parties.(3) Trace Plots. A trace plot, also known as “a time-series plot,” is a representation of the iteration number versus the value of the parameter drawn at each iteration. Because the plots do not show long-term increasing or decreasing trends but rather resemble a horizontal band in Figure 23, we can conclude that the chains have converged.Figure 23
The trace plots for alpha, eta, and kappa parameters.


## 8.1. Convergence Diagnostics

The proposed model is built with the goal of calculating Bayesian estimates for GLL parameters using the McMC method. Due to the Ergodic property of the Markov chain, all inferences are based on the assumption that it will converge. Hence, the McMC convergence diagnostic is crucial. If the simulated sample gives an acceptable approximation for the posterior density, the inferences are correct. Several convergence diagnostic analyses are used to determine whether the chains have converged, including the following.


### 8.1.1. Geweke’s Convergence Diagnostic

Geweke’s diagnostic, also called Geweke’sz-score diagnostic, focuses on comparing the first and last parts of a chain. It is, in fact, a frequentist comparison, of means, with 95 percent of the values falling between −2 and 2, as proposed by [42]. All three values of the three parameters for the three chains in Figure 16 are between −2 and 2.Figure 16
Geweke’s diagnostic plot for alpha, eta, and kappa parameters.


### 8.1.2. Autocorrelation Diagnostics

The autocorrelation plot for the parameters is shown in Figure17.Figure 17
Autocorrelation plot for the alpha, eta, and kappa parameters.


### 8.1.3. Heidelberger and Welch’s Convergence Diagnostic

Schruben [43] and Schruben et al. [44] proposed detecting nonstationarity in simulation output using a spectral analysis approach to estimate the sample mean variance. They applied the Cramer-von Mises statistic and Brownian bridge theory to test the null hypothesis of stationarity of the Markov chain.Heidelberger and Welch [45] applied the aforementioned test to introduce a comprehensive method for generating a confidence interval of a predetermined width for the mean of a parameter when the chain has an initial transient (a state when the algorithm has not reached stationarity yet). They computed a test statistic (based on the Cramer-von Mises test statistic) to reject or accept the null hypothesis that the Markov chain belongs to a stationary distribution. A single chain was subjected to diagnostic.


### 8.1.4. Raftery and Lewis’s Diagnostic

Raftery and Lewis [46, 47] proposed “a method for a single chain that tests for chain convergence to the target distribution and estimates the run-lengths required to properly estimate quantiles of functions of the parameters.”In this study, we applied a quantile of interest (0.025), the desired level of accuracy of ±0.0005, and a probability of 0.95 to attain the indicated degree of accuracy.


### 8.1.5. Brooks–Gelman–Rubin (BGR) Convergence Diagnostic

The fact that the lines for all of the parameters are close to 1 indicates convergence from BGR plots as shown in Figure18.Figure 18
BGR plots for alpha, eta, and kappa parameters.In this section, a summary of some common statistical convergence diagnostics tests is provided in Table10.Table 10
Summary of some statistical convergence diagnostic tests.
ParameterGeweke’s diagnosticRaftery and LewisHeidelberger-WelchPr > |z|Total no. of samp.p valueStationarity testHalfwidth testAlpha−1.199238230.072PassedPassedEta−0.571143380.690PassedPassedKappa0.414441060.980PassedPassed


### 8.1.6. Ergodic Mean (Running Mean) Plot

The running mean, also known as the ergodic mean, is the average of all samples up to and including a specific iteration. It is used to observe the McMC chains’ convergence pattern. Figure19 shows a time-series graph of each parameter and it displays the running mean (or ergodic mean) plots for the three parameters of the GLL distribution. The running mean plots of alpha, eta, and kappa show that the chains converge to the values in Table 11 after N iterations.Figure 19
The ergodic mean plots for alpha, eta, and kappa.Table 11
Numerical summaries of posterior properties for the GLL model with gamma priors based on an McMC sample.
CharacteristicsChain 1Chain 2Chain 3αηkαηkαηkMean1.4440.0940.1441.4370.0930.1441.4410.0930.143SD0.1750.0410.0180.1720.0400.0180.1740.0400.018Naïve SE0.0020.0010.00020.0020.0010.00020.0020.0010.0002Time-series SE0.0030.00030.00030.0030.00030.00040.0030.0010.0003MC error0.0010.00040.00010.0010.00030.00020.0020.0020.0001Minimum0.9670.0030.0900.9330.0010.0900.9650.0010.0882.5th percentile1.1390.0270.1121.1340.0250.1121.1340.0240.112Q11.3190.0650.1311.3160.0640.1311.3160.0640.131Medium (Q2)1.4320.0900.1421.4250.0880.1431.4250.0900.142Q31.5500.1180.1551.5480.1170.1551.5530.1180.15597.5th percentile1.8250.1840.1821.8200.1800.1811.8200.1800.183Maximum2.2400.3070.2382.0870.2812.3922.2330.2900.233Mode1.4500.0900.1451.4500.0900.1451.3500.0900.145Variance0.0310.0020.00031.6360.0020.00030.0300.0020.0003Skewness0.4630.6550.4840.3720.5140.3750.3650.5240.423Kurtosis0.3720.8210.6120.1530.3440.3010.0440.4200.41295% credible interval(1.139, 1.825)(0.027, 0.184)(0.112, 0.182)(1.134, 1.820)(0.025, 0.180)(0.112, 0.181)(1.134, 1.820)(0.024, 0.180)(0.112, 0.183)95% HPD interval(1.113, 1.784)(0.021, 0.174)(0.112, 0.181)(1.104, 1.764)(0.021, 0.171)(0.110, 0.179)(1.107, 1.779)(0.019, 0.172)(0.108, 0.178)


## 8.1.1. Geweke’s Convergence Diagnostic

Geweke’s diagnostic, also called Geweke’sz-score diagnostic, focuses on comparing the first and last parts of a chain. It is, in fact, a frequentist comparison, of means, with 95 percent of the values falling between −2 and 2, as proposed by [42]. All three values of the three parameters for the three chains in Figure 16 are between −2 and 2.Figure 16
Geweke’s diagnostic plot for alpha, eta, and kappa parameters.


## 8.1.2. Autocorrelation Diagnostics

The autocorrelation plot for the parameters is shown in Figure17.Figure 17
Autocorrelation plot for the alpha, eta, and kappa parameters.


## 8.1.3. Heidelberger and Welch’s Convergence Diagnostic

Schruben [43] and Schruben et al. [44] proposed detecting nonstationarity in simulation output using a spectral analysis approach to estimate the sample mean variance. They applied the Cramer-von Mises statistic and Brownian bridge theory to test the null hypothesis of stationarity of the Markov chain.Heidelberger and Welch [45] applied the aforementioned test to introduce a comprehensive method for generating a confidence interval of a predetermined width for the mean of a parameter when the chain has an initial transient (a state when the algorithm has not reached stationarity yet). They computed a test statistic (based on the Cramer-von Mises test statistic) to reject or accept the null hypothesis that the Markov chain belongs to a stationary distribution. A single chain was subjected to diagnostic.


## 8.1.4. Raftery and Lewis’s Diagnostic

Raftery and Lewis [46, 47] proposed “a method for a single chain that tests for chain convergence to the target distribution and estimates the run-lengths required to properly estimate quantiles of functions of the parameters.”In this study, we applied a quantile of interest (0.025), the desired level of accuracy of ±0.0005, and a probability of 0.95 to attain the indicated degree of accuracy.


## 8.1.5. Brooks–Gelman–Rubin (BGR) Convergence Diagnostic

The fact that the lines for all of the parameters are close to 1 indicates convergence from BGR plots as shown in Figure18.Figure 18
BGR plots for alpha, eta, and kappa parameters.In this section, a summary of some common statistical convergence diagnostics tests is provided in Table10.Table 10
Summary of some statistical convergence diagnostic tests.
ParameterGeweke’s diagnosticRaftery and LewisHeidelberger-WelchPr > |z|Total no. of samp.p valueStationarity testHalfwidth testAlpha−1.199238230.072PassedPassedEta−0.571143380.690PassedPassedKappa0.414441060.980PassedPassed


## 8.1.6. Ergodic Mean (Running Mean) Plot

The running mean, also known as the ergodic mean, is the average of all samples up to and including a specific iteration. It is used to observe the McMC chains’ convergence pattern. Figure19 shows a time-series graph of each parameter and it displays the running mean (or ergodic mean) plots for the three parameters of the GLL distribution. The running mean plots of alpha, eta, and kappa show that the chains converge to the values in Table 11 after N iterations.Figure 19
The ergodic mean plots for alpha, eta, and kappa.Table 11
Numerical summaries of posterior properties for the GLL model with gamma priors based on an McMC sample.
CharacteristicsChain 1Chain 2Chain 3αηkαηkαηkMean1.4440.0940.1441.4370.0930.1441.4410.0930.143SD0.1750.0410.0180.1720.0400.0180.1740.0400.018Naïve SE0.0020.0010.00020.0020.0010.00020.0020.0010.0002Time-series SE0.0030.00030.00030.0030.00030.00040.0030.0010.0003MC error0.0010.00040.00010.0010.00030.00020.0020.0020.0001Minimum0.9670.0030.0900.9330.0010.0900.9650.0010.0882.5th percentile1.1390.0270.1121.1340.0250.1121.1340.0240.112Q11.3190.0650.1311.3160.0640.1311.3160.0640.131Medium (Q2)1.4320.0900.1421.4250.0880.1431.4250.0900.142Q31.5500.1180.1551.5480.1170.1551.5530.1180.15597.5th percentile1.8250.1840.1821.8200.1800.1811.8200.1800.183Maximum2.2400.3070.2382.0870.2812.3922.2330.2900.233Mode1.4500.0900.1451.4500.0900.1451.3500.0900.145Variance0.0310.0020.00031.6360.0020.00030.0300.0020.0003Skewness0.4630.6550.4840.3720.5140.3750.3650.5240.423Kurtosis0.3720.8210.6120.1530.3440.3010.0440.4200.41295% credible interval(1.139, 1.825)(0.027, 0.184)(0.112, 0.182)(1.134, 1.820)(0.025, 0.180)(0.112, 0.181)(1.134, 1.820)(0.024, 0.180)(0.112, 0.183)95% HPD interval(1.113, 1.784)(0.021, 0.174)(0.112, 0.181)(1.104, 1.764)(0.021, 0.171)(0.110, 0.179)(1.107, 1.779)(0.019, 0.172)(0.108, 0.178)


## 8.2. Posterior Analysis

In this section, we present numerical and visual summaries of the posterior distribution for each of the three chains. The joint posterior distribution for the proposed model was estimated using the JAGS software [48]. For each proposed model, we ran three parallel chains with 50,000 iterations and a burn-in of 5,000. Chains were thinned by storing every fifth iteration to reduce autocorrelation in the sample. The use of various convergence diagnostic tools ensured convergence to the joint posterior.


### 8.2.1. Numerical Summary

We have considered different quantities of interest and their numeric data based on an McMC sample of posterior properties for generalized log-logistic distribution. The McMC simulation results include the results of of the posterior mean, posterior standard deviation, naïve standard error, time-series standard error, Markov chain error, the posterior five-point summary statistics (minimum, lower quartile (Q1), median (Q2), upper quartile (Q3), and maximum), the posterior skewness, posterior kurtosis, 2.5th percentile, 97.5th percentile, and the credible interval followed by the highest probability density (HPD).The naïve standard error is defined as a measure of simulation error in the mean rather than posterior uncertainty.(79)naive SE=posterior SDn.The time-series SE adjusts the “naïve” SE for autocorrelation.


### 8.2.2. Visual Summary

In this subsection, we have considered different graphs for a visual summary of the posterior properties; those include the box plot, density strip plots, histogram, and trace plots for the parameters. These graphs and plots provide a nearly complete picture of the parameters’ posterior uncertainty [49]. We applied the posterior sample αj,kjandηj, j=1,…,15000, to draw these graphs.(1) Box Plots. The boxes in Figure 20 represent interquartile ranges, and the line in the middle of each box is the median; the arms of each box extend to encompass the central 95 percent of the distribution, and their ends thus correspond to the 2.5 percent and 97.5 percent quartiles, respectively.Figure 20
The box plots for the alpha, eta, and kappa parameters.(2) Density and Histogram Plots. Histogram can provide information about the behaviour in the tails, skewness, data outliers, and the presence of multimodal behaviour. The graphs in Figure 21 can provide us with a nearly complete picture of the posterior uncertainty about the GLL parameters, while the graphs in Figure 22 show a comparison of the full density and partial density of the parameters.Figure 21
Kernel density estimate and the histogram plots for alpha, eta, and kappa parameters.Figure 22
Density plots for the parameters comparing the whole chains with their last parties.(3) Trace Plots. A trace plot, also known as “a time-series plot,” is a representation of the iteration number versus the value of the parameter drawn at each iteration. Because the plots do not show long-term increasing or decreasing trends but rather resemble a horizontal band in Figure 23, we can conclude that the chains have converged.Figure 23
The trace plots for alpha, eta, and kappa parameters.


## 8.2.1. Numerical Summary

We have considered different quantities of interest and their numeric data based on an McMC sample of posterior properties for generalized log-logistic distribution. The McMC simulation results include the results of of the posterior mean, posterior standard deviation, naïve standard error, time-series standard error, Markov chain error, the posterior five-point summary statistics (minimum, lower quartile (Q1), median (Q2), upper quartile (Q3), and maximum), the posterior skewness, posterior kurtosis, 2.5th percentile, 97.5th percentile, and the credible interval followed by the highest probability density (HPD).The naïve standard error is defined as a measure of simulation error in the mean rather than posterior uncertainty.(79)naive SE=posterior SDn.The time-series SE adjusts the “naïve” SE for autocorrelation.


## 8.2.2. Visual Summary

In this subsection, we have considered different graphs for a visual summary of the posterior properties; those include the box plot, density strip plots, histogram, and trace plots for the parameters. These graphs and plots provide a nearly complete picture of the parameters’ posterior uncertainty [49]. We applied the posterior sample αj,kjandηj, j=1,…,15000, to draw these graphs.(1) Box Plots. The boxes in Figure 20 represent interquartile ranges, and the line in the middle of each box is the median; the arms of each box extend to encompass the central 95 percent of the distribution, and their ends thus correspond to the 2.5 percent and 97.5 percent quartiles, respectively.Figure 20
The box plots for the alpha, eta, and kappa parameters.(2) Density and Histogram Plots. Histogram can provide information about the behaviour in the tails, skewness, data outliers, and the presence of multimodal behaviour. The graphs in Figure 21 can provide us with a nearly complete picture of the posterior uncertainty about the GLL parameters, while the graphs in Figure 22 show a comparison of the full density and partial density of the parameters.Figure 21
Kernel density estimate and the histogram plots for alpha, eta, and kappa parameters.Figure 22
Density plots for the parameters comparing the whole chains with their last parties.(3) Trace Plots. A trace plot, also known as “a time-series plot,” is a representation of the iteration number versus the value of the parameter drawn at each iteration. Because the plots do not show long-term increasing or decreasing trends but rather resemble a horizontal band in Figure 23, we can conclude that the chains have converged.Figure 23
The trace plots for alpha, eta, and kappa parameters.


## 9. Conclusions

This work introduced and presented results on the mathematical and statistical properties of the generalized log-logistic distribution. The GLL model contains several parametric survival submodels that could be used in a variety of statistics and probability applications. Statistical properties such as quantile function and their related results, moments and their related results,rth central moments, and residual and reversed residual life were derived. We have also considered the Bayesian and classical inference of the unknown parameters of the proposed distribution when the data is uncensored or complete. The Bayesian estimates are obtained using the Gibbs sampling method under the assumption of independent gamma priors on the shape and scale parameters. It is worth noting that when prior information is available, Bayes estimates clearly outperform maximum likelihood estimates. To assess the behaviour of the estimators, Monte Carlo simulations are run. The proposed distribution was also applied to a real-world data set and provided a better fit than its submodels and other common parametric survival distributions based on goodness-of-fit statistics, log-likelihood function, and information criterion values. As a result, we conclude that the GLL is the most appropriate model among the distributions considered and it is a very competitive model for explaining lifetime phenomena.This work has numerous potential extensions. In practice, for example, the presence of explanatory variables and long-term survivals is common. Furthermore, a regression model for both complete and incomplete (or censored) data could be beneficial. As a result, our framework can be further researched in these contexts. The GLL distribution could also be useful in studies comprising survival models such as accelerated failure time, competing risks, mixture cure, frailty, multiple states, and joint survival models, as well as longitudinal data.



---

*Source: 5820435-2021-10-11.xml*