Parameter Estimation of the Lomax Lifetime Distribution Based on Middle-Censored Data: Methodology, Applications, and Comparative Analysis

Abstract

The Lomax distribution has important applications in survival analysis, reliability engineering, insurance, finance, and other fields. Middle-censoring is an important censoring scheme, and data with middle-censoring will produce censoring in random intervals. This paper studies the parameter estimation of the Lomax distribution based on middle-censored data. The expectation–maximization algorithm is employed to compute the maximum likelihood estimates of the two unknown parameters of the Lomax distribution. After processing the data using the midpoint approach estimation, the parameter estimates are obtained by two computational methods: the Newton–Raphson iteration method and the fixed-point method. Moreover, the calculation methods for the asymptotic confidence intervals of the two parameters are provided, with the confidence interval coverage rate serving as one of the criteria for evaluating the estimation performance. In the Bayesian estimation aspect, the shape parameter is estimated using a Gamma prior distribution, and the Gibbs sampling method is employed for the solution. Finally, both simulation data and real data are used to compare the accuracy of the various estimation methods.

1. Introduction

The Lomax distribution, due to its flexibility in describing extreme events, is widely applied in the fields of survival analysis, reliability engineering, insurance, and finance. Middle-censoring is an important mechanism for data censoring and has extensive applications in various fields, including medical research, public health and epidemiology, economics and finance, and more. Although there is a substantial amount of literature on the parameter estimation of the Lomax distribution and middle-censored data separately, their combination has not yet been sufficiently explored. The following is a detailed elaboration on the relevant content.

In lifetime testing, data censoring often occurs, where the exact lifetime data of observed samples may only be partially available. Depending on the different types of censoring that occur in the samples, censoring types can be categorized into right censoring, left censoring, double censoring, and so on. In a study by Leiderman et al. [1], there was a double deletion model of how long it took infants to learn to complete specific tasks in the first year. Groeneboom and Wellner [2] proposed the interval-censorship model, and Schick and Yu [3] proposed the mixed case IC model.

The type of censoring we focus on in this paper is middle-censoring. Middle-censoring is a type of phenomenon in survival analysis, and this type of data was successively proposed by Huang [4] and Jammalamadaka [5] from two different perspectives. It refers to situations in practical research where continuous observation is not possible. If an event occurs in an interval not observed by researchers, only the fact that the subject experienced the event within that interval is known, but not the exact time of occurrence. In the middle-censoring scheme, for n sample units, the lifetimes of the sample units are denoted as $T_1,T_2,\dots ,T_n$. For the i-th unit (i = 1, 2, …, n), there exists a random censoring interval $[L_i,R_i]$, which follows an unknown bivariate distribution and is independent of the lifetime. The exact failure time $T_i$ of the i-th unit can only be observed if $T_i\notin [L_i,R_i]$. Otherwise, only the censored interval $[L_i,R_i]$ is observed.

Many studies have been conducted on middle-censoring for parametric models. Jammalamadaka and Iyer [6] considered the approximate self-consistency of middle-censored data. Yu, Wong and Li [7] established strong consistency, asymptotic normality, and asymptotic efficiency of the self-consistent estimator and the generalized MLE with middle-censored data. Srikanth K. Iyer [8] studied the maximum likelihood estimation of exponential distributions with middle-censored data and investigated its consistency and asymptotic normality properties. Bennett et al. [9] provided a proof that the EM algorithm converges in such cases. Wang et al. [10] studied the bivariate Weibull distribution and considered the statistical analysis of related competing risks models, obtaining the maximum likelihood estimation (MLE), midpoint approximation (MPA) estimation, and approximate confidence intervals (ACI) for the unknown parameters. Jammalamadaka et al. [11] analyzed discrete lifetime data that follow a geometric distribution that is subject to middle-censoring.

The Lomax distribution, often referred to as the second type of Pareto distribution, was first proposed by K.S. Lomax [12] in 1954. It plays an important role in the analysis of lifetime data in medicine, biological sciences, and engineering. In recent years, many articles have studied the extensions of this distribution and its related statistical properties. For example, Lun [13] conducted research and pointed out the connection between the Lomax lifetime distribution and the exponential distribution under specific conditions. Aljohani et al. [14] investigated the statistical inference of the stress–strength reliability $R=P(X>Y)$ of the Lomax distribution under accelerated life testing, using both classical and Bayesian methods for point estimation of R. Alsuhabi et al. [15] extended the Lomax lifetime distribution by proposing a new continuous lifetime model with four parameters, estimating the model parameters using classical methods (e.g., maximum likelihood estimation (MLE) and maximum product of spacings) and non-classical methods (e.g., Bayesian analysis).

Many researchers have studied the estimation methods for the parameters of the Lomax distribution under censoring. Cramer et al. [16] discussed the competing risks model of the Lomax distribution under progressive Type-II censoring, establishing the maximum likelihood estimation of the distribution parameters. Alfaer et al. [17] studied the Lomax lifetime distribution under a balanced joint progressive Type-I censoring scheme, providing point estimates and asymptotic confidence intervals for the model parameters using maximum likelihood estimation. Bayesian estimates and related credible intervals were obtained using MCMC techniques under independent Gamma priors. El-Sherpieny [18] studied the Power Lomax distribution under progressive Type-II hybrid censoring. The authors discussed parameter estimation for this distribution using maximum likelihood estimation (MLE) and the maximum product spacing method. Hasaballah et al. [19] conducted research on parameter estimation for the Lomax lifetime distribution under censored data, employing both Bayesian and non-Bayesian methods to estimate the parameters.

There is an extensive body of research on the Lomax distribution and middle-censoring, and the existing research results are shown in Figure 1.

However, there has been no research to date that establishes a connection between these two concepts. Therefore, this paper aims to address this issue by studying parameter estimation of the Lomax distribution based on middle-censored data and by comparing and analyzing different estimation methods.

In Section 2 of our paper, we introduce the Lomax distribution as a continuous probability distribution and conduct modeling for middle-censored data. We provide a detailed expression for the maximum likelihood estimation under middle-censoring for the Lomax distribution and discuss its properties, which are crucial for the subsequent estimation methods.

In Section 3, we delve into frequentist methods for parameter estimation. We propose the EM algorithm to address the maximum likelihood estimation (MLE) problem for the two unknown parameters of the Lomax distribution. We also provide two alternative solutions for midpoint approach estimation (MPA). Additionally, we calculate the approximate confidence intervals for these parameters, which are necessary for statistical inference.

In Section 4, we explore Bayesian estimation methods, which offer a different perspective on parameter estimation compared to frequentist methods. We use a Gamma distribution as the prior for the shape parameter of the Lomax distribution and apply Bayesian techniques to estimate the parameters. This approach allows us to incorporate prior knowledge and provides a complete posterior distribution for the parameters.

In Section 5, we present the results of simulation experiments and experiments with real datasets. These experiments are designed to assess the accuracy of the various parameter estimation methods developed in the previous sections and to demonstrate their practical utility in real-life scenarios.

Finally, we summarize the main findings of this paper and discuss the significance of this research. In addition, we propose directions for future research and potential solutions.

2. Model Description

2.1. Lomax Distribution

The exponential distribution, commonly encountered in daily life, has a thin tail and is suitable for describing memoryless random events. In contrast, the Lomax distribution is a heavy-tailed distribution, making it more appropriate for modeling data with heavy-tailed characteristics and for describing extreme events. Additionally, compared with the Weibull distribution, whose shape parameter can control the thickness of the tail, the Lomax distribution may be more suitable for describing certain specific types of heavy-tailed data. In practical applications, the Lomax distribution is widely used in areas such as survival analysis (e.g., the lifetime of power transformers), wealth distribution, and queuing service times. Suppose the random variable T follows Lomax distribution, then its probability density function (PDF) and cumulative distribution function (CDF) are given by the Functions (1) and (2).

(1)$f\left(t\right)={\displaystyle \frac{\theta }{\lambda }}{\left(1+{\displaystyle \frac{t}{\lambda }}\right)}^{-\left(\theta +1\right)},t>0,$

(2)$F\left(t\right)=1-{\left(1+{\displaystyle \frac{t}{\lambda }}\right)}^{-\theta },t>0.$

where $\mathsf{\lambda }$ is called the scale parameter, $\mathit{\theta }$ is called the shape parameter, $\mathsf{\lambda }$ > 0, $\mathit{\theta }$ > 0.

The probability density functions of the Lomax distribution with different parameters are shown in Figure 2.

The Lomax distribution can be regarded as a heavy-tailed version of the exponential distribution. While the tail of the exponential distribution decays exponentially, the tail of the Lomax distribution decays according to a power law, specifically ${(x+\lambda )}^{-(\theta +1)}$. When $\theta$ is small, the Lomax distribution has a heavier tail, making it more suitable for describing extreme events. Compared to the exponential distribution, the shape parameter $\theta$ and the scale parameter $\lambda$ of the Lomax distribution provide sufficient flexibility to adapt to different types of data distributions. Therefore, the Lomax distribution has broader applications than the exponential distribution in some cases.

2.2. The Middle-Censoring Model of Lomax

The type of censoring we focus on in this paper is middle-censoring. If a data point falls within a random interval and is not observable, middle-censoring occurs. The middle-censoring scheme can be described as follows.

From the Figure 3, we can see the schematic diagram of the middle-censored data type. For some data, we can only observe the interval form and cannot observe the exact values.

For middle-censored data, there are two modeling approaches. One, proposed by Huang [4], specifically describes that observers selectively observe at certain times. When the lifetime data ends at the observation point, the data is observed; when the end occurs between two observation times, the data can only be recorded as an interval with the two observation times as its endpoints. The other, proposed by Jammalamadaka [5], specifically describes that during the continuous observation of data, there will be intervals that are not observed. When the lifetime data ends within this interval, the data will not be observable. The two think from different angles but the final effect is the same, both are descriptions of middle-censored data. This paper adopts the second modeling approach for middle-censored data.

Assume that n items are tested, and the lifetimes of these items are $T_1$, …, $T_n$. For item i, there exists a random censored interval $[L_i,R_i]$, which follows some unknown bivariate distribution. $T_i$ can be observed only when $T_i\notin [L_i,R_i]$, otherwise it cannot be observed. Suppose ${\delta }_i=I(T_i\notin [L_i,R_i])$, where I (·) denotes the indicator function. Therefore, when ${\delta }_i=1$, the observed results are not censored, and we observe the actual value $t_i=T_i$. In this case, we did not observe $[L_i,R_i]$. On the other hand, when ${\delta }_i=0$, we only observe the truncated interval $[L_i,R_i]$. For item i, we observe the following Formula (3):

(3)$\left\{\begin{array}{c}\left(t_i,1\right)if{T}_i\notin \left[L_i,R_i\right]\hfill \\ \left(\left[L_i,R_i\right],0\right)otherwise.\hfill \end{array}\right.$

In this section, the distribution of the two ends of the censored interval is unknown. In the subsequent parameter estimation test, the exponential distribution with different mean values is used to generate the two ends of the censored data to manually censor the data.

3. Frequentist Estimation

In this section, we use the traditional frequency method and the approximation method to represent the point and interval estimation of the unknown parameters.

3.1. Maximum Likelihood Estimation

Assuming that the number of actual observed data size is ${\mathit{n}}_{\mathbf{1}}$, the number of middle-censored sample is ${\mathit{n}}_{\mathbf{2}}$, and the total number of sample is n. Let $t_i$ be the observed data point, and $L_i$ and $R_i$ be the endpoints of the observed middle-censoring interval.

(4)$\begin{array}{cc}\hfill & \left\{\left(t_1,1\right),\left(t_2,1\right),\cdots ,\left(t_{n_1},1\right),\right.\hfill \\ \hfill & \left.\left(\left[L_{n_1+1},R_{n_1+1}\right],0\right),\left(\left[L_{n_1+2},R_{n_1+2}\right],0\right),\cdots ,\left(\left[L_{n_1+n_2},R_{n_1+n_2}\right],0\right)\right\}.\hfill \end{array}$

The likelihood function of the shape parameter $\theta$ and the scale parameter $\lambda$ based on middle-censored data of the Lomax life distribution [20] is

(5)$L\left(\theta ,\lambda \right)=\prod _{i=1}^{n_1}f\left(t_i\right)\prod _{i=n_1+1}^{n_1+n_2}\left[F\left(R_i\right)-F\left(L_i\right)\right].$

Substituting into Functions (1) and (2), we can obtain Function (6).

(6)$L\left(\theta ,\lambda \right)={\displaystyle \frac{{\theta }^{n_1}}{{\lambda }^{n_1}}}\prod _{i=1}^{n_1}{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\left(\theta +1\right)}\prod _{i=n_1+1}^{n_1+n_2}\left[{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }-{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }\right].$

Taking the logarithm of Function (6), we obtain the log-likelihood function.

(7)$\begin{array}{cc}\hfill lnL\left(\theta ,\lambda \right)=& n_{1}ln\theta -n_{1}ln\lambda -\left(\theta +1\right)\sum _{i=1}^{n_1}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)\hfill \\ & +\sum _{i=n_1+1}^{n_1+n_2}ln\left[{\left(1+\frac{L_i}{\lambda }\right)}^{-\theta }-{\left(1+{\displaystyle \frac{R_i}{\lambda }}\right)}^{-\theta }\right].\hfill \end{array}$

By differentiating Function (7) with respect to the two parameters, we obtain Functions (8) and (9).

(8)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial lnL\left(\theta ,\lambda \right)}{\partial \theta }}=& {\displaystyle \frac{n_1}{\theta }}-\sum _{i=1}^{n_1}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)\hfill \\ & +\sum _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{{\left(\lambda +L_i\right)}^{\theta }ln\left(1+\frac{R_i}{\lambda }\right)-{\left(\lambda +R_i\right)}^{\theta }ln\left(1+\frac{L_i}{\lambda }\right)}{{\left(\lambda +R_i\right)}^{\theta }-{\left(\lambda +L_i\right)}^{\theta }}}.\hfill \end{array}$

(9)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial lnL\left(\theta ,\lambda \right)}{\partial \lambda }}=& -{\displaystyle \frac{n_1}{\lambda }}+\left(\theta +1\right)\sum _{i=1}^{n_1}{\displaystyle \frac{t_i}{{\lambda }^2+\lambda t_i}}+\hfill \\ & \sum _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{\frac{\theta }{\lambda }\left[R_i{\left(\lambda +L_i\right)}^{\theta +1}-L_i{\left(\lambda +R_i\right)}^{\theta +1}\right]}{{\left(\lambda +R_i\right)}^{\theta +1}\left(\lambda +L_i\right)-{\left(\lambda +L_i\right)}^{\theta +1}\left(\lambda +R_i\right)}}.\hfill \end{array}$

Equation (10) is obtained by setting the two functions to zero and joining them.

(10)$\left\{\begin{array}{c}{\displaystyle \frac{n_1}{\theta }}-\sum _{i=1}^{n_1}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)+\sum _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{{\left(\lambda +L_i\right)}^{\theta }ln\left(1+\frac{R_i}{\lambda }\right)-{\left(\lambda +R_i\right)}^{\theta }ln\left(1+\frac{L_i}{\lambda }\right)}{{\left(\lambda +R_i\right)}^{\theta }-{\left(\lambda +L_i\right)}^{\theta }}}=0\hfill \\ -{\displaystyle \frac{n_1}{\lambda }}+\left(1+\theta \right)\sum _{i=1}^{n_1}{\displaystyle \frac{t_i}{{\lambda }^2+\lambda t_i}}+\sum _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{\frac{\theta }{\lambda }\left[R_i{\left(\lambda +L_i\right)}^{\theta +1}-L_i{\left(\lambda +R_i\right)}^{\theta +1}\right]}{{\left(\lambda +R_i\right)}^{\theta +1}\left(\lambda +L_i\right)-{\left(\lambda +L_i\right)}^{\theta +1}\left(\lambda +R_i\right)}}=0.\hfill \end{array}\right.$

The strong consistency, asymptotic normality, and asymptotic efficiency of the maximum likelihood estimation under middle-censoring have been proven [7]. When the data are complete, we can directly solve for the shape parameter $\theta$ and the scale parameter $\lambda$ using this expression, and the existence and uniqueness of the solution can be proven [21]. MLE with middle-censored data also possesses similar properties, but due to the presence of middle-censoring data, it is not feasible to directly derive the maximum likelihood estimates of the parameters. To solve this problem, we will use the EM algorithm to process the missing data.

3.2. Expectation–Maximization Algorithm

The expectation–maximization (EM) algorithm [22] is an iterative optimization algorithm commonly used for models with latent variables or missing data. In the parameter estimation of middle-censored data, the EM algorithm has significant advantages, as it can effectively handle missing data issues and gradually improve the accuracy of parameter estimation through iterative optimization. In most cases, the algorithm converges to a local optimal solution [23]. For many practical problems, the local optimal solution is sufficiently good. Its simplicity, flexibility, and numerical stability make it the preferred method for many practical problems. In particular, for exponential family distributions (such as the Lomax distribution) [22], the E-step and M-step calculations usually have closed-form solutions.

The core idea of the EM algorithm is to gradually improve the estimation of parameters by constructing an expectation step (E-step) and a maximization step (M-step). At the same time, let $t_i^{*}$ ($i=n_1+1,n_1+2,\dots ,n$) be the failure time in the interval $[L_i,R_i]$ at iteration step k, then the complete data can be obtained.

(11)$\left\{\left(t_1,1\right),\left(t_2,1\right),\dots ,\left(t_{n_1},1\right),\left(t_{n_1+1}^{*},0\right),\dots ,\left(t_{n_1+n_2}^{*},0\right)\right\}.$

By completing the complete data, the likelihood function can be obtained as Function (12).

(12)$L_c\left(\theta ,\lambda \right)=\prod _{i=1}^{n_1}f\left(t_i\right)\prod _{i=n_1+1}^{n_1+n_2}f\left(t_i^{*}\right).$

The log-likelihood function is

(13)$ln{L}_c\left(\theta ,\lambda \right)=nln\theta -nln\lambda -\left(\theta +1\right)\sum _{i=1}^{n_1}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)-\left(\theta +1\right)\sum _{i=n_1+1}^{n_1+n_2}ln\left(1+{\displaystyle \frac{t_i^{*}}{\lambda }}\right).$

The expected values of the shape parameter and scale parameter are expressed as Equation (14) and Equation (15), respectively, where $f(t;\theta ,\lambda )$ and $F(t;\theta ,\lambda )$ represent the probability density function and cumulative distribution function of the Lomax distribution under specific parameters.

(14)$E_{\theta ,\lambda }\left[ln\left(1+{\displaystyle \frac{T_i}{\lambda }}\right)\left|T_i\in \right.\left(L_i,R_i\right)\right]={\displaystyle \frac{\underset{L_i}{\overset{R_i}{\int }}ln\left(1+\frac{t}{\lambda }\right)f\left(t;\theta ,\lambda \right)dt}{F\left(R_i;\theta ,\lambda \right)-F\left(L_i;\theta ,\lambda \right)}},$

(15)$E_{\theta ,\lambda }\left[{\displaystyle \frac{T_i}{{\lambda }^2+\lambda T_i}}\left|T_i\in \right.\left(L_i,R_i\right)\right]={\displaystyle \frac{\underset{L_i}{\overset{R_i}{\int }}\frac{t}{{\lambda }^2+\lambda T_i}f\left(t;\theta ,\lambda \right)dt}{F\left(R_i;\theta ,\lambda \right)-F\left(L_i;\theta ,\lambda \right)}}.$

The EM algorithm is frequently used to handle parameter estimation problems under middle-censored data [24], and the corresponding algorithm steps are as follows (Algorithm 1):

Algorithm 1: EM Algorithm for parameter estimation.

1:. Initialize the parameter values ${\lambda }^{\left(0\right)}$ and ${\theta }^{\left(0\right)}$ using an exponential distribution, and set the tolerance threshold $\epsilon =0.05$. 2:. E-step: During the k-th iteration, calculate the expected values. $\begin{array}{c}{E}_{1i}^{\left(k\right)}=E_{{\theta }^{(k-1)},{\lambda }^{(k-1)}}\left[ln\left(1+{\displaystyle \frac{T_i}{\lambda }}\right)\mid T_i\in (L_i,R_i)\right]\hfill \\ E_{2i}^{\left(k\right)}=E_{{\theta }^{(k-1)},{\lambda }^{(k-1)}}\left[ln\left(1+{\displaystyle \frac{T_i}{\lambda }}\right)\mid T_i=t_i\right]\hfill \\ E_{3i}^{\left(k\right)}=E_{{\theta }^{(k-1)},{\lambda }^{(k-1)}}\left[{\displaystyle \frac{T_i}{{\lambda }^2+\lambda T_i}}\mid T_i\in (L_i,R_i)\right]\hfill \\ E_{4i}^{\left(k\right)}=E_{{\theta }^{(k-1)},{\lambda }^{(k-1)}}\left[{\displaystyle \frac{T_i}{{\lambda }^2+\lambda T_i}}\mid T_i=t_i\right]\hfill \end{array}$ 3:. M-step: Update the parameter estimates by maximizing the likelihood function $L_c(\theta ,\lambda )$ based on the complete data. $\begin{array}{c}{\theta }^{\left(k\right)}={\displaystyle \frac{n}{{\sum }_{i=1}^{n_1}{E}_{1i}^{\left(k\right)}+{\sum }_{i=n_1+1}^{n_1+n_2}{E}_{2i}^{\left(k\right)}}}\hfill \\ {\lambda }^{\left(k\right)}={\displaystyle \frac{n}{({\theta }^{\left(k\right)}+1)\left({\sum }_{i=1}^{n_1}{E}_{3i}^{\left(k\right)}+{\sum }_{i=n_1+1}^{n_1+n_2}{E}_{4i}^{\left(k\right)}\right)}}\hfill \end{array}$ 4:. Convergence Check: If the following conditions are met, stop the iteration; otherwise, return to Step 2. $\begin{array}{c}|{\theta }^{\left(k\right)}-{\theta }^{(k-1)}|<\epsilon \\ |{\lambda }^{\left(k\right)}-{\lambda }^{(k-1)}|<\epsilon \end{array}$

The value of $\epsilon$ in the algorithm needs to be adjusted according to the actual situation. When higher precision for estimation is required, a smaller $\epsilon$ can be chosen. When higher computational efficiency is desired, a smaller $\epsilon$ can also be chosen.

By iterating through the above EM algorithm, the estimated values of $\lambda$ and $\theta$ can be obtained.

3.3. Midpoint Approach Estimation

Besides the EM estimation, the midpoint approach estimation (MPA) is also a method used for handling censored data [25], which is especially common in survival analysis and reliability engineering. This method replaces the censored data with the midpoint of the censored interval, thereby performing parameter estimation.

For censored data, we use the average of the two endpoints of the censored interval to replace the censored data.

(16)$\widehat{t_i}={\displaystyle \frac{L_i+R_i}{2}}.$

(17)$\left\{\left(t_1,1\right),\left(t_2,1\right),\dots ,\left(t_{n_1},1\right),\left(\widehat{t_{n_1+1}},0\right),\dots ,\left(\widehat{t_{n_1+n_2}},0\right)\right\}.$

After obtaining the complete data using the midpoint estimation method, the estimates of the shape and scale parameters can be obtained by solving the equations using the Newton–Raphson iteration method and the fixed-point method. For the complete data, the likelihood function is given by Function (18).

(18)$L_M\left(\theta ,\lambda \right)=\prod _{i=1}^{n_1}f\left(t_i\right)\prod _{i=n_1+1}^{n_1+n_2}f\left(\widehat{t_i}\right)=\prod _{i=1}^{n_1}{\displaystyle \frac{\theta {\lambda }^{\theta }}{{\left(t_i+\lambda \right)}^{\theta +1}}}\prod _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{\theta {\lambda }^{\theta }}{{\left(\widehat{t_i}+\lambda \right)}^{\theta +1}}}.$

3.3.1. The Newton–Raphson Iteration Method for Solving Under the MPA Estimation

Since the censored interval data have been imputed by the midpoint values of the intervals, the statistical software can directly compute the solutions for the two parameters. Therefore, the Newton–Raphson iteration method can be employed to obtain the estimates of the scale parameter and the shape parameter.

The negative log-likelihood function is given by Function (19).

(19)$\begin{array}{cc}\hfill neg\text{\_}log\text{\_}likelihood& =-ln{L}_M\left(\theta ,\lambda \right)\hfill \\ & =-\left(\sum _{i=1}^{n_1}ln\left({\displaystyle \frac{\theta {\lambda }^{\theta }}{{\left(t_i+\lambda \right)}^{\theta +1}}}\right)+\sum _{i=n_1+1}^{n_1+n_2}ln\left({\displaystyle \frac{\theta {\lambda }^{\theta }}{{\left(\widehat{t_i}+\lambda \right)}^{\theta +1}}}\right)\right).\hfill \end{array}$

Using R program (version 4.3.2) (All the codes in the text are run with this device. The device models are as follows: LAPTOP-90NAADHP, equipped with an Intel(R) Core(TM) i7-10510U CPU running at 1.80 GHz with a boost speed of 2.30 GHz. It has 16.0 GB of RAM, with 15.8 GB available. The device ID is B0B39209-7A3F-4006-81CE-AB49F27A02E1, and the product ID is 00342-35692-61862-AAOEM. The system is a 64-bit operating system based on an x64 processor.) for calculation: Select appropriate initial values ${\lambda }_0$ and ${\theta }_0$, we can utilize the general optimization function optim to find the minimum of the negative log-likelihood function, thereby obtaining the estimates of the shape and scale parameters.

3.3.2. The Fixed-Point Method for Solving Under the MPA Estimation

After taking the logarithm of the likelihood Function (18) for the complete data and differentiating with respect to the two parameters, the resulting system of equations, where the derivatives are set to zero, is given as Equation (20).

(20)$\left\{\begin{array}{c}\frac{n}{\theta }-\left(\sum _{i=1}^{n_1}ln\left(1+\frac{t_i}{\lambda }\right)+\sum _{i=n_1+1}^{n_1+n_2}ln\left(1+\frac{t_i}{\lambda }\right)\right)=0\hfill \\ -\frac{n}{\lambda }+\left(1+\theta \right)\left(\sum _{i=1}^{n_1}\frac{t_i}{{\lambda }^2+\lambda t_i}+\sum _{i=n_1+1}^{n_1+n_2}\frac{\widehat{t_i}}{{\lambda }^2+\lambda \widehat{t_i}}\right)=0.\hfill \end{array}\right.$

The expression for the shape parameter $\theta$ in terms of the scale parameter $\lambda$ is obtained through transformation by Equation (20).

(21)$\theta \left(\lambda \right)={\displaystyle \frac{n}{{\sum }_{i=1}^{n_1}ln\left(1+\frac{t_i}{\lambda }\right)+{\sum }_{i=n_1+1}^{n_1+n_2}ln\left(1+\frac{\widehat{t_i}}{\lambda }\right)}}.$

Substituting Function (21) into Equation (20) results in

(22)${\displaystyle \frac{n}{\lambda }}-\left(1+\theta \left(\lambda \right)\right)\left[\sum _{i=1}^{n_1}{\displaystyle \frac{t_i}{{\lambda }^2+\lambda t_i}}+\sum _{i=n_1+1}^{n_1+n_2}{\displaystyle \frac{\widehat{t_i}}{{\lambda }^2+\lambda \widehat{t_i}}}\right]=0.$

To solve Equation (22), we use the fixed-point method. First, we define the function $g\left(\lambda \right)$ as

(23)$g\left(\lambda \right)={\displaystyle \frac{1}{n}}\left(1+{\displaystyle \frac{n}{{\sum }_{i=1}^{n_1}{Q}_i+{\sum }_{i=n_1+1}^{n_1+n_2}{\widehat{Q}}_i}}\right)\left(\sum _{i=1}^{n_1}{P}_i+\sum _{i=n_1+1}^{n_1+n_2}{\widehat{P}}_i\right).$

The symbols in this equation are defined as

(24)$Q_i=ln\left(1+{\displaystyle \frac{t_{\mathbf{i}}}{\lambda }}\right),{\widehat{Q}}_i=ln\left(1+{\displaystyle \frac{\widehat{t_{\mathbf{i}}}}{\lambda }}\right),P_i={\displaystyle \frac{t_i}{1+\frac{t_{\mathbf{i}}}{\lambda }}},\widehat{P_i}={\displaystyle \frac{\widehat{t_i}}{1+\frac{\widehat{t_{\mathbf{i}}}}{\lambda }}}.$

Select an appropriate initial value ${\lambda }_0$, and use the fixed-point iteration method in the R program (version 4.3.2) to solve the equation iteratively. This will yield a proper solution $\lambda$. Substituting $\lambda$ into Function (21) allows us to obtain the estimate of $\theta$.

3.4. The Estimation of Asymptotic Confidence Interval

Censored data, such as the introduction of missing data in the middle, may lead to biased parameter estimates. The censoring mechanism can affect the sample distribution, making it impossible to directly use standard asymptotic theory because the sample distribution is no longer a simple known distribution. Meanwhile, with censored data, increasing the sample size does not necessarily result in more reliable estimates, especially when the degree of censoring is high, as the sample may fail to reflect the true characteristics of the population.

The asymptotic confidence interval (ACI) is a statistical method used for estimating the confidence interval of a parameter. Unlike the classical confidence interval, the ACI is derived based on the properties when the sample size approaches infinity.

For the parameters $\theta$ and $\lambda$, it is impossible to obtain an exact confidence interval. Therefore, we use the asymptotic properties of maximum likelihood estimators to construct the ACI [26].

The observed Fisher information matrix for parameters $\theta$ and $\lambda$ is represented as

(25)$I\left(\lambda ,\theta \right)=\left[\begin{array}{cc}{I}_{11}& I_{12}\\ I_{21}& I_{22}\end{array}\right],$

(26)$I_{11}=-{\displaystyle \frac{{\partial }^{2}lnL\left(\theta ,\lambda \right)}{\partial {\lambda }^2}},I_{22}=-{\displaystyle \frac{{\partial }^{2}lnL\left(\theta ,\lambda \right)}{\partial {\theta }^2}},$

(27)$I_{12}=I_{21}=-{\displaystyle \frac{{\partial }^{2}lnL\left(\theta ,\lambda \right)}{\partial \lambda \partial \theta }}.$

The logarithm of the likelihood is $lnL(\theta ,\lambda )$. We calculate $I_{11}$, $I_{12}$, $I_{21}$, and $I_{22}$, and substitute the estimated values of $\widehat{\lambda }$ and $\widehat{\theta }$ into the calculation formula for the Fisher matrix. The inverse of the matrix yields Function (28).

(28)$W\left(\lambda ,\theta \right)={\left[\begin{array}{cc}{I}_{11}& I_{12}\\ I_{21}& I_{22}\end{array}\right]}^{-1}{|}_{\lambda =\widehat{\lambda },\theta =\widehat{\theta }}=\left[\begin{array}{cc}\widehat{W_{11}}& \widehat{W_{12}}\\ \widehat{W_{21}}& \widehat{W_{22}}\end{array}\right].$

Through the asymptotic distribution properties of the likelihood function, we know that $\left(\lambda -\widehat{\lambda }\right)/\sqrt{\widehat{W_{11}}},\left(\theta -\widehat{\theta }\right)\sqrt{\widehat{W_{22}}}$ both follow a standard normal distribution. Under a confidence level of $\gamma$, the asymptotic confidence intervals for the parameters $\theta$ and $\lambda$ can be obtained.

(29)$\left[\widehat{\lambda }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}},\widehat{\lambda }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}}\right],$

(30)$\left[\widehat{\theta }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}},\widehat{\theta }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}}\right].$

4. Bayes Estimation

4.1. The Posterior Distribution of the Parameters of the Lomax Distribution Under Middle-Censoring

When the scale parameter $\lambda$ is known, we choose the Gamma distribution with parameters a and b as the prior distribution for the shape parameter $\theta$ [27], and denote it as $\pi \left(\theta \right)$.

(31)$\pi \left(\theta \right)={\displaystyle \frac{b^a}{\mathsf{\Gamma }\left(a\right)}}{\theta }^{a-1}{e}^{-b\theta }.$

Since the Gamma distribution is the conjugate prior for the shape parameter of the Lomax distribution, it has desirable properties when selected as the prior distribution. The proof that the Gamma distribution is the conjugate prior for the shape parameter of the Lomax distribution has been included in the Appendix A to ensure the rigor of the paper.

Suppose the experimental data with middle-censoring are denoted as data; then, the likelihood function of data can be represented as $p(\mathrm{data}\mid \theta )$, and the joint distribution of data and $\theta$ is $h(\mathrm{data},\theta )$. By integrating $h(\mathrm{data},\theta )$ over the parameter space $\Theta$, we can obtain the marginal density of data, denoted as $m\left(\mathrm{data}\right)$.

(32)$\begin{array}{cc}\hfill \mathit{data}=& \left\{\left(t_1,1\right),\left(t_2,1\right),\cdots ,\left(t_{n_1},1\right),\right.\hfill \\ \hfill & \left.\left(\left[L_{n_1+1},R_{n_1+1}\right],0\right),\left(\left[L_{n_1+2},R_{n_1+2}\right],0\right),\cdots ,\left(\left[L_{n_1+n_2},R_{n_1+n_2}\right],0\right)\right\}.\hfill \end{array}$

(33)$p\left(data\right|\theta )=L(\theta ,\lambda )={\displaystyle \frac{{\theta }^{n_1}}{{\lambda }^{n_1}}}\prod _{i=1}^{n_1}{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\left(\theta +1\right)}\prod _{i=n_1+1}^{n_1+n_2}\left[{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }-{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }\right].$

(34)$\begin{array}{cc}\hfill h(data,\theta )& =p\left(data\right|\theta \left)\pi \right(\theta )\hfill \\ & ={\displaystyle \frac{{\theta }^{n_1}}{{\lambda }^{n_1}}}\prod _{i=1}^{n_1}{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\left(\theta +1\right)}\prod _{i=n_1+1}^{n_1+n_2}\left[{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }-{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }\right]{\displaystyle \frac{b^a}{\mathsf{\Gamma }\left(a\right)}}{\theta }^{a-1}{e}^{-b\theta }.\hfill \end{array}$

(35)$m\left(data\right)={\int }_{\Theta }h\left(data,\theta \right)d\theta ={\int }_{\Theta }p\left(data|\theta \right)\pi \left(\theta \right)d\theta .$

From Functions (34) and (35), we can derive the posterior distribution of the parameter $\theta$ given data.

(36)$\begin{array}{cc}\hfill \pi \left(\theta |data\right)& ={\displaystyle \frac{h\left(data,\theta \right)}{m\left(data\right)}}\hfill \\ & \propto {\theta }^n\prod _{i=1}^{n_1}{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }{\theta }^{a-1}{e}^{-b\theta }\prod _{i=n_1+1}^{n_1+n_2}\left[{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }-{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }\right].\hfill \end{array}$

Using the absolute error loss function, we can obtain the Bayesian estimate of the shape parameter $\theta$.

(37)$\mathrm{MSE}={\displaystyle \frac{1}{M}}\sum _{i=1}^M{(\theta -{\widehat{\theta }}_i)}^2.$

4.2. The Bayesian Estimation of the Shape Parameter Based on Gibbs Sampling When the Scale Parameter Is Known

Directly solving the posterior distribution by integration poses some difficulties. Therefore, we resort to Gibbs sampling for the relevant computations of the posterior distribution. Gibbs sampling is a sampling technique within the Markov Chain Monte Carlo methods and is widely used in Bayesian inference. It generates approximate samples from a target distribution by iteratively sampling from conditional distributions [28].

To utilize Gibbs sampling, we first need to establish the condition that, when the data are complete and the scale parameter is known, the conditional posterior density of the shape parameter follows a Gamma distribution. Let the complete data be denoted as ${\mathrm{data}}^{*}$.

Theorem 1. 

When the data are complete and the scale parameters are known, the posterior condition density of the shape parameters is Gamma distribution.

Proof. 

Through the likelihood function of the Lomax distribution, we can obtain the conditional posterior density of the shape parameter when the data are complete.    □

(38) $\pi \left(\theta |dat{a}^{*}\right)\propto {\theta }^n\prod _{i=1}^n{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }{\theta }^{a-1}{e}^{-b\theta }.$

Thus, when the scale parameter is known, the conditional posterior distribution of the shape parameter is

(39)$Gamma\left(n+a,b+\sum _{i=1}^{n}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)\right).$

Therefore, the conclusion is proved. The steps to obtain the shape parameter using Gibbs sampling when the scale parameter is known are shown in Algorithm 2.

Algorithm 2: Gibbs sampling method for estimating shape parameter $\theta$.

1:. Determine the initial value of the shape parameter ${\theta }_0$. Define the number of Gibbs sampling iterations as M and the number of initial samples to be discarded as Q. 2:. In the k-th iteration, the missing data in the interval are replaced by conditional expectations, and the replaced values are denoted as $t_i^{k\#}$, for $i=n_1+1,n_1+2,\dots ,n_1+n_2$. $\begin{array}{c}\hfill t_i^{k\#}=E\left(T\right|T\in (L_i,R_i);\lambda ,{\theta }_{k-1})={\int }_{L_i}^{R_i}{\displaystyle \frac{Tf(T;\lambda ,{\theta }_{k-1})dT}{F(R_i;\lambda ,{\theta }_{k-1})-F(L_i;\lambda ,{\theta }_{k-1})}}\end{array}$ 3:. Substitute the imputed data into the distribution, and then randomly draw a sample from the distribution to obtain ${\theta }_k$. $\begin{array}{c}\hfill Gamma\left(n+a,b+\sum _{i=1}^{n_1}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)+\sum _{i=n_1+1}^{n_1+n_2}ln\left(1+{\displaystyle \frac{t_i^{k\#}}{\lambda }}\right)\right).\end{array}$ 4:. Repeat Steps 2 and 3 M times to obtain M posterior conditional distribution-related sample data ${\theta }_k$. The estimate of the shape parameter $\theta$ based on the sum of squared loss function is $\begin{array}{c}\hfill {\widehat{\theta }}_{\mathrm{bayes}}={\displaystyle \frac{{\sum }_{k=Q+1}^{M-Q}{\theta }_k}{M-Q}}.\end{array}$

By following the above Gibbs sampling steps, the Bayesian estimate of the shape parameter can be obtained through programming when the scale parameter is known.

5. Simulation Experiments and Real Data Analysis

5.1. Analysis of Simulation Experiments

5.1.1. Analysis of the Solution Methods for MPA Estimation

Using the R program (version 4.3.2), we generate data with three different sets of parameters, where $\lambda$ is the scale parameter and $\theta$ is the shape parameter. For each data point, censoring intervals are randomly generated using an exponential distribution with different parameters to determine whether the data point is censored. We then perform parameter estimation using the aforementioned methods on the middle-censored data.

The estimated parameters are compared with the true parameters used for data generation. The accuracy of the parameter estimates is evaluated using bias and mean absolute error (MAE). The relevant formulas are as follows:

(40)$\begin{array}{cc}\hfill & L_{MAE}={\displaystyle \frac{1}{M}}\sum _{i=1}^M\left|\lambda -\widehat{{\lambda }_i}\right|,L_{Bias}={\displaystyle \frac{1}{M}}\sum _{i=1}^M\left(\lambda -\widehat{{\lambda }_i}\right),\hfill \\ & T_{MAE}={\displaystyle \frac{1}{M}}\sum _{i=1}^M\left|\theta -\widehat{{\theta }_i}\right|,T_{Bias}={\displaystyle \frac{1}{M}}\sum _{i=1}^M\left(\theta -\widehat{{\theta }_i}\right).\hfill \end{array}$

We also calculate the coverage probability of the asymptotic confidence intervals for the parameter estimates to assess the reliability of these intervals in containing the true parameters. The coverage probability is defined as the proportion of intervals that include the true parameter value.

(41)$I_{ACI}\left(\lambda \right)=\left\{\begin{array}{c}1,\lambda \in \left[\widehat{\lambda }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}},\widehat{\lambda }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}}\right]\\ 0,\lambda \notin \left[\widehat{\lambda }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}},\widehat{\lambda }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{11}}}\right]\end{array}\right.,$

(42)$I_{ACI}\left(\theta \right)=\left\{\begin{array}{c}1,\theta \in \left[\widehat{\theta }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}},\widehat{\theta }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}}\right]\\ 0,\theta \notin \left[\widehat{\theta }-Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}},\widehat{\theta }+Z_{\frac{\gamma }{2}}\sqrt{\widehat{W_{22}}}\right]\end{array}\right..$

(43)$A{C}_{\theta }=E\left(I_{ACI}\left(\theta \right)\right),A{C}_{\lambda }=E\left(I_{ACI}\left(\lambda \right)\right).$

In subsequent simulations and experiments with real data, the significance level $\gamma$ of the confidence interval was set to 0.05.

We test the selection of initial values for the simulation experiments and find that the choice of initial values does not affect the convergence of these three estimation methods, which still provide good estimation results. For the fixed-point estimation method, the iterative initial values of the scale parameters are generated using an exponential distribution with an average of 1. For Newton–Raphson algorithm, the iterative initial values of shape parameters are generated using an exponential distribution with a mean of 1, and the iterative initial values of scale parameters are generated using an exponential distribution with a mean of 1.

During the estimation process, the nature of randomly generated middle-censored data may be poor, leading to a compromise in the convergence of the algorithm, resulting in a significant discrepancy between the simulated estimated parameters and the actual set parameters. This can affect the assessment of the accuracy of the estimation method. Therefore, for the parameter estimation results in the looped simulation process, the following outlier removal scheme is adopted: If a data point exceeds the range of $P_1-1.5\times IQR$ or $P_2+1.5\times IQR$, it is considered an outlier and is removed before assessing the accuracy of the estimation results. $P_1$ represents the first quartile of the estimated parameter results from all iterations, $P_2$ represents the third quartile of the estimated parameter results from all iterations, and $\mathrm{IQR}=P_2-P_1$. In the simulation results table, we use OR to denote Outlier Ratio, which means the proportion of simulation experiments where the estimation results are outliers out of the total number of experiments during the iterative simulation process.

By comparing the Newton–Raphson iterative method and the fixed-point estimation method, we obtain the accuracy and efficiency of the two methods under different parameter combinations. The results are shown in the following Table 1 and Table 2.

The figure below presents the line charts of the Newton–Raphson iteration method and the fixed-point estimation method in MPA estimation. With the legend given in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 present line graphs drawn based on Table 1 and Table 2. Through these line charts, we can more intuitively make horizontal comparisons of the linear loss, mean absolute error, and confidence interval coverage rate of the parameter estimation results between the two methods, as well as vertical comparisons of the changes in estimation accuracy under different sample sizes. This intuitive comparison allows us to comprehensively evaluate the strengths and weaknesses of the two methods.

Considering that parameter estimation in real scenarios requires not only the accuracy of computational results but also computational efficiency, we introduce a new metric: computation time per cycle. The computation time per cycle is obtained by dividing the total calculation time by the total number of simulation cycles H. This metric is used to compare the computational efficiency of different methods, and the corresponding formula is Formula (44).

(44)$\mathrm{Computation}\mathrm{time}\mathrm{per}\mathrm{cycle}={\displaystyle \frac{\mathrm{Total}\mathrm{cycle}\mathrm{time}}{H}}.$

The computational efficiency of the two estimation methods under different parameter combinations and data scales is shown in Table 3. Here, FP represents the computation time per cycle of the fixed-point estimation method, and NR represents the computation time per cycle of the Newton–Raphson iteration method.

In this simulation data estimation, after generating the data, we performed censoring on it. The endpoints of the censoring intervals follow exponential distributions with parameters 0.5 and 1, respectively. Censoring was applied to each data point; if the data fell within the censoring interval, it was replaced by the interval; otherwise, no censoring occurred. The number of estimation repetitions, H, was set to 1000. The initial shape and scale parameters for the first iteration followed exponential distributions with parameters 0.5 and 1, respectively.

Comparing the results of the MPA parameter estimation using the fixed-point estimation method and the Newton–Raphson iteration method, we can draw the following conclusions:

Computational Efficiency: The Newton–Raphson iteration method optimizes the two parameters by iterating the negative log-likelihood function, yielding the estimation results. The fixed-point estimation method, by incorporating the relationship between the two parameters of the Lomax distribution, simplifies the optimization process. In terms of computational efficiency, the fixed-point method is relatively more efficient.

Point Estimation Results:

(1). As the sample size increases, the linear and absolute errors of both estimation methods decrease correspondingly (applicable to both the shape and scale parameters). When the sample size reaches a certain level, the trend of error reduction slows down. Therefore, it can be seen that in practical applications, a moderate number of sample points can yield satisfactory estimation results. Sampling too many points increases the workload while providing limited improvement in the accuracy of the estimation results.

(2). Compared with the results of the three sets of parameters, the Newton–Raphson iteration method of the MPA estimation methods overall outperforms the fixed-point solution, especially in the estimation of the shape parameter $\theta$, where its advantage is more pronounced.

(3). The estimation of scale parameter $\lambda$ is not as good as that of shape parameter $\theta$.

Interval Estimation Results:

(1). The coverage rate of the confidence intervals for the two parameters in the Newton–Raphson iteration method increases with the increase of the initial sample size. The coverage rate of the confidence interval for the shape parameter is particularly satisfactory. However, for the fixed-point estimation method, there is a situation where the confidence interval coverage rate decreases as the amount of data increases, which is speculated to be due to the interval estimation itself not being suitable for application to the fixed-point estimation method.

(2). The Newton–Raphson iteration method of the MPA estimation methods generally yields better coverage rates for the confidence intervals compared to the fixed-point estimation method.

In addition, we also notice that the proportion of outliers decreases when the data volume is large. This indicates that increasing the data volume can enhance the convergence of the estimation results in the simulation experiments and reduce the occurrence of outliers.

Considering the comparison results of the two MPA estimation methods mentioned above, the method with better overall performance, namely the Newton–Raphson solution, is selected for further comparison with the EM algorithm in the following text.

5.1.2. Comparison of the EM Algorithm and MPA Estimation

In this simulation data estimation, the same data generation and censoring schemes are employed, and the same outlier removal method is applied to the results, thus obtaining the relevant results in Table 4 and Table 5.

In the EM algorithm experiment above, we use the exponential distribution with the mean of 2 to generate the iterative initial values of the shape parameters, and the exponential distribution with the mean of 1 to generate the iterative initial values of the scale parameters.

Similarly, for the convenience of intuitive comparison, with Figure 11 serves as the legend of Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17, the relevant line charts are drawn as Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17 present line graphs drawn based on tables.

From the tables, we can draw the following conclusions:

Computational Efficiency: As can be seen from Table 6, the Newton–Raphson solution for MPA estimation runs faster than the EM algorithm. This is because MPA does not require iteration over the samples, whereas the EM algorithm needs to perform sample iteration and integration operations.

Point Estimation Results:

(1). Both the maximum likelihood estimates obtained by the EM algorithm and the MPA estimates obtained by the Newton–Raphson method show significant bias when the sample size is small, but achieve higher accuracy when the sample size is large.

(2). In terms of estimating the shape parameter, the two estimation methods have similar estimation errors overall. However, when estimating the scale parameter, the EM algorithm provides more accurate results with smaller errors.

Interval Estimation Results:

(1). In terms of confidence interval coverage, the confidence interval coverage of shape parameter $\theta$ by both algorithms is inferior to the confidence interval coverage of shape parameter $\theta$ estimates. It is speculated that the nature of the scale parameter causes this phenomenon.

(2). As the amount of data increases, the confidence interval coverage of the estimated results of both methods increases.

Similarly, we can see from the results that the proportion of outliers decreases when the data volume is large.

5.1.3. Bayesian Estimation Results for the Shape Parameter

Before using Gibbs sampling for estimation, the convergence of the sampling results should first be verified. We selected the parameter combination $(\theta ,\lambda )=(0.7,0.7)$ and plotted the trace plots of the sampling results when the sample sizes were 1000 and 1500, respectively.

From the Figure 18 and Figure 19, we can observe that the sampling results have become stable, with relatively fixed ranges of fluctuation. Most of the values (black line) fluctuate around the initially set parameter results, which are represented by the red line. Moreover, the range of fluctuation becomes narrower when the sample size is larger, indicating that the parameters have converged. After conducting tests with other parameter combinations, we obtained the same results. Therefore, we can demonstrate that the Bayesian estimation of the Lomax shape parameter under middle-censoring using Gibbs sampling is both scientific and reliable.

After numerically verifying the convergence of the Gibbs sampling estimates, we also experiment with the choice of hyperparameters for the prior distribution. The results show that the impact of different hyperparameter choices on the estimation error is very subtle. The inferred reason is that during the sampling process, the hyperparameters are iterated continuously, and the results gradually converge around the correct estimates. The values influenced by the initial hyperparameters have been discarded and thus do not significantly affect the results. In the subsequent simulation experiments, we estimated the Bayesian estimation accuracy and confidence interval coverage of the shape parameter $\theta$ under different numbers of Gibbs sampling iterations.

In the simulation experiments, the number of periods N is set to 1000. From the trace plots of the sampling samples, we can see that the samples converge quickly. Therefore, we set the number of initial samples to be discarded Q to 100 to ensure the accuracy of the estimation results and the efficiency of the computational process. To compare the impact of the number of samples on the estimation results, we conducted experiments with 200, 500, and 1000 samples, respectively. The hyperparameters in the prior distribution are chosen as $a=2.5$ and $b=5$.

From Table 7, it can be concluded that the Bayesian estimation of the shape parameter, under the condition of a known scale parameter, achieves good fitting results. As the sample size increases, the accuracy significantly improves. When the number of initial samples discarded is 100, taking 500 sampling counts versus 1000 sampling counts does not significantly affect the accuracy of the results. Therefore, using an appropriate number of samples can yield good estimation results while saving computational power.

5.2. Real Data Analysis

When testing the estimation method for the following two sets of real data sets, we still used the same method of setting the initial iteration values in the simulation experiment and the same hyperparameter selection.

5.2.1. COVID-19 Datasets

According to the daily COVID-19 death counts divided by the number of new cases reported in [15], we obtain a set of real data samples that conform to the Lomax lifetime distribution.

Using the truncated interval of the exponential distribution with endpoints following parameters 10 and 30, we perform middle-censoring for each data point and obtain the data shown in Table 8.

Estimation results for complete data: $\theta =8.9602$, $\lambda =1.1554$. The estimated results of the truncated data are shown in Table 9 and Table 10, which are in line with theoretical expectations.

5.2.2. Global Extreme Weather and Climate Disaster Event Loss Datasets

The Global Extreme Weather and Climate Disaster Loss Dataset for 2000–2021 [29] covers data related to extreme weather and climate disaster events worldwide during the period from 2000 to 2021. We selected the number of fatalities caused by these disasters from 2017 to 2021 as the subject of our analysis. After detailed data processing and statistical analysis, we found that the data essentially conform to the Lomax distribution.

Using the truncated interval of the exponential distribution with endpoints following parameters 0.002 and 0.001, we perform middle-censoring for each data point and obtain the data shown in Table 11.

Estimation results for complete data: $\theta =2.9346$, $\lambda =653.5948$. The estimated results of the truncated data are shown in Table 12 and Table 13. In this practical application, the advantages of the EM algorithm are better demonstrated.

6. Summary and Future Prospects for Parameter Estimation of Lomax Distribution Under Middle-Censoring Data

In the field of Lomax parameter estimation with middle-censored data, this paper makes the following contributions:

In this paper, the details of EM algorithm for estimating two unknown parameters of Lomax distribution in middle-censored data are presented, and the results are verified by simulation experiments.

We propose two computational methods for parameter estimation after processing data using the Midpoint Approximation Method (MPA): the Newton–Raphson iteration method and the fixed-point method.

We calculate the ACIs for the two parameters in the middle-censored data, and take the covering probability of the confidence interval as one of the criteria to evaluate the estimation performance. The effectiveness of each estimation method is tested, and the applicability of the asymptotic confidence interval under intermediate censored data is also evaluated.

In the Bayesian estimation framework, we specifically estimate the shape parameter under middle-censoring using the Gamma prior distribution and the Gibbs sampling method, achieving satisfactory estimation results.

Through experiments with both simulated and real data, we validate the effectiveness and applicability of the proposed methods.

Based on the methods and achievements presented in this paper, the following future prospects are proposed, along with potential solutions and a ranking according to their importance:

(1). This paper only experiments with data where the middle-censoring endpoints follow an exponential distribution, without discussing the impact of different distributions of the censoring interval endpoints on estimation, such as the Weibull distribution, uniform distribution, etc. Most of the current research on middle-censoring focuses on endpoints that follow an exponential distribution. Not only for the Lomax distribution, but also for other distributions, the study of different distributions of the endpoints of middle-censoring intervals holds research value.

(2). In Bayesian estimation, we only discuss the estimation of the shape parameter $\theta$ when the scale parameter $\lambda$ is known. However, in practice, it is rare to encounter situations where only one parameter needs to be estimated. Therefore, the Bayesian estimation of both parameters is an important issue that needs to be addressed in the future. Additionally, our Bayesian estimation only considers the Gamma distribution as the prior distribution, without discussing other types of prior distributions, such as non-informative priors or Jeffreys priors.

(3). In simulation experiments, the estimation of confidence intervals is affected by the randomness of generating and censoring data, leading to situations where confidence intervals cannot be calculated. Therefore, the calculation of confidence interval coverage is based on the premise of generating middle-censored data for which progressive confidence intervals can be calculated. Another problem is that the asymptotic confidence interval coverage of scale parameters in simulation experiments is not high, which may be related to the nature of the parameters and the insufficient amount of data required for asymptotic confidence interval calculation. To address these interval estimation-related issues, Monte Carlo simulation can be considered to estimate confidence intervals in the future, or modified asymptotic confidence intervals can be adopted.

(4). Regarding the existence and uniqueness of the maximum likelihood estimation solution for the parameters of the Lomax distribution under middle-censored data, this paper does not provide theoretical proof due to the complexity of the formulas. Future research needs to continue exploring the parameter estimation problems for related distributions under middle-censored data to further improve research in this area.

(5). The EM algorithm has achieved good results in parameter estimation for the Lomax distribution under middle-censored data. However, during its operation, the issue of long computational time has arisen. In future research, the EM algorithm can be optimized through methods such as optimizing initial parameters, changing convergence criteria, and parallel computing. This will increase its computational speed while ensuring its accuracy.

Author Contributions

Conceptualization: P.R. and W.G.; Methodology: P.R. and W.G.; Software: P.R. and S.L.; Investigation: P.R. and S.L.; Writing—Original Draft: P.R.; Writing—Review and Editing: W.G.; Supervision: W.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are openly available in [15,29].

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Research about the Lomax lifetime distribution and middle-censoring.

Figure 2 Density curves of Lomax distribution with different parameters.

Figure 3 Example of middle-censored data.

Figure 4 The legend of the estimated result line chart in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.

Figure 5 Estimation results of the shape parameter $\theta$ under the parameter $(0.7,0.7)$. The legend is shown in Figure 4.

Figure 6 Estimation results of the scale parameter $\lambda$ under the parameter $(0.7,0.7)$. The legend is shown in Figure 4.

Figure 7 Estimation results of the shape parameter $\theta$ under the parameter $(0.6,0.9)$. The legend is shown in Figure 4.

Figure 8 Estimation results of the scale parameter $\lambda$ under the parameter $(0.6,0.9)$. The legend is shown in Figure 4.

Figure 9 Estimation results of the shape parameter $\theta$ under the parameter $(0.8,1.1)$. The legend is shown in Figure 4.

Figure 10 Estimation results of the scale parameter $\lambda$ under the parameter $(0.8,1.1)$. The legend is shown in Figure 4.

Figure 11 The legend of the estimated result line chart in Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17.

Figure 12 Estimation results of the shape parameter $\theta$ under the parameter $(0.7,0.7)$. The legend is shown in Figure 11.

Figure 13 Estimation results of the scale parameter $\lambda$ under the parameter $(0.7,0.7)$. The legend is shown in Figure 11.

Figure 14 Estimation results of the shape parameter $\theta$ under the parameter $(0.6,0.9)$. The legend is shown in Figure 11.

Figure 15 Estimation results of the scale parameter $\lambda$ under the parameter $(0.6,0.9)$. The legend is shown in Figure 11.

Figure 16 Estimation results of the shape parameter $\theta$ under the parameter $(0.8,1.1)$. The legend is shown in Figure 11.

Figure 17 Estimation results of the scale parameter $\lambda$ under the parameter $(0.8,1.1)$. The legend is shown in Figure 11.

Figure 18 Parameter $(\theta ,\lambda )=(0.7,0.7)$ data size $n=50$.

Figure 19 Parameter $(\theta ,\lambda )=(0.7,0.7)$ data size $n=200$.

Table 1

Comparison of shape parameter $\theta$ estimates between the fixed-point estimation and Newton–Raphson iteration methods.

$(\mathit{\theta },\mathit{\lambda })$	n	Fixed-Point Estimation				Newton–Raphson Iteration
		Bias	MAE	${\mathbf{AC}}_{\theta }$	OR	Bias	MAE	${\mathbf{AC}}_{\theta }$	OR
(0.7, 0.7)	30	0.2269	0.2520	0.722	0.060	0.1161	0.1922	0.847	0.048
	50	0.1435	0.1617	0.808	0.057	0.0458	0.1201	0.934	0.044
	100	0.0945	0.1111	0.855	0.037	−0.0022	0.0694	1.000	0.037
	200	0.0640	0.0793	0.883	0.021	−0.0300	0.0520	1.000	0.025
(0.6, 0.9)	30	0.2249	0.2935	0.726	0.076	0.1369	0.1713	0.835	0.090
	50	0.1530	0.2304	0.738	0.046	0.0955	0.1200	0.889	0.056
	100	0.0820	0.1481	0.792	0.039	0.0459	0.0647	0.962	0.053
	200	0.0558	0.1076	0.810	0.020	0.0234	0.0388	1.000	0.057
(0.8, 1.1)	30	0.2181	0.3698	0.672	0.083	0.0882	0.1978	0.802	0.072
	50	0.1741	0.3003	0.670	0.057	−0.0085	0.1285	0.912	0.061
	100	0.1114	0.2210	0.681	0.037	−0.0594	0.1022	0.976	0.029
	200	0.0744	0.1459	0.669	0.033	−0.0905	0.1007	1.000	0.016

Table 2

Comparison of scale parameter $\lambda$ estimates between the fixed-point estimation and Newton–Raphson iteration methods.

$(\mathit{\theta },\mathit{\lambda })$	n	Fixed-Point Estimation				Newton–Raphson Iteration
		Bias	MAE	${\mathbf{AC}}_{\lambda }$	OR	Bias	MAE	${\mathbf{AC}}_{\lambda }$	OR
(0.7, 0.7)	30	0.5895	0.6341	0.539	0.065	0.2003	0.4101	0.609	0.080
	50	0.3961	0.4285	0.594	0.050	0.1106	0.3002	0.694	0.040
	100	0.2666	0.2926	0.610	0.036	−0.0152	0.1708	0.826	0.033
	200	0.2037	0.2189	0.577	0.026	−0.0655	0.1252	0.876	0.015
(0.6, 0.9)	30	0.3997	0.4677	0.613	0.080	0.6912	0.8105	0.597	0.086
	50	0.2740	0.3412	0.617	0.063	0.5002	0.5937	0.612	0.059
	100	0.1974	0.2498	0.567	0.027	0.2670	0.3272	0.768	0.049
	200	0.1375	0.1704	0.567	0.028	0.1984	0.2361	0.806	0.043
(0.8, 1.1)	30	0.4165	0.5693	0.535	0.075	0.2823	0.6048	0.610	0.076
	50	0.3160	0.4418	0.492	0.059	0.0670	0.4076	0.631	0.054
	100	0.2361	0.3253	0.476	0.051	−0.0836	0.2831	0.646	0.030
	200	0.1899	0.2386	0.398	0.037	−0.1890	0.2345	0.665	0.037

Table 3

Computation time per cycle of fixed-point estimation method and Newton–Raphson iteration method.

$(\mathit{\theta },\mathit{\lambda })$	(0.7, 0.7)		(0.6, 0.9)		(0.8, 1.1)
$\mathit{n}$	FP	NR	FP	NR	FP	NR
30	0.00143	0.00259	0.00264	0.00266	0.00295	0.00235
50	0.00166	0.00374	0.00168	0.00365	0.00200	0.00267
100	0.00190	0.00508	0.00183	0.00563	0.00248	0.00396
200	0.00247	0.00748	0.00272	0.00935	0.00277	0.00713

Table 4

Comparison of shape parameter $\theta$ estimates between the EM Algorithm and Newton–Raphson iteration methods.

$(\mathit{\theta },\mathit{\lambda })$	n	EM Algorithm				Newton–Raphson Iteration
		Bias	MAE	${\mathbf{AC}}_{\theta }$	OR	Bias	MAE	${\mathbf{AC}}_{\theta }$	OR
(0.7, 0.7)	30	0.1764	0.2043	0.809	0.071	0.1161	0.1922	0.847	0.048
	50	0.1318	0.1570	0.848	0.046	0.0458	0.1201	0.934	0.044
	100	0.0703	0.0910	0.926	0.032	−0.0022	0.0694	1.000	0.037
	200	0.0363	0.0609	0.952	0.022	−0.0300	0.0520	1.000	0.025
(0.6, 0.9)	30	0.1953	0.2047	0.799	0.074	0.1369	0.1713	0.835	0.090
	50	0.1369	0.1459	0.839	0.044	0.0955	0.1200	0.889	0.056
	100	0.0901	0.0960	0.921	0.038	0.0459	0.0647	0.962	0.053
	200	0.0610	0.0642	0.969	0.026	0.0234	0.0388	1.000	0.057
(0.8, 1.1)	30	0.2008	0.2572	0.707	0.074	0.0882	0.1978	0.802	0.072
	50	0.1041	0.1713	0.778	0.058	−0.0085	0.1285	0.912	0.061
	100	0.0503	0.1129	0.830	0.037	−0.0594	0.1022	0.976	0.029
	200	0.0204	0.0775	0.874	0.025	−0.0905	0.1007	1.000	0.016

Table 5

Comparison of scale parameter $\lambda$ estimates between the EM Algorithm and Newton–Raphson iteration methods.

$(\mathit{\theta },\mathit{\lambda })$	n	EM Algorithm				Newton–Raphson Iteration
		Bias	MAE	${\mathbf{AC}}_{\lambda }$	OR	Bias	MAE	${\mathbf{AC}}_{\lambda }$	OR
(0.7, 0.7)	30	0.3543	0.4383	0.662	0.077	0.2003	0.4101	0.609	0.080
	50	0.2572	0.3274	0.662	0.049	0.1106	0.3002	0.694	0.040
	100	0.1261	0.1852	0.759	0.034	−0.0152	0.1708	0.826	0.033
	200	0.0794	0.1340	0.785	0.013	−0.0655	0.1252	0.876	0.015
(0.6, 0.9)	30	0.6013	0.6691	0.608	0.070	0.6912	0.8105	0.597	0.086
	50	0.4150	0.4776	0.642	0.057	0.5002	0.5937	0.612	0.059
	100	0.2701	0.3090	0.701	0.029	0.2670	0.3272	0.768	0.049
	200	0.1723	0.2028	0.761	0.027	0.1984	0.2361	0.806	0.043
(0.8, 1.1)	30	0.5249	0.7043	0.542	0.079	0.2823	0.6048	0.610	0.076
	50	0.2648	0.4598	0.607	0.060	0.0670	0.4076	0.631	0.054
	100	0.1244	0.3020	0.610	0.038	−0.0836	0.2831	0.646	0.030
	200	0.0637	0.2073	0.630	0.020	−0.1890	0.2345	0.665	0.037

Table 6

Computation time per cycle of EM algorithm and MPA estimation method.

$(\mathit{\theta },\mathit{\lambda })$	(0.7, 0.7)		(0.6, 0.9)		(0.8, 1.1)
$\mathit{n}$	EM	MPA	EM	MPA	EM	MPA
30	0.05542	0.00259	0.04147	0.00266	0.09993	0.00235
50	0.06408	0.00374	0.05607	0.00365	0.08504	0.00267
100	0.12139	0.00508	0.10925	0.00563	0.15334	0.00396
200	0.20563	0.00748	0.17732	0.00935	0.25766	0.00713

Table 7

Performance metrics for different settings.

$(\mathit{\theta },\mathit{\lambda })$	n	M = 200			M = 500			M = 1000
		Bias	MAE	${\mathbf{AC}}_{\theta }$	Bias	MAE	${\mathbf{AC}}_{\theta }$	Bias	MAE	${\mathbf{AC}}_{\theta }$
(0.7, 0.7)	30	0.0420	0.0860	0.996	0.0401	0.0907	0.998	0.0372	0.0855	0.994
	50	0.0259	0.0699	1.000	0.0344	0.0708	0.998	0.0269	0.0697	0.998
	100	0.0155	0.0503	0.999	0.0129	0.0489	0.998	0.0105	0.0496	1.000
	200	0.0078	0.0364	1.000	0.0072	0.0349	1.000	0.0058	0.0355	1.000
(0.6, 0.9)	30	0.0754	0.0918	0.999	0.0659	0.0901	0.999	0.0690	0.0888	0.996
	50	0.0574	0.0719	0.998	0.0558	0.0702	0.999	0.0560	0.0724	0.996
	100	0.0446	0.0548	0.999	0.0375	0.0495	0.999	0.0417	0.0527	1.000
	200	0.0290	0.0362	1.000	0.0303	0.0361	1.000	0.0283	0.0356	1.000
(0.8, 1.1)	30	−0.0050	0.1029	0.993	−0.0035	0.0954	0.997	0.0009	0.0965	0.993
	50	0.0085	0.0845	0.993	−0.0082	0.0808	0.998	0.0021	0.0770	0.998
	100	−0.0072	0.0585	0.996	−0.0079	0.0603	0.998	−0.0063	0.0606	0.996
	200	−0.0054	0.0441	0.999	−0.0045	0.0418	0.998	−0.0054	0.0427	0.998

Table 8

COVID-19 medical data after censoring treatment.

	0.0149	0.0096	0.0322	0.0558	0.0636	0.0634	0.0467	0.0541
	0.0610	0.0282	0.1021	0.0859	0.1266	0.1718	0.1896	0.1052
	0.1886	0.1052	0.1394	0.1653	0.1865	0.1513	0.1607	0.1019
${\delta }_i=0$	0.1426	0.0835	1.0842	0.1117	0.1192	0.0726	0.0723	0.1573
$n_1=70$	0.0683	0.1842	0.1524	0.1242	0.1079	0.0590	0.1293	0.1068
	0.0745	0.2047	0.1998	0.1340	0.3541	0.1960	0.1919	0.0975
	0.2443	0.1311	0.0238	0.0393	0.0352	0.0157	0.3040	0.2362
	0.0717	0.0558	0.1155	0.2163	0.1155	0.2163	0.0978	0.0888
	0.0953	0.0491	0.2164	0.1310	0.1333	0.1124
	(0.0190, 0.1672)		(0.0035, 0.0657)		(0.0193, 0.4777)		(0.0403, 0.3831)
	(0.0828, 0.6752)		(0.0715, 0.8918)		(0.0029, 0.3830)		(0.1430, 0.5425)
	(0.0803, 0.7193)		(0.0105, 0.1338)		(0.0315, 0.2710)		(0.0057, 0.2977)
	(0.0296, 0.1984)		(0.0698, 0.2709)		(0.0391, 0.6342)		(0.0577, 0.2009)
${\delta }_i=1$	(0.0965, 0.3380)		(0.0480, 0.0703)		(0.0692, 0.2260)		(0.0299, 0.1601)
$n_2=37$	(0.0206, 0.4020)		(0.0573, 0.1095)		(0.0243, 0.2168)		(0.0061, 0.3026)
	(0.0314, 0.1729)		(0.0419, 0.2356)		(0.1373, 0.4450)		(0.0262, 0.2354)
	(0.0949, 0.2687)		(0.0567, 0.5827)		(0.0516, 0.1286)		(0.0099, 0.3902)
	(0.0424, 0.2731)		(0.0486, 0.4021)		(0.0976, 0.2715)		(0.0375, 0.3590)
	(0.0385, 0.3161)

Table 9

Estimation results of double parameters for COVID-19 real data.

	Estimation Results $(\mathit{\theta },\mathit{\lambda })$	Bias
EM	(8.4261, 1.1237)	(−0.5340, −0.0317)
MPA	(9.5718, 1.2946)	(0.6115, 0.1391)

Table 10

Estimation results of shape parameters for COVID-19 real data.

	Estimation Results $\mathit{\theta }$	Bias
M = 200	9.7846	0.8244
M = 500	9.6675	0.7073
M = 1000	9.6198	0.6596

Table 11

Global extreme weather and climate disaster event loss data after censoring treatment.

	78	36	3059	167	184	264	284	514
	253	251	88	97	142	212	112	53
	182	504	246	84	136	17	143	119
	148	53	59	179	300	206	72	5
${\delta }_i=0$	628	128	400	188	356	1900	115	99
$n_1=75$	567	868	3	162	257	212	1460	22
	107	1922	82	285	410	1924	2556	110
	400	111	97	21	197	0	116	260
	3	352	1282	201	23	226	87	13
	0	0	457
	(132.98, 1068.8)		(290.67, 307.32)		(137.82, 600.24)		(838.38, 1264.8)
	(80.467, 186.90)		(120.93, 640.36)		(2.3398, 1565.2)		(155.36, 907.70)
${\delta }_i=1$	(58.948, 620.32)		(62.411, 803.27)		(295.84, 941.48)		(244.63, 382.45)
$n_2=25$	(67.397, 191.28)		(57.946, 1565.5)		(163.48, 516.84)		(75.673, 241.76)
	(118.45, 165.72)		(21.115, 298.93)		(74.457, 212.94)		(47.828, 588.44)
	(696.11, 1202.0)		(184.87, 603.93)		(209.73, 255.58)		(399.54, 752.28)
	(37.435, 486.19)

Table 12

Estimation results of double parameters for global extreme weather and climate disaster event loss data.

	Estimation Results $(\mathit{\theta },\mathit{\lambda })$	Bias
EM	(2.9602, 688.3042)	(0.0256, 34.7454)
MPA	(2.0521, 394.2962)	(−0.8825, −259.8826)

Table 13

Estimation results of shape parameters for global extreme weather and climate disaster event loss data.

	Estimation Results $\mathit{\theta }$	Bias
M = 200	2.8004	−0.1342
M = 500	2.8081	−0.1265
M = 1000	2.7967	−0.1379

Appendix A

Theorem A1. 

The Gamma distribution is a conjugate prior distribution of the shape parameters of the Lomax life distribution.

Proof. 

Suppose the sample data is denoted as data, consisting of n sample data points, which are $t_1,t_2,\dots ,t_n$, respectively.

The prior distribution of shape parameters is the Gamma distribution of parameters a and b, which is shown as follows:(A1)$$\pi \left(\theta \right)={\displaystyle \frac{b^a}{\mathsf{\Gamma }\left(a\right)}}{\theta }^{a-1}{e}^{-b\theta }.$$

Suppose the experimental data as data; then, the likelihood function of data can be represented as $p(\mathrm{data}\mid \theta )$, and the joint distribution of data and $\theta$ is $h(\mathrm{data},\theta )$. By integrating $h(\mathrm{data},\theta )$ over the parameter space $\Theta$, we can obtain the marginal density of data, denoted as $m\left(\mathrm{data}\right)$.(A2)$$p\left(data\right|\theta )={\displaystyle \frac{{\theta }^n}{{\lambda }^n}}\prod _{i=1}^n{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\left(\theta +1\right)}$$(A3)$$\begin{array}{c}\hfill h(data,\theta )=p\left(data\right|\theta )\pi \left(\theta \right)={\displaystyle \frac{{\theta }^n}{{\lambda }^n}}\prod _{i=1}^n{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\left(\theta +1\right)}{\displaystyle \frac{b^a}{\mathsf{\Gamma }\left(a\right)}}{\theta }^{a-1}{e}^{-b\theta }.\end{array}$$

The posterior distribution of the shape parameters obtained by computational simplification is as follows:(A4)$$\pi \left(\theta |data\right)\propto {\theta }^n\prod _{i=1}^n{\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)}^{-\theta }{\theta }^{a-1}{e}^{-b\theta }.$$

Thus, the conditional posterior distribution of the shape parameter is(A5)$$Gamma\left(n+a,b+\sum _{i=1}^{n}ln\left(1+{\displaystyle \frac{t_i}{\lambda }}\right)\right).$$

To sum up, we prove that the posterior distribution obtained when the Gamma distribution is selected as the prior distribution of the shape parameter is also the Gamma distribution. Therefore, the Gamma distribution is a conjugate prior distribution of the Lomax shape parameters.

□

References

1. Leiderman, P.H.; Babu, B.; Kagia, J.; Kraemer, H.C.; Leiderman, G.F. African Infant Precocity and Some Social Influences during the First Year. Nature; 1973; 242, pp. 247-249. [DOI: https://dx.doi.org/10.1038/242247a0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4716589]

2. Groeneboom, P.; Wellner, J.A. Information Bounds and Nonparametric Maximum Likelihood Estimation; Birkhäuser: Basel, Switzerland, 1992; Volume 19.

3. Schick, A.; Yu, Q. Consistency of the GMLE with Mixed Case Interval-Censored Data. Scand. J. Stat.; 2000; 27, pp. 45-55. [DOI: https://dx.doi.org/10.1111/1467-9469.00177]

4. Huang, J. Asymptotic properties of nonparametric estimation based on partly interval-censored data. Stat. Sin.; 1999; 9, pp. 501-520.

5. Jammalamadaka, S.R.; Mangalam, V. Nonparametric estimation for middle-censored data. J. Nonparametr. Stat.; 2003; 15, pp. 253-265. [DOI: https://dx.doi.org/10.1080/1048525031000089266]

6. Jammalamadaka, S.R.; Iyer, S.K. Approximate self-consistency for middle-censored data. J. Stat. Plan. Inference; 2004; 124, pp. 75-86. [DOI: https://dx.doi.org/10.1016/S0378-3758(03)00199-X]

7. Yu, Q.Q.; Wong, G.Y.C.; Li, L.X. Asymptotic properties of self-consistent estimators with mixed interval-censored data. Ann. Inst. Stat. Math.; 2001; 53, pp. 469-486. [DOI: https://dx.doi.org/10.1023/A:1014656726982]

8. Iyer, S.K.; Jammalamadaka, S.R.; Kundu, D. Analysis of middle-censored data with exponential lifetime distributions. J. Stat. Plan. Inference; 2008; 138, pp. 3550-3560. [DOI: https://dx.doi.org/10.1016/j.jspi.2007.03.062]

9. Bennett, N.A. Some Contributions to Middle-Censoring. Ph.D. Thesis; University of California: Santa Barbara, CA, USA, 2011; 3473719.Available online: https://www.proquest.com/dissertations-theses/some-contributions-middle-censoring/docview/896352567/se-2 (accessed on 1 June 2011).

10. Yan, W.; Yimin, S.; Min, W. Statistical inference for dependence competing risks model under middle censoring. J. Syst. Eng. Electron.; 2019; 30, pp. 209-222. [DOI: https://dx.doi.org/10.21629/JSEE.2019.01.20]

11. Jammalamadaka, S.R.; Leong, E. Analysis of discrete lifetime data under middle-censoring and in the presence of covariates. J. Appl. Stat.; 2015; 42, pp. 905-913. [DOI: https://dx.doi.org/10.1080/02664763.2014.993364]

12. Lomax, K.S. Business failure: Another example of the analysis of the failure data. J. Am. Stat. Assoc.; 1954; 49, pp. 847-852. [DOI: https://dx.doi.org/10.1080/01621459.1954.10501239]

13. Lun, Z.; Khattree, R. An overlooked Lomax-exponential connection. Commun. Stat.–Theory Methods; 2022; 51, pp. 26-28. [DOI: https://dx.doi.org/10.1080/03610926.2020.1734836]

14. Aljohani, H.M. Estimation for the P(X > Y) of Lomax distribution under accelerated life tests. Heliyon; 2024; 10, e25802. [DOI: https://dx.doi.org/10.1016/j.heliyon.2024.e25802]

15. Alsuhabi, H.; Alkhairy, I.; Almetwally, E.M.; Almongy, H.M.; Gemeay, A.M.; Hafez, E.H.; Aldallal, R.A.; Sabry, M. A superior extension for the Lomax distribution with application to COVID-19 infections real data. Alex. Eng. J.; 2022; 61, pp. 11077-11090. [DOI: https://dx.doi.org/10.1016/j.aej.2022.03.067]

16. Cramer, E.; Schmiedt, A.B. Progressively Type-II censored competing risks data from Lomax distributions. Comput. Stat. Data Anal.; 2011; 55, pp. 1285-1303. [DOI: https://dx.doi.org/10.1016/j.csda.2010.09.017]

17. Alfaer, N.M.; Aljohani, H.M. Balanced joint progressively hybrid Type-I censoring samples in estimating the lifetime Lomax distributions. Complexity; 2021; 2021, 9929691. [DOI: https://dx.doi.org/10.1155/2021/9929691]

18. El-Sherpieny, E.S.A.; Almetwally, E.M.; Muhammed, H.Z. Progressive Type-II hybrid censored schemes based on maximum product spacing with application to Power Lomax distribution. Phys. A Stat. Mech. Its Appl.; 2020; 553, 124251. [DOI: https://dx.doi.org/10.1016/j.physa.2020.124251]

19. Hasaballah, M.M.; Balogun, O.S.; Bakr, M.E. Non-Bayesian and Bayesian estimation for Lomax distribution under randomly censored with application. AIP Adv.; 2024; 14, 025318. [DOI: https://dx.doi.org/10.1063/5.0176326]

20. Davarzani, N.; Parsian, A. Statistical inference for discrete middle-censored data. J. Stat. Plan. Inference; 2011; 141, pp. 1455-1462. [DOI: https://dx.doi.org/10.1016/j.jspi.2010.10.012]

21. Giles, D.E.A.; Feng, H.; Godwin, R.T. On the Bias of the Maximum Likelihood Estimator for the Two-Parameter Lomax Distribution. Commun.-Stat.–Theory Methods; 2013; 42, pp. 1934-1950. [DOI: https://dx.doi.org/10.1080/03610926.2011.600506]

22. Dempster, A.P.; Laird, N.M.; Rubin, D.B. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. Ser. B (Methodol.); 1977; 39, pp. 1-22. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1977.tb01600.x]

23. Wu, C.F.J. On the convergence properties of the EM algorithm. Ann. Stat.; 1983; 11, pp. 95-103. [DOI: https://dx.doi.org/10.1214/aos/1176346060]

24. Asl, M.N.; Belaghi, R.A.; Bevrani, H. Classical and Bayesian inferential approaches using Lomax model under progressively Type-I hybrid censoring. J. Comput. Appl. Math.; 2018; 343, pp. 397-412. [DOI: https://dx.doi.org/10.1016/j.cam.2018.04.028]

25. Turkson, A.; Ayiah-Mensah, F.; Nimoh, V. Handling censoring and censored data in survival analysis: A standalone systematic literature review. Int. J. Math. Math. Sci.; 2021; 2021, 9307475. [DOI: https://dx.doi.org/10.1155/2021/9307475]

26. Buzaridah, M.M.; Ramadan, D.A.; El-Desouky, B.S. Estimation of some lifetime parameters of flexible reduced logarithmic-inverse Lomax distribution under progressive Type-II censored data. J. Math.; 2022; 2022, 1690458. [DOI: https://dx.doi.org/10.1155/2022/1690458]

27. Okasha, H.; Lio, Y.; Albassam, M. On reliability estimation of Lomax distribution under adaptive Type-I progressive hybrid censoring scheme. Mathematics; 2021; 9, 2903. [DOI: https://dx.doi.org/10.3390/math9222903]

28. Kundu, D. Bayesian inference and life testing plan for the Weibull distribution in presence of progressive censoring. Technometrics; 2008; 50, pp. 144-154. [DOI: https://dx.doi.org/10.1198/004017008000000217]

29. Yang, S.N. Global Extreme Weather and Climate Disaster Loss Dataset (2000–2021). 2023; Available online: https://data.casearth.cn/dataset/65388536819aec0f26f4c93a (accessed on 27 May 2023).