Properties and Maximum Likelihood Estimation of the Novel Mixture of Fréchet Distribution

Abstract

In recent decades, there have been numerous endeavors to develop a novel category of survival distributions possessing enhanced flexibility through the extension of existing distributions. This article constructs and validates the statistical properties of a novel survival distribution in order to obtain an alternative distribution that is suitable for analyzing survival data by presenting the novel mixture of the Fréchet distribution along with statistical properties such as the probability density function (PDF), cumulative distribution function (CDF), $r^{th}$ ordinary moment, skewness, kurtosis, moment-generating function, mean, variance, mode, survival function, hazard function, and asymptotic behavior, as well as constructing the estimators of the unknown parameter by employing the expectation-maximization (EM) algorithm, and simulated annealing. Additionally, the performance of the proposed estimators was compared with bias, mean squared errors (MSE), and simulated variances, and given an illustrative example of the proposed distribution to the survival data set in order to show that the proposed distribution is appropriate for the right-skewed data. This will be extremely advantageous in survival analysis.

1. Introduction

Survival analysis, a branch of statistics pertaining to death or failure, encompasses various types of statistical methods to draw conclusions. These methods include (1) nonparametric statistics, such as the Kaplan–Meier estimator and the log-rank test; (2) semi-parametric statistics, exemplified by the Cox proportional hazards model; and (3) parametric statistics, which focus on simulating survival time probabilities. Analysts may deduce that the survival function has a parametric distribution. For instance, if the survival time adheres to an exponential distribution, the hazard rate will be constant. Conversely, if the survival time conforms to a log-normal distribution, the hazard rate varies with time. Consequently, estimation of the survival function, calculation of the confidence interval, and assessment of the relative risk ensue. The utilization of a parametric survival function proves highly effective when appropriate distributions and parameter values are selected. The parametric survival distribution serves as a comprehensive representation of various types of survival data.

Hundreds of univariate continuous distributions exist. Mixture models play a crucial role in numerous applications, including survival analysis, such as Farewell [1], Hunsberger et al. [2], and Joudaki et al. [3]. These models involve the combination of two or more statistical distributions to create a new distribution, thereby addressing various challenges encountered in the field. Recognizing the evident necessity for mixture distributions, extensive efforts have been devoted to integrating multiple well-established distributions and utilizing them to tackle relevant issues. In the context of complete samples, Niyomdecha and Srisuradetchai [4] introduce a novel continuous three-parameter survival distribution referred to as the Complementary Gamma Zero-Truncated Poisson distribution. The traits of the maximum value in a series of independently identical gamma-distributed random variables are combined with those of zero-truncated Poisson random variables in this distribution. Abdullahi and Phaphan [5] present a mixture of Nakagami distribution, accompanied by statistical properties and a comparative analysis of the efficacy of estimators utilizing the quasi-Newton method and simulated annealing. Nanuwong et al. [6] proposed the mixture Pareto distribution by combining a Pareto distribution and a length-biased Pareto distribution. This distribution was formulated based on the concept of a weighted two-component distribution. Further investigation pertaining to the mixture models can be found in the references [7,8].

The Fréchet distribution, alternatively referred to as the inverse Weibull distribution, holds extensive application in the field of survival modeling. Fréchet [9] initially introduced the Fréchet distribution, which subsequently underwent further exploration by Fisher and Tippett [10] as well as Gumbel [11]. Furthermore, Abbas and Yincai [12] conducted a comparative analysis of the scale parameter estimation for the Fréchet distribution, employing maximum likelihood, probability-weighted moments, and Bayes estimations. Nasir and Aslam [13] utilized a Bayesian technique to estimate the parameter of the Fréchet distribution. Reyad et al. [14] established QE-Bayes and E-Bayes estimates for the scale parameters associated with the Fréchet distribution. Recent developments have introduced various extensions to the Fréchet distribution. Notably, Mead et al. [15] proposed the beta exponential Fréchet distribution.

Consequently, this article paid special attention to developing a new survival distribution by employing the notion of a mixture distribution, which is based on the Fréchet distribution, to obtain a new alternative distribution with the value of the time-varying hazard rate and investigating the statistical properties of the new distribution, such as the probability density function, cumulative distribution function, $r^{th}$ ordinary moment, skewness, kurtosis, moment-generating function, mean, variance, mode, survival function, hazard function, asymptotic behavior, comparison of the estimators with several methods, and samples of applying to real data, which will be extremely useful in survival analysis.

2. The Fréchet Distribution

The Fréchet distribution, being a specific case of the generalized extreme value distribution [16], finds extensive application in the field of hydrology [17]. This distribution is commonly employed to model extreme events, including daily rainfall [18] and river discharges [19]. Moreover, the Fréchet distribution holds considerable significance in survival analysis utilizing experimental data from clinical research. Given its status as the inverse Weibull distribution, the Fréchet distribution exhibits properties akin to the Weibull distribution, such as time-varying hazard rates. As a result, the Fréchet distribution has been a subject of widespread discussion in the field of survival analysis.

Afify et al. [20] provides the probability density function (PDF), cumulative distribution function (CDF), and mean of the Fréchet distribution. The PDF of the Fréchet distribution described by

$$\begin{array}{c}\hfill g\left(x\right)=\delta {\lambda }^{\delta }{x}^{-\left(\delta +1\right)}{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }},x>0.\end{array}$$

(1)

Given that $\lambda >0$ represents a scale parameter and $\delta >0$ represents a shape parameter, the CDF associated with these parameters can be expressed as follows:

$$\begin{array}{c}\hfill G\left(x\right)=e^{-{\left(\frac{\lambda }{x}\right)}^{\delta }},x>0.\end{array}$$

(2)

Furthermore, the mean of the distribution can be determined as follows:

$$\begin{array}{c}\hfill E\left(X\right)=\left\{\begin{array}{c}\lambda \Gamma \left(1-\frac{1}{\delta }\right),\delta >1\\ \infty ;otherwise,\end{array}\right.\end{array}$$

(3)

where $\Gamma \left(1-\frac{1}{\delta }\right)$ represents a gamma function: $\Gamma \left(1-\frac{1}{\delta }\right)=\underset{0}{\overset{\infty }{\int }}{y}^{\left(1-\frac{1}{\delta }\right)-1}{e}^{-y}dy$.

3. The Length-Biased Fréchet Distribution

Within the framework presented by Hesham et al. [21], a length-biased Fréchet distribution was introduced along with its associated CDF, PDF, and mean. The specific form of the CDF can be expressed using Equation (4).

$$\begin{array}{c}\hfill G_L\left(x\right)=\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\gamma \left(1-\frac{1}{\delta },{\left(\frac{\lambda }{x}\right)}^{\delta }\right),x>0,\end{array}$$

(4)

where $\lambda >0$, $\delta >1$, $\Gamma$ represents a gamma function, and $\gamma$ represents an incomplete gamma function. The associated PDF can be expressed as follows:

$$\begin{array}{c}\hfill g_L\left(x\right)=\frac{\delta {\lambda }^{\delta -1}}{\Gamma \left(1-\frac{1}{\delta }\right)}{x}^{-\delta }{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}.\end{array}$$

(5)

Additionally, the distribution’s mean can be determined using the formula below:

$$\begin{array}{c}\hfill E_L\left(X\right)=\frac{\lambda \Gamma \left(1-\frac{2}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)},\delta >2.\end{array}$$

(6)

4. Theoretical Result

4.1. The Probability Density Function of the Novel Mixture Fréchet (NMF) Distribution

This subsection aims to construct a novel distribution by employing the notion of a mixture distribution. The proposed distribution will be a combination of two distinct distributions, namely the Fréchet distribution and the length-biased Fréchet distribution, with the probabilities weighted between the two distributions. Nevertheless, the inclusion of the weighted parameter, represented by “p” with $0\le p\le 1$, would lead to a more complex PDF, thereby increasing the difficulty of implementing it. Consequently, this article opts to employ the function of parameter $\lambda$ as a weighted parameter, thus rendering the PDF of the newly developed distribution with two parameters for enhanced flexibility in practical application. Hence, the PDF of the novel mixture Fréchet (NMF) distribution is defined as follows:

$$\begin{array}{c}\hfill f_{NMF}\left(x\right)=\frac{1}{\lambda +1}g\left(x\right)+\frac{\lambda }{\lambda +1}{g}_L\left(x\right),x>0,\end{array}$$

(7)

where $\lambda >0$ and $\frac{1}{\lambda +1}+\frac{\lambda }{\lambda +1}=1$. By substituting Equations (1) and (5) into Equation (7), the resulting expression is denoted as

$$\begin{array}{c}\hfill f_{NMF}\left(x\right)=\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta }{e}^{{-\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right),x>0,\lambda >0,\delta >1.\end{array}$$

(8)

Therefore, Equation (8) represents the PDF of the NMF distribution.

4.2. Validity Check of the NMF Distribution for a Proper Density Function

A PDF is considered valid if it satisfies the following conditions:

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }f\left(x\right)dx=1.\end{array}$$

(9)

In order to demonstrate the validity of the proposed NMF distribution as a PDF, the following steps are undertaken:

$$\begin{array}{c}\hfill {\int }_0^{\infty }\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta }{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)dx=1,\end{array}$$

(10)

let

$$\begin{array}{c}\hfill u={\left(\frac{\lambda }{x}\right)}^{\delta }⟹x=\frac{k^{\frac{1}{\beta }}}{\lambda },\end{array}$$

(11)

and

$$\begin{array}{c}\hfill dx=\frac{dk}{\lambda \beta {\left(\lambda x\right)}^{\beta -1}}.\end{array}$$

(12)

By substituting Equations (11) and (12) into Equation (10), the resulting expression can be obtained.

$$\begin{array}{ccc}\hfill {\int }_0^{\infty }{f}_{NMF}\left(x\right)dx& =& \frac{1}{\lambda +1}{\int }_0^{\infty }{k}^{\alpha -1+\frac{1}{\beta }}{e}^{-k}\left(\frac{\lambda }{k^{\frac{1}{\beta }}\Gamma \left(\alpha \right)}+\frac{1}{\Gamma \left(\alpha +\frac{1}{\beta }\right)}\right)dk\hfill \\ & =& \frac{1}{\lambda +1}\left(\frac{\lambda }{\Gamma \left(\alpha \right)}\Gamma \left(\alpha \right)+\frac{1}{\Gamma \left(\alpha +\frac{1}{\beta }\right)}\Gamma \left(\alpha +\frac{1}{\beta }\right)\right)\hfill \\ & =& \frac{\lambda }{(\lambda +1)}+\frac{1}{(\lambda +1)}\hfill \\ & =& 1.\hfill \end{array}$$

This demonstrates that the PDF defined in Equation (8) conforms to the properties of a valid probability density distribution. Figure 1 depicts the PDF of the novel mixture of Fréchet distribution for different parameter values. The displayed variety of shapes demonstrates the right-skewed nature of the NMF distribution. Additionally, being a family of asymmetric distributions, the NMF distribution proves to be valuable for analyzing skewed data, particularly data with a right-skewed distribution, such as survival data.

4.3. The Cumulative Density Function of the NMF Distribution

Let $G\left(x\right)$ and $G_L\left(x\right)$ represent the cumulative density function (CDF) of the Fréchet distribution and the length-biased Fréchet distribution, respectively. Consider a random variable X following the novel mixture Fréchet (NMF) distribution. The CDF for X in this instance can be written as follows:

$$\begin{array}{ccc}\hfill F_{NMF}\left(x\right)& =& {\int }_0^{x}f\left(t\right)dt\hfill \\ & =& {\int }_0^x\left(\frac{1}{\lambda +1}g\left(t\right)+\frac{\lambda }{\lambda +1}{g}_L\left(t\right)\right)dt\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& =& \frac{1}{\lambda +1}G\left(x\right)+\frac{\lambda }{\lambda +1}{G}_L\left(x\right).\hfill \end{array}$$

(14)

From Equation (14), the CDF of the novel mixture of Fréchet distribution can be expressed as

$$\begin{array}{ccc}\hfill F\left(x\right)& =& \frac{1}{\lambda +1}\left[e^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\right]+\frac{\lambda }{\lambda +1}\left[\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\gamma \left(1-\frac{1}{\delta },{\left(\frac{\lambda }{x}\right)}^{\delta }\right)\right]\hfill \\ & =& \frac{1}{\lambda +1}\left[e^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}+\frac{\lambda }{\Gamma \left(1-\frac{1}{\delta }\right)}\gamma \left(1-\frac{1}{\delta },{\left(\frac{\lambda }{x}\right)}^{\delta }\right)\right].\hfill \end{array}$$

(15)

4.4. The $r^{th}$ Ordinary Moment of the NMF Distribution

The NMF distribution’s $r^{th}$ ordinary moment is expressed as follows:

$$\begin{array}{c}\hfill {\mu }_r^{\prime }=E_M\left(X^r\right)={\int }_0^{\infty }{x}^{r}f\left(x\right)dx.\end{array}$$

(16)

Equation (19) gives the explicit expression for the $r^{th}$ ordinary moment of the NMF distribution upon inserting Equation (8) into Equation (16) and performing integration with respect to x.

$$\begin{array}{c}\hfill {\int }_0^{\infty }\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta +r}{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)dx.\end{array}$$

(17)

$$\begin{array}{c}\hfill {\mu }_r^{\prime }=\frac{1}{\lambda +1}E\left(X^r\right)+\frac{\lambda }{\lambda +1}{E}_L\left(X^r\right),\end{array}$$

(18)

where

$$\begin{array}{c}\hfill E\left(X^r\right)={\lambda }^r\Gamma \left(1-\frac{r}{\delta }\right),\end{array}$$

and

$$\begin{array}{c}\hfill E_L\left(X\right)=\frac{{\lambda }^r\Gamma \left(1-\frac{r+1}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}.\end{array}$$

$$\begin{array}{c}\hfill {\mu }_r^{\prime }=\frac{1}{\lambda +1}{\lambda }^r\Gamma \left(1-\frac{r}{\delta }\right)+\frac{1}{\lambda +1}\frac{{\lambda }^{r+1}\Gamma \left(1-\frac{r+1}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)},r=1,2,3,\cdots \end{array}$$

(19)

The following is the mathematical expression for the mean of the NMF distribution:

$$\begin{array}{c}\hfill E_{NMF}\left(X\right)=\frac{\lambda }{(\lambda +1)}\left[\Gamma \left(1-\frac{1}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{2}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right],\delta >2.\end{array}$$

(20)

The second moment of the NMF distribution, denoted as $E\left(X^2\right)$, can be derived from Equation (19) by setting the value of $r=2$.

$$\begin{array}{c}\hfill E_{NMF}\left(X^2\right)=\frac{{\lambda }^2}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{2}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{3}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right],\delta >3.\end{array}$$

(21)

The third moment of the NMF distribution, denoted as $E\left(X^3\right)$, can be obtained from Equation (19) by substituting $r=3$.

$$\begin{array}{c}\hfill E_{NMF}\left(X^3\right)=\frac{{\lambda }^3}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{3}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{4}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right],\delta >4.\end{array}$$

(22)

The fourth moment of the NMF distribution, denoted as $E\left(X^4\right)$, can be calculated by substituting $r=4$ in Equation (19).

$$\begin{array}{c}\hfill E_{NMF}\left(X^4\right)=\frac{{\lambda }^4}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{4}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{5}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right],\delta >5.\end{array}$$

(23)

Equation (19) at $r=1$ and $r=2$ and substituting into Equation (24) yields the variance of the NMF distribution.

$$\begin{array}{cc}\hfill Va{r}_{NMF}\left(X\right)=& E_M\left(X^2\right)-{\left[E_M\left(X\right)\right]}^2,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill Va{r}_{NMF}\left(X\right)=& \frac{{\lambda }^2}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{2}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{3}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]-\hfill \\ \hfill & {\left\{\frac{\lambda }{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{1}{\delta }\right)+\frac{\delta \Gamma \left(1-\frac{2}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]\right\}}^2,\delta >3.\hfill \end{array}$$

(25)

Hence, the standard deviation of the NMF distribution is:

$$\begin{array}{c}\hfill S{D}_{NMF}\left(X\right)=\sqrt{Va{r}_{NMF}\left(X\right)}.\end{array}$$

(26)

4.5. The Skewness and Kurtosis of the NMF Distribution

The novel mixture Fréchet (NMF) distribution’s skewness and kurtosis coefficients are provided as follows, respectively:

$$\begin{array}{c}\hfill {\Phi }_1=\frac{E_{NMF}\left(X^3\right)}{{\left(E_{NMF}\left(X^2\right)\right)}^{\frac{3}{2}}}=\frac{\frac{{\lambda }^3}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{3}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{4}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]}{{\left\{\frac{{\lambda }^2}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{2}{\delta }\right)+\frac{\delta \Gamma \left(1-\frac{3}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]\right\}}^{\frac{3}{2}}},\end{array}$$

(27)

and

$$\begin{array}{c}\hfill {\Phi }_2=\frac{E_{NMF}\left(X^4\right)}{{\left(E_{NMF}\left(X^2\right)\right)}^2}=\frac{\frac{{\lambda }^4}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{4}{\delta }\right)+\frac{\lambda \Gamma \left(1-\frac{5}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]}{{\left\{\frac{{\lambda }^2}{\left(\lambda +1\right)}\left[\Gamma \left(1-\frac{2}{\delta }\right)+\frac{\delta \Gamma \left(1-\frac{3}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right]\right\}}^2}.\end{array}$$

(28)

4.6. The Moment Generating Function of the NMF Distribution

The NMF distribution’s moment-generating function is provided by

$$\begin{array}{c}\hfill E_{NMF}\left(e^{Xt}\right)=M_X\left(t\right)=\sum _{r=0}^{\infty }\frac{t^r{E}_{NMF}\left(X^r\right)}{r!}.\end{array}$$

(29)

By substituting Equation (19) into (29), the NMF distribution’s moment-generating function is derived as presented in Equation (30).

$$\begin{array}{c}\hfill M_X\left(t\right)=\sum _{r=0}^{\infty }\frac{t^r}{r!}\left\{\frac{1}{\lambda +1}{\lambda }^r\Gamma \left(1-\frac{r}{\delta }\right)+\frac{1}{\lambda +1}\frac{{\lambda }^{r+1}\Gamma \left(1-\frac{r+1}{\delta }\right)}{\Gamma \left(1-\frac{1}{\delta }\right)}\right\}.\end{array}$$

(30)

4.7. The Mode of the NMF Distribution

By computing the derivative of the natural logarithm of Equation (8) with respect to x, setting it equal to zero, and solving for x, one is able to determine the mode of the NMF distribution. In this subsection, a nonlinear equation is obtained in Equation (32).

$$\begin{array}{cc}\hfill logf\left(x\right)=& log\left[\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}\right]-\delta log\left(x\right)-{\left(\frac{\lambda }{x}\right)}^{\delta }+log\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & \frac{{\left(\frac{\lambda }{x}\right)}^{\delta }\delta }{x}-\frac{\delta }{x}-\frac{1}{x^2\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)}=0.\hfill \end{array}$$

(32)

4.8. The Survival Function and the Hazard Rate Function of the NMF Distribution

Consider a continuous random variable, X, whose cumulative density function, $F\left(x\right)$, is specified on the range, $[0,\infty )$. The following is an expression for the survival function of X:

$$\begin{array}{c}\hfill S\left(x\right)=1-F\left(x\right).\end{array}$$

(33)

The survival function of the NMF distribution is obtained by inserting Equation (15) into Equation (33):

$$\begin{array}{c}\hfill S\left(x\right)=1-\frac{1}{\lambda +1}\left[e^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}+\frac{\lambda }{\Gamma \left(1-\frac{1}{\delta }\right)}\gamma \left(1-\frac{1}{\delta },{\left(\frac{\lambda }{x}\right)}^{\delta }\right)\right].\end{array}$$

(34)

Theoretically possible to define the hazard rate function of X as:

$$\begin{array}{c}\hfill hrf\left(x\right)=\frac{f\left(x\right)}{S\left(x\right)}.\end{array}$$

Consequently, the NMF distribution’s hazard rate function is given by

$$\begin{array}{c}\hfill hrf\left(x\right)=\frac{\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta }{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)}{1-\frac{1}{\lambda +1}\left[e^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}+\frac{\lambda }{\Gamma \left(1-\frac{1}{\delta }\right)}\Gamma \left(1-\frac{1}{\delta }{\left(\frac{\lambda }{x}\right)}^{\delta }\right)\right]}.\end{array}$$

(35)

4.9. Asymptotic Behavior of the NMF Distribution

The NMF distribution exhibits zero asymptotic behavior as x approaches infinity.

$$\begin{array}{c}\hfill \underset{x\to \infty }{lim}f\left(x\right)=\underset{x\to \infty }{lim}\left[\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta }{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)\right]=0.\end{array}$$

(36)

As x approaches $\lambda$:

$$\begin{array}{ccc}\hfill \underset{x\to \lambda }{lim}f\left(x\right)& =& \underset{x\to \infty }{lim}\left[\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}^{-\delta }{e}^{-{\left(\frac{\lambda }{x}\right)}^{\delta }}\left(\frac{1}{x}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)\right]\hfill \\ & =& \frac{\delta }{e\left(\lambda +1\right)}\left(\frac{1}{\lambda }+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right).\hfill \end{array}$$

(37)

4.10. Maximum Likelihood Estimation of the NMF Distribution

Maximum likelihood estimators will be utilized in this subsection to estimate the NMF distribution’s parameters. The likelihood function of the NMF distribution is defined as follows if $x_1$,…,$x_n$ represent a random sample of size n taken from the NMF distribution:

$$\begin{array}{ccc}\hfill L\left(x\right)& =& \prod _{i=1}^n\left\{\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}_i^{-\delta }{e}^{-{\left(\frac{\lambda }{x_i}\right)}^{\delta }}\left(\frac{1}{x_i}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)\right\},\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill \mathsf{\ell }(x)& =& log\prod _{i=1}^n\left\{\frac{\delta {\lambda }^{\delta }}{\left(\lambda +1\right)}{x}_i^{-\delta }{e}^{-{\left(\frac{\lambda }{x_i}\right)}^{\delta }}\left(\frac{1}{x_i}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)\right\}.\hfill \end{array}$$

(39)

Equation (38)’s natural logarithm has been employed to derive the log-likelihood function shown in Equation (40).

$$\begin{array}{ccc}\hfill \mathsf{\ell }(x)& =& nlog\left(\delta \right)+n\delta log\left(\lambda \right)-nlog\left(\lambda +1\right)-\delta \sum _{i=1}^{\infty }log\left(x_i\right)-\sum _{i=1}^{\infty }{\left(\frac{\lambda }{x_i}\right)}^{\delta }\hfill \\ & +& \sum _{i=1}^{\infty }log\left(\frac{1}{x_i}+\frac{1}{\Gamma \left(1+\frac{1}{\delta }\right)}\right).\hfill \end{array}$$

(40)

By taking the derivative of Equation (40) with respect to $\lambda$ and $\delta$ and then solving for each of those values, one can obtain the maximum likelihood estimators (MLEs).

$$\begin{array}{ccc}\hfill \frac{\partial \mathsf{\ell }(x)}{\partial \lambda }& =& \frac{n\delta }{\lambda }-\frac{n}{\lambda +1}-\sum _{i=1}^{\infty }\frac{{\left(\frac{\lambda }{x_i}\right)}^{\delta }\delta }{\lambda },\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill \frac{\partial \mathsf{\ell }(x)}{\partial \delta }& =& \frac{n}{\delta }+nlog\left(\lambda \right)-\sum _{i=1}^{\infty }log\left(x_i\right)-\sum _{i=1}^{\infty }{\left(\frac{\lambda }{x_i}\right)}^{\delta }log\left(\frac{\lambda }{x_i}\right)\hfill \\ & -& \frac{\Psi \left(1-\frac{1}{\delta }\right)}{{\delta }^2\Gamma \left(1-\frac{1}{\delta }\right)\left(\frac{1}{x_i}+\frac{1}{\Gamma \left(1-\frac{1}{\delta }\right)}\right)},\hfill \end{array}$$

(42)

where

$$\begin{array}{cc}\hfill \Psi \left(\delta \right)& =\frac{\partial \mathrm{l}\mathrm{n}(\Gamma (\delta \left)\right)}{\partial \delta }.\hfill \end{array}$$

Due to the nonlinearity of these equations, analytical solutions are not feasible, but iterative methods can be used to solve these numerically. This article proposes the utilization of the expectation-maximization (EM) algorithm and the simulated annealing to construct the MLEs for the NMF distribution.

4.10.1. Maximum Likelihood Estimation Employing the Simulated Annealing Algorithm

This article examines the MLEs for the unknown parameters of the NMF distribution. Analytical solutions for the MLEs are not attainable in Section 4.10. Therefore, in this part, the R optimization function, particularly the “optim” function, is employed for maximum likelihood estimation (MLE) using the simulated annealing. The steps of the Simulated Annealing Algorithm are as follows:

Step 1:

Give a initial value $x^{(k=0)}$, temperature T, number of iterations n, and desired accuracy $\epsilon$.

Step 2:

Pick a random value $x^{(k+1)}$ in the vicinity of $x^{\left(k\right)}$.

Step 3:

If $\Delta E<0$, where $\Delta E=f\left(x^{(k+1)}\right)-f\left(x^{\left(k\right)}\right)$, and $f\left(x\right)$ represents the objective function, then accept $x^{(k+1)}$. Otherwise, generate a random number $\alpha$ such that $\alpha \in (0,1)$. If $\alpha \le exp(-\Delta E/KT)$, where K is the Boltzmann constant, then accept $x^{(k+1)}$. Otherwise, return to Step 2.

Step 4:

If $|x^{(k+1)}-x^{\left(k\right)}|<\epsilon$ and T is sufficiently small, terminate the iterations. Otherwise, if the number of random number generations reaches n, decrease the value of T, let $k=k+1$, and go to Step 2. Otherwise, give $k=k+1$ and go to Step 2.

4.10.2. Maximum Likelihood Estimation Employing the EM-Algorithm

An EM algorithm is an iterative method employed to estimate unknown parameters in incomplete statistical models. The application of the EM algorithm encompasses two primary scenarios. The first arises when the data is incomplete due to observational process issues or limitations. The second arises when optimizing the likelihood function becomes challenging. The procedure for implementing the EM algorithm for the NMF distribution is outlined as follows:

Steps involved in the Expectation (E)-Step

• Calculate the log-likelihood function for an NMF distribution.

  $$\begin{array}{c}\hfill lnL\left(x\right)=\sum _{i=1}^{n}ln\left\{\frac{1}{\lambda +1}g\left(x_i\right)+\frac{\lambda }{\left(\lambda +1\right)}{g}_L\left(x_i\right)\right\}.\end{array}$$

  (43)
• Compute a complete log-likelihood function by assigning a missing value ${\kappa }_i$ in the function $lnL\left(x\right)$. The missing values ${\kappa }_i$ can take either 0 or 1. Thus, the complete random variable is denoted as $Y=(X;K)$, where $y_1,y_2,\cdots ,y_n$ represent the observations with $y_i=(t_i,{\kappa }_i)$ for $i=1,2,...,n$. Consequently, a complete log-likelihood function is written in:

  $$\begin{array}{cc}\hfill l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)& =\sum _{i=1}^n{\kappa }_{i}ln\left[\frac{1}{\lambda +1}g\left(t\right)\right]+\sum _{i=1}^n\left(1-{\kappa }_i\right)ln\left[\frac{\lambda }{\lambda +1}{g}_L\left(t\right)\right],\hfill \\ \hfill & =\sum _{i=1}^n{\kappa }_{i}ln\left[\frac{1}{\lambda +1}g\left(t\right)\right]+\sum _{i=1}^n\left(1-{\kappa }_i\right)ln\left[\frac{x_i}{(\lambda +1)\Gamma \left(1-\frac{1}{\delta }\right)}g\left(t\right)\right],\hfill \end{array}$$

  (44)

  where $\Theta =\{\lambda ,\delta \}$. The Equation (44) can be simplified by substituting Equation (1), resulting in the complete log-likelihood function, denoted as $l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)$, which is expressed as follows:

  $$\begin{array}{cc}\hfill l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)& =\sum _{i=1}^n\ln \left(x_i\right)-n\ln (\lambda +1)-n\ln \Gamma \left(1-\frac{1}{\delta }\right)\hfill \\ & -\sum _{i=1}^n{\kappa }_i\ln \left(x_i\right)+n\overline{\kappa }\Gamma \left(1-\frac{1}{\delta }\right)+n\ln \left(\delta \right)+n\delta \ln \left(\lambda \right)\hfill \\ & -(\delta +1)\sum _{i=1}^n\ln \left(x_i\right)-\sum _{i=1}^n{\left(\frac{\lambda }{x_i}\right)}^{\delta },\hfill \end{array}$$

  (45)

  where $\overline{\kappa }=\frac{1}{n}{\sum }_{i=1}^n{\kappa }_i$.
• Formulate the new complete log-likelihood function by eliminating constant expressions, resulting in the following expression:

  $$\begin{array}{cc}\hfill l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)& =-n\ln (\lambda +1)-n\ln \Gamma \left(1-\frac{1}{\delta }\right)+n\overline{\kappa }\Gamma \left(1-\frac{1}{\delta }\right)\hfill \\ & \hfill +n\ln \left(\delta \right)+n\delta \ln \left(\lambda \right)-(\delta +1)\sum _{i=1}^n\ln \left(x_i\right)-\sum _{i=1}^n{\left(\frac{\lambda }{x_i}\right)}^{\delta }.\end{array}$$

  (46)

A pseudo-log-likelihood function is derived at an E-step of an EM algorithm by replacing missing values with their respective expectations. Hence, the pseudo-log-likelihood function at the $k^{th}$ stage can be expressed as follows:

$$\begin{array}{cc}\hfill l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)& =-n\mathrm{l}\mathrm{n}(\lambda +1)-n\mathrm{l}\mathrm{n}\Gamma \left(1-\frac{1}{\delta }\right)+n{a}^{\left(k\right)}\Gamma \left(1-\frac{1}{\delta }\right)\hfill \\ & \hfill +n\mathrm{l}\mathrm{n}\left(\delta \right)+n\delta \mathrm{l}\mathrm{n}\left(\lambda \right)-(\delta +1)\sum _{i=1}^n\mathrm{l}\mathrm{n}\left(x_i\right)-\sum _{i=1}^n{\left(\frac{\lambda }{x_i}\right)}^{\delta },\end{array}$$

(47)

where $a^{\left(k\right)}=\frac{1}{n}{\sum }_{i=1}^n{a}_i^{\left(k\right)}$, and $a_i^{\left(k\right)}$ is given by

$$\begin{array}{c}\hfill a_i^{\left(k\right)}=\frac{\frac{1}{\lambda +1}g(x;{\lambda }^{\left(k\right)},{\delta }^{\left(k\right)})}{\frac{1}{\lambda +1}g(x;{\lambda }^{\left(k\right)},{\delta }^{\left(k\right)})+\frac{\lambda }{\lambda +1}{g}_L(x;{\lambda }^{\left(k\right)},{\delta }^{\left(k\right)})}.\end{array}$$

(48)

Steps involved in the Maximization (M)-Step

The M-step process involves iteratively increasing the number of function expressions. With each iteration, the values of $a^{\left(k\right)}$ and the estimated parameters ${\lambda }^{(k+1)},$ and ${\delta }^{(k+1)}$ will adjust. The process continues until the estimated values remain unchanged. Consequently, the MLEs for $\lambda ,$ and $\delta$ obtained via an EM algorithm are ${\lambda }^{(k+1)}$, and ${\delta }^{(k+1)}$, respectively, achieved by maximizing Equation (53). The initial values suggested in this article for the EM algorithm are ${\lambda }^{\left(0\right)}$ and ${\delta }^{\left(0\right)}$, which are as follows:

$$\begin{array}{cc}\hfill {\lambda }^{\left(0\right)}& ={\left(\frac{n}{t}\right)}^{\left(\frac{1}{\delta }\right)},\mathrm{where}t=\sum _{i=1}^n{\left(\frac{1}{t_i}\right)}^{\delta },\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\delta }^{\left(0\right)}& =2\mathrm{for}\mathrm{sample}\mathrm{size}\mathrm{is}\mathrm{small},\mathrm{and}\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\delta }^{\left(0\right)}& =1.5\mathrm{for}\mathrm{sample}\mathrm{size}\mathrm{is}\mathrm{large}.\hfill \end{array}$$

(51)

Steps of EM-Algorithm:

Step 1:

Generate a random sample $t_1,t_2,\cdots ,t_n$ according to the NMF distribution.

Step 2:

Set $k=0$ and compute the initial values ${\lambda }^{\left(0\right)}$ and ${\delta }^{\left(0\right)}$ as specified in Equations (49)–(51).

Step 3:

Calculate $a^{\left(k\right)}=\frac{1}{n}{\sum }_{i=1}^n{a}_i^{\left(k\right)}$ for $i=0,1,2,\cdots ,n$, when $a_i^{\left(k\right)}$ was given by Equation (52). For example, when $k=0$, we obtain the following:

$$\begin{array}{c}\hfill a_i^{\left(0\right)}=\frac{\frac{1}{\lambda +1}g(x;{\lambda }^{\left(0\right)},{\delta }^{\left(0\right)})}{\frac{1}{\lambda +1}g(x;{\lambda }^{\left(0\right)},{\delta }^{\left(0\right)})+\frac{\lambda }{\lambda +1}{g}_L(x;{\lambda }^{\left(0\right)},{\delta }^{\left(0\right)})}.\end{array}$$

(52)

Step 4:

Obtain the values of ${\lambda }^{(k+1)}$ and ${\delta }^{(k+1)}$ by maximizing Equation (53). For instance, when $k=0$, we obtain the following values:

$$\begin{array}{cc}\hfill l_{\mathrm{complete}}\left(\Theta \setminus y_1,y_2,\cdots ,y_n\right)& =-n\mathrm{l}\mathrm{n}(\lambda +1)-n\mathrm{l}\mathrm{n}\Gamma \left(1-\frac{1}{\delta }\right)+n{a}^{\left(0\right)}\Gamma \left(1-\frac{1}{\delta }\right)\hfill \\ & \hfill +n\mathrm{l}\mathrm{n}\left(\delta \right)+n\delta \mathrm{l}\mathrm{n}\left(\lambda \right)-(\delta +1)\sum _{i=1}^n\mathrm{l}\mathrm{n}\left(x_i\right)-\sum _{i=1}^n{\left(\frac{\lambda }{x_i}\right)}^{\delta }.\end{array}$$

(53)

Step 5:

If ${\lambda }^{(k+1)}={\lambda }^{\left(k\right)}$ and ${\delta }^{(k+1)}={\delta }^{\left(k\right)}$, then the algorithm stops. Otherwise, update $k=k+1$ and proceed to Step 3 and Step 4.

4.10.3. Assessment of the Efficacy of the Parameter Estimation

In this subsection, a series of simulations were performed to compare the outcomes of maximum likelihood estimators obtained using EM algorithms and simulated annealing. The utilization of Equations (49) and (50) as the initial value for the simulated annealing via “optim” function is favored in this context. The random number generator employed for generating samples from the NMF distribution followed an acceptance-rejection algorithm, utilizing a Fréchet distribution from a VGAM package in R program version 4.3.0. Each model was subjected to 500 repetitions. Sample sizes of $n=5,10,30,50$ were generated for the NMF distribution with parameters $\lambda =1.5,2.5$ and $\delta =2,3,4$. The resulting computations yielded six models for each method and sample size, as presented in

Table 1. The average estimations, biases, MES, and variance estimation simulation of EM estimators $\widehat{\lambda }$ and $\widehat{\delta }$ with a sample size of $n=5$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\delta }}$	Bias $\left(\widehat{\mathit{\lambda }}\right)$	Bias $\left(\widehat{\mathit{\delta }}\right)$	MSE $\left(\widehat{\mathit{\lambda }}\right)$	MSE $\left(\widehat{\mathit{\delta }}\right)$	VarSim $\left(\widehat{\mathit{\lambda }}\right)$	VarSim $\left(\widehat{\mathit{\delta }}\right)$
1.5	2	1.7470	2.7426	0.2470	0.7426	0.4008	2.1382	0.3398	1.5867
	3	1.5637	4.0928	0.0637	1.0928	0.0814	4.3311	0.0773	3.1369
	4	1.5621	5.3009	0.0621	1.3009	0.0455	7.8207	0.0416	6.1283
2.5	2	2.9805	2.7295	0.4805	0.7295	1.2316	1.9732	1.0007	1.4410
	3	2.6579	3.9894	0.1579	0.9894	0.2846	3.6459	0.2597	2.6669
	4	2.6156	5.7329	0.1156	1.7329	0.1347	13.7042	0.1214	10.7012

,

Table 2. The average estimations, biases, MES, and variance estimation simulation of simulated annealing estimators $\tilde{\lambda }$ and $\tilde{\delta }$ with a sample size of $n=5$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\tilde{\mathit{\lambda }}$	$\tilde{\mathit{\delta }}$	Bias $\left(\tilde{\mathit{\lambda }}\right)$	Bias $\left(\tilde{\mathit{\delta }}\right)$	MSE $\left(\tilde{\mathit{\lambda }}\right)$	MSE $\left(\tilde{\mathit{\delta }}\right)$	VarSim $\left(\tilde{\mathit{\lambda }}\right)$	VarSim $\left(\tilde{\mathit{\delta }}\right)$
1.5	2	3.1557	14.6394	1.6557	12.6394	16.4554	231.1034	13.7139	71.3496
	3	2.0213	15.5284	0.5213	12.5284	2.1113	225.2154	1.8395	68.2533
	4	1.9652	15.0794	0.4652	11.0794	1.4879	188.2035	1.2715	65.4498
2.5	2	5.6882	14.0020	3.1882	12.0020	36.7693	213.4220	26.6046	69.3742
	3	3.9382	14.6874	1.4382	11.6874	8.9652	220.4475	6.8969	83.8532
	4	3.7100	13.9510	1.2100	9.9510	5.6971	180.6055	4.2331	81.5837

,

Table 3. The average estimations, biases, MES, and variance estimation simulation of EM estimators $\widehat{\lambda }$ and $\widehat{\delta }$ with a sample size of $n=10$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\delta }}$	Bias $\left(\widehat{\mathit{\lambda }}\right)$	Bias $\left(\widehat{\mathit{\delta }}\right)$	MSE $\left(\widehat{\mathit{\lambda }}\right)$	MSE $\left(\widehat{\mathit{\delta }}\right)$	VarSim $\left(\widehat{\mathit{\lambda }}\right)$	VarSim $\left(\widehat{\mathit{\delta }}\right)$
1.5	2	1.5660	2.2744	0.0660	0.2744	0.0973	0.2963	0.0929	0.2210
	3	1.5444	3.4037	0.0444	0.4037	0.0381	0.9369	0.0361	0.7740
	4	1.5279	4.6493	0.0279	0.6493	0.0197	2.3606	0.0190	1.9390
2.5	2	2.6520	2.3669	0.1520	0.3669	0.3017	0.4242	0.2786	0.2896
	3	2.5847	3.3550	0.0847	0.3550	0.1212	0.9211	0.1141	0.7950
	4	2.5446	4.5805	0.0446	0.5805	0.0539	1.6688	0.0519	1.3318

,

Table 4. The average estimations, biases, MES, and variance estimation simulation of simulated annealing estimators $\tilde{\lambda }$ and $\tilde{\delta }$ with a sample size of $n=10$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\tilde{\mathit{\lambda }}$	$\tilde{\mathit{\delta }}$	Bias $\left(\tilde{\mathit{\lambda }}\right)$	Bias $\left(\tilde{\mathit{\delta }}\right)$	MSE $\left(\tilde{\mathit{\lambda }}\right)$	MSE $\left(\tilde{\mathit{\delta }}\right)$	VarSim $\left(\tilde{\mathit{\lambda }}\right)$	VarSim $\left(\tilde{\mathit{\delta }}\right)$
1.5	2	3.2963	16.4428	1.7963	14.4428	18.4601	273.3670	15.2334	64.7724
	3	2.2299	16.4617	0.7299	13.4617	5.0049	240.0290	4.4722	58.8125
	4	1.9610	17.9110	0.4610	13.9110	2.0165	248.4057	1.8039	54.8905
2.5	2	7.4968	15.4376	4.9968	13.4376	75.0960	254.3734	50.1281	73.8031
	3	5.3225	16.0267	2.8225	13.0267	28.4780	252.5156	20.5112	82.8207
	4	4.3101	15.7719	1.8101	11.7719	12.9169	227.9040	9.6403	89.3260

,

Table 5. The average estimations, biases, MES, and variance estimation simulation of EM estimators $\widehat{\lambda }$ and $\widehat{\delta }$ with a sample size of $n=30$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\delta }}$	Bias $\left(\widehat{\mathit{\lambda }}\right)$	Bias $\left(\widehat{\mathit{\delta }}\right)$	MSE $\left(\widehat{\mathit{\lambda }}\right)$	MSE $\left(\widehat{\mathit{\delta }}\right)$	VarSim $\left(\widehat{\mathit{\lambda }}\right)$	VarSim $\left(\widehat{\mathit{\delta }}\right)$
1.5	2	1.5302	2.1075	0.0302	0.1075	0.0302	0.0595	0.0293	0.0479
	3	1.5117	3.1312	0.0117	0.1312	0.0108	0.1637	0.0106	0.1465
	4	1.5042	4.1778	0.0042	0.1778	0.0058	0.3782	0.0057	0.3466
2.5	2	2.5583	2.1745	0.0583	0.1745	0.0856	0.0743	0.0822	0.0439
	3	2.5240	3.1078	0.0240	0.1078	0.0309	0.1500	0.0303	0.1384
	4	2.5174	4.1438	0.0174	0.1438	0.0158	0.3076	0.0155	0.2869

,

Table 6. The average estimations, biases, MES, and variance estimation simulation of simulated annealing estimators $\tilde{\lambda }$ and $\tilde{\delta }$ with a sample size of $n=30$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\tilde{\mathit{\lambda }}$	$\tilde{\mathit{\delta }}$	Bias $\left(\tilde{\mathit{\lambda }}\right)$	Bias $\left(\tilde{\mathit{\delta }}\right)$	MSE $\left(\tilde{\mathit{\lambda }}\right)$	MSE $\left(\tilde{\mathit{\delta }}\right)$	VarSim $\left(\tilde{\mathit{\lambda }}\right)$	VarSim $\left(\tilde{\mathit{\delta }}\right)$
1.5	2	5.0867	16.3140	3.5867	14.3140	45.6847	263.9548	32.8206	59.0629
	3	3.6042	17.4744	2.1042	14.4744	22.9393	272.4161	18.5115	62.9077
	4	2.6454	18.8461	1.1454	14.8461	8.7214	286.8080	7.4095	66.4026
2.5	2	13.3191	16.1026	10.8191	14.1026	213.7684	278.4768	96.7157	79.5939
	3	9.2326	15.6701	6.7326	12.6701	105.1613	244.4895	59.8331	83.9577
	4	7.6995	14.3182	5.1995	10.3182	64.7645	201.8541	37.7293	95.3889

,

Table 7. The average estimations, biases, MES, and variance estimation simulation of EM estimators $\widehat{\lambda }$ and $\widehat{\delta }$ with a sample size of $n=50$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\widehat{\mathit{\lambda }}$	$\widehat{\mathit{\delta }}$	Bias $\left(\widehat{\mathit{\lambda }}\right)$	Bias $\left(\widehat{\mathit{\delta }}\right)$	MSE $\left(\widehat{\mathit{\lambda }}\right)$	MSE $\left(\widehat{\mathit{\delta }}\right)$	VarSim $\left(\widehat{\mathit{\lambda }}\right)$	VarSim $\left(\widehat{\mathit{\delta }}\right)$
1.5	2	1.5185	2.0694	0.0185	0.0694	0.0167	0.0293	0.0164	0.0245
	3	1.5071	3.0634	0.0071	0.0634	0.0066	0.0885	0.0065	0.0845
	4	1.5038	4.0746	0.0038	0.0746	0.0033	0.1818	0.0033	0.1763
2.5	2	2.6443	2.1385	0.1443	0.1385	2.6075	0.0505	2.5867	0.0313
	3	2.5367	3.0881	0.0367	0.0881	0.0511	0.0968	0.0497	0.0890
	4	2.5067	4.0758	0.0067	0.0758	0.0087	0.1513	0.0087	0.1456

and

Table 8. The average estimations, biases, MES, and variance estimation simulation of simulated annealing estimators $\tilde{\lambda }$ and $\tilde{\delta }$ with a sample size of $n=50$.

$\mathit{\lambda }$	$\mathit{\delta }$	$\tilde{\mathit{\lambda }}$	$\tilde{\mathit{\delta }}$	Bias $\left(\tilde{\mathit{\lambda }}\right)$	Bias $\left(\tilde{\mathit{\delta }}\right)$	MSE $\left(\tilde{\mathit{\lambda }}\right)$	MSE $\left(\tilde{\mathit{\delta }}\right)$	VarSim $\left(\tilde{\mathit{\lambda }}\right)$	VarSim $\left(\tilde{\mathit{\delta }}\right)$
1.5	2	7.9441	16.3235	6.4441	14.3235	109.6222	273.1470	68.0953	67.9839
	3	5.0498	17.6010	3.5498	14.6010	44.2769	290.7988	31.6761	77.6082
	4	4.2309	17.1473	2.7309	13.1473	28.2579	254.7019	20.8003	81.8496
2.5	2	19.1666	15.4000	16.6666	13.4000	397.1666	255.1433	119.3912	75.5836
	3	14.1514	13.6071	11.6514	10.6071	224.4617	202.7200	88.7069	90.2098
	4	12.3039	12.5652	9.8039	8.5652	168.4808	178.4291	72.3640	105.0656

.

Upon reviewing all the results from Figure 2. The performance of the EM algorithm was remarkable, with estimated values for most parameters closely resembling the actual values. Moreover, the proposed EM algorithm demonstrated higher precision compared to the maximum likelihood estimates obtained through simulated annealing, as evidenced by reduced bias, lower mean squared error (MSE), and decreased variance estimation simulation.

5. Illustrative Example

The proposed distribution is applied to an actual dataset in this part. The dataset used in this analysis was collected from a clinical trial conducted by Freireich et al. [22], where patients received a placebo to evaluate the efficacy of 6-mercaptopurine (6-MP) in maintaining remission. Following the completion of the trial after a year, the following remission times were recorded and are expressed in weeks: 1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23.

Based on the results shown in Figure 3, the remission times of patients who received a placebo had a right-skewed distribution. In order to compare the goodness of fit, four right-skewed distributions—the Fréchet distribution, the length-biased Fréchet distribution, the mixture of Nakagami distribution [5], and the proposed mixture Fréchet distribution—are chosen.

While the parameters of the other candidate distributions are determined using maximum likelihood estimation utilizing simulated annealing, the parameters of the novel mixture Fréchet (NMF) distribution are estimated using the EM algorithm. The best model is the one that provides the smallest Akaike information criterion (AIC) value, which is used as the evaluation criterion.

Based on the findings presented in

Table 9. The MLE of the model’s parameters for patients who received a placebo’s times of remission.

Fitting Distribution	Estimate Parameters			AIC Statistics
	$\mathit{\lambda }$	$\mathit{\delta }$	p
Fréchet Distribution	15.50508	12.18451	-	5.58502
Length-biased Fréchet Distribution	30.18082	1.5	-	11.4393
NMF Distribution	2.050587	3.5	-	2.951101
The mixture of Nakagami Distribution	1.773775	1.452312	0.7	4.524097

, it is evident that the NMF distribution yields the lowest value of the AIC. This indicates that the NMF distribution outperforms the other potential distributions when using an AIC statistic as a measure of goodness-of-fit for this example data. Therefore, as indicated by Equations (20) and (26), the mean and standard deviation of the remission times observed in a group of 21 patients who received a placebo are 3.091147 weeks and 2.792774 weeks, respectively.

6. Conclusions and Discussion

This article presents the introduction of a novel survival distribution known as the novel mixture Fréchet (NMF) distribution. This distribution is characterized by its right-skewed distribution. The study explores various statistical properties of this newly proposed distribution and estimates its two parameters using both EM algorithms and simulated annealing. To assess the performance of both methods, a simulation study is conducted, involving twenty-four different combination scenarios. The illustrative examples of the proposed distribution are implemented using patient remission times data. The results reveal that the EM estimators exhibit greater efficiency compared to the simulated annealing estimators. Additionally, the NMF distribution demonstrates a better fit when compared to other candidate distributions, as indicated by the Akaike information criterion (AIC). Consequently, this article presents a novel right-skewed distribution that holds potential application in diverse areas, including extreme value analysis, survival analysis, and reliability analysis.

In future research, it is advisable to investigate interval estimation using different methods, such as [23,24], to further enhance the accuracy of the estimations.