On asymptotics of solutions for superdiffusion and subdiffusion equations with the Riemann-Liouville fractional derivative

1.   Introduction

Lifetime models have found widespread application in statistical modeling across various scientific and engineering domains. Lindley distribution, as one of the classical distributions, was first proposed by Lindley ^[1]. Lindley distribution is highly flexible and has a wide range of applications, such as in the field of medicine, astrophysics, and reliability engineering. However, it has strong limitations in processing complex data, such as skewed and multi-peak data. Based on this, many researchers improved the original Lindley distribution by adding parameters. With the help of the power exponentiated family of distributions, Rajitha and Akhilnath ^[2] added two parameters to the original Lindley distribution. They called it PEL distribution and notes its higher flexiblility compared to the original model. Fatehi and Chhaya ^[3] extended the Lindley distribution into the extended odd Weibull-Lindley family. Ashour and Eltehiwy ^[4] introduced a three-parameter exponentiated power Lindley (EPL) distribution by extending the two-parameter power Lindley distribution. Alizadeh et al. ^[5] defined a four-parameter exponentiated power Lindley power series distribution on the EPL distribution and found that the newly proposed model provided a better fit than the original Lindley distribution to real datasets.

In most real-life testing and reliability experiments, it is seldom possible to wait until all test samples fail; in other words, it is difficult for investigators to observe the lifetime of all items under test, so experimental data obtained often contain censored data. The development and replacement of censored samples have been the focus of many researchers ^[6,7,8]. Type-Ⅰ and type-Ⅱ censored schemes, as classical methods in right examination, can only move the unit point at the end of the experiment, which lacks some flexibility ^[9]. The progressively type-Ⅱ censored scheme is popular for its flexibility, whose various properties and applications have been extensively studied. Balakrishnan et al. ^[10] discussed the maximum likelihood estimation and the corresponding interval estimation of extreme value distributions under progressively type-Ⅱ censored samples. Based on progressively type-Ⅱ censored samples, Seo et al. ^[11] studied the hierarchical Bayesian estimation of the unknown parameters of a lifetime distribution with a bathtub-shaped failure rate function. Alshenawy et al. ^[12] used the maximum likelihood estimation method and the maximum product spacing method to estimate the parameters of the extended odd Weibull exponential distribution under progressively type-Ⅱ censored samples. Their study further delved into the construction of both asymptotic and bootstrap confidence intervals for the said parameters.

This paper aims to introduce a variant of the Lindley distribution, referred to as the extension of the generalized Lindley (NGL) distribution. This extension is developed under progressively type-Ⅱ censored samples with the objective of broadening the applicability of the traditional Lindley distribution. The NGL distribution is highly flexible, featuring many variants of the Lindley distribution and exponential function. Another purpose of this paper is to evaluate the estimator performance of the NGL distribution in preparation for the subsequent processing of real datasets.

The remainder of this paper is organized as follows: In Section 2, we introduce the NGL distribution and its basic properties. The numerical characteristics of the proposed distribution are investigated in Section 3. We study three estimators of this distribution in Section 4, obtain the corresponding point and interval estimates, and perform Monte Carlo simulations in Section 5. In Section 6, we examine the practical application of the proposed distribution using a real dataset. In Section 7, the findings of this paper are summarized, and future research priorities are indicated.

2.   Generalized Lindley distribution

In this section, we will introduce the NGL distribution, whose probability density function (PDF) and cumulative distribution function (CDF) are

Here, $\lambda$ and $\theta$ are the scale and shape parameters of the $NGL(\lambda, \theta)$, respectively.

● If $\theta = 1/(\lambda +1)$, then the CDF is the Lindley distribution, i.e.,

● If $\theta = \beta /(\lambda +\beta)$, then the CDF is a two-parameter Lindley distribution ^[13], i.e.,

● If $\theta \to {0}^{+}$, then the CDF is exponential distribution, i.e.,

● If $\theta = 1/(\eta +1)$, then the CDF is a two-parameter Lindley distribution ^[14], i.e.,

The survival function (SF) and hazard rate function (HRF) of the NGL distribution are:

There are two main models for the NGL distribution: the inverted J-type and the unimodal and left-leaning (see Figure 1). This type of model is very suitable for the study of product life and species abundance distribution. $\theta$, the shape parameter, has a large influence on the probability density distribution: when $\theta$ increases, the PDF image changes and shows a constant tendency to shift to the right. $\lambda$ also has a certain impact on the NGL distribution; with an increase in $\lambda$, the peak value of PDF increases, and the image presents a left-leaning trend.

The HRF curves also have three shapes: increasing, constant, and upside-down bathtub (see Figure 2). With an increase in $\lambda$ and $\theta$, the slope of (0, 2) will increase continuously along with the decrease of peak value.

3.   Statistical properties of NGL distribution

In this section, we discuss some important statistical properties of NGL distribution, such as moments, kurtosis and skewness, quantile functions, order statistical functions, and other statistical properties.

3.1.   Moments and related measures for the NGL distribution

As one of the most important digital features, the $r$-$th$ moment plays an important role in both application and theory. It can be represented by the following formula

The first four moments can be easily calculated:

Also, the variance of $X$ is

3.2.   Coefficient of skewness and kurtosis for the NGL distribution

The coefficient of kurtosis and skewness are important for describing the tail shape, peak degree, and asymmetry of probability distributions ^[15]. Let NCS stand for coefficient of skewness and NCK for coefficient of kurtosis. According to Eqs (10) and (11), NCS and NCK can be obtained as:

The coefficient of kurtosis and skewness of the NGL distribution show a certain regularity. In the NGL distribution, the coefficient of kurtosis and skewness are usually negatively correlated with $\theta$ and positively correlated with $\lambda$. This statistical property suggests that the thickness and asymmetry of the NGL distribution tail tend to decrease as $\theta$ increases, while an increase in $\lambda$ increases these characteristics. Table 1 and Figure 3 illustrate different values of the coefficient of kurtosis and skewness.

3.3.   Moment-generating function for the NGL distribution

Using Eq (9), the moment-generating function of the NGL distribution can be obtained as

3.4.   Quantile function for the NGL distribution

The quantile function can be obtained from the inverse function of CDF. That is, $x = {F}^{-1}\left(R\right)$. By applying Eq (2), we have

Let $t = (1-R){e}^{\lambda x}$, then $x = \frac{1}{\lambda }\mathit{ln}\frac{t}{1-R}$. Hence

Then

Here, $W$ is Lambert $W$ function. According to the basic properties of the Lambert $W$ function, the quantile function of the NGL distribution can be obtained as

3.5.   Order statistics of the NGL distribution

Order statistics are an important analytical tool for identifying outliers. The earliest failure time of a product can be expressed by the minimum order statistic (the smallest observed value in the sample), and the longest life of a product can be estimated by the maximum order statistic (the largest observed value in the sample). Let $({X}_{1}, {X}_{2}, ..., {X}_{n})$ be a random sample from the NGL distribution, which is sorted from smallest to largest as $({X}_{\left(1\right)}, {X}_{\left(2\right)}, ..., {X}_{\left(n\right)})$, where the $i$-$th$ order statistic is ${X}_{\left(i\right)}$. The PDF of the $i$-$th$ order statistic ${X}_{\left(i\right)}$ is ${f}_{{X}_{\left(i\right)}}\left(x\right) = \frac{n!}{(i-1)!(n-i)!}F(x)^{i-1}[1-F\left(x\right)]^{n-i}f\left(x\right)$. By inserting Eqs (1) and (2), the PDF of $i$-$th$ order statistic ${X}_{\left(i\right)}$ of the NGL distribution is

The PDF of the minimum order statistics ${X}_{\left(1\right)}$ of the NGL distribution can be obtained as

The PDF of the maximum order statistics ${X}_{\left(n\right)}$ of the NGL distribution can be obtained as

4.   Parameter estimation of the NGL distribution

4.1.   Maximum product spacing estimation

Maximum product spacing (MPS) estimation is a robust method. It uses an optimization algorithm to find the corresponding parameter values that maximize the product spacing of the parametric functions ^[16]. Let $({X}_{1:m:n}, {X}_{2:m:n}, \cdots, {X}_{m:m:n})$ be the progressively type-Ⅱ censored sample of the NGL distribution, and ${x}_{1:m:n}, {x}_{2:m:n}, \cdots, {x}_{m:m:n}$ is the observation of the sample $({X}_{1:m:n}, {X}_{2:m:n}, \cdots, {X}_{m:m:n})$. For convenience, in the following discussion, we always set ${x}_{i} = {x}_{i:m:n}$. The product spacing function under the progressively type-Ⅱ censored scheme is ^[17]:

Using Eq (2), the product spacing function of the NGL distribution can be obtained as

and the log-product spacing function is given by

where

MPS estimators of $\lambda$ and $\theta$, denoted by ${\widehat{\lambda }}_{MPS}$ and ${\widehat{\theta }}_{MPS}$, are the solutions of Eqs (25) and (26).

Here,

${\widehat{\lambda }}_{MPS}$ and ${\widehat{\theta }}_{MPS}$ cannot be obtained directly from the above equations, so we used the Newton Raphson algorithm to obtain the approximate maximum product spacing estimates of these parameters. The steps are as follows:

(ⅰ) Establish the corresponding log-product spacing function.

(ⅱ) Compute the gradient vector and Hessian matrix for the log-product spacing function:

(ⅲ) Select the appropriate initial value $({\lambda }^{\left(0\right)}, {\theta }^{\left(0\right)}{)}^{T}.$

(ⅳ) The parameters are updated by Newton iteration formula:

(ⅴ) Repeat (ⅳ) until the change in parameter values between the two iterations is less than a very small threshold: $|{\lambda }^{(k+1)}-{\lambda }^{\left(k\right)}| < \epsilon, |{\theta }^{(k+1)}-{\theta }^{\left(k\right)}| < \epsilon$

(ⅵ) The final estimated parameters can be obtained and are noted by ${\widehat{\lambda }}_{MPS} = {\lambda }^{(k+1)}, {\widehat{\theta }}_{MPS} = {\theta }^{(k+1)}.$

4.2.   Maximum likelihood estimation

Under the progressively type-Ⅱ censored sample, the likelihood function is (^[18])

Here, $c = n(n-{P}_{1}-1)(n-{P}_{1}-{P}_{2}-2)\cdots (n-{\sum }_{i = 1}^{m-1}({P}_{i}+1))$. Using Eqs (1) and (7), the likelihood function of the NGL distribution can be represented as

The log-likelihood function of the NGL distribution is

The ML estimator of $\lambda$ and $\theta$ can be obtained by solving the following equations:

Newton Raphson algorithm to obtain the corresponding result. Due to the complexity of Eqs (35) and (36), it is difficult to judge the existence and uniqueness of ${\widehat{\lambda }}_{ML}$ and ${\widehat{\theta }}_{ML}$ by conventional numerical methods. In this paper, a graphical tool, the counter diagram, is employed, as referenced in Alotaibi et al. ^[19]. Setting the real values $(\lambda, \theta) = \left(\mathrm{0.5, 0.75}\right)$, $(m, n) = \left(\mathrm{50,100}\right)$, $P = \left(0*\mathrm{49, 50}\right)$ and generating progressively type-Ⅱ censored samples, the counter diagram of ${\widehat{\lambda }}_{ML}$ and ${\widehat{\theta }}_{ML}$, which is obtained from Eq (34), is shown in Figure 4. It shows that the likelihood function obtained has only one obvious peak, that is, there is uniqueness.

Figure 4 shows the existence and uniqueness of ML estimates of $\lambda$ and $\theta$, with $({\widehat{\lambda }}_{ML}, {\widehat{\theta }}_{ML}) = \left(\mathrm{0.855, 0.625}\right)$. The ML estimation of $S\left(t\right)$ and $h\left(t\right)$ can be obtained by substituting ${\widehat{\lambda }}_{ML}$ and ${\widehat{\theta }}_{ML}$ ^[20]:

4.3.   Asymptotic confidence interval

In the above analysis, we have obtained point estimates for $\lambda$, $\theta$, $S\left(t\right)$, and $h\left(t\right)$. Point estimates provide a singular numerical approximation for unknown parameters, yet they fall short of encapsulating the full spectrum of uncertainty associated with these parameters. At this point, we need to transition from point estimation to interval estimation to get a more comprehensive understanding. The ML estimator is asymptotically normal ^[21]. That is, $\sqrt{n}(\widehat{\Theta }-\Theta)\stackrel{d}{\to }N(0, {I}^{-1}(\Theta \left)\right)$, where $\Theta = (\lambda, \theta)$, and ${I}^{-1}\left(\Theta \right)$ as the inverse of the information matrix for unknown parameters, and it can be obtained as follows:

with the following elements

The $100(1-\alpha)\%$ asymptotic confidence intervals (ACIs) for $\lambda$ and $\theta$ are $\left(\widehat{\lambda }\pm {z}_{\alpha /2}\sqrt{Var\left(\widehat{\lambda }\right)}\right)$ and $\left(\widehat{\theta }\pm {z}_{\alpha /2}\sqrt{Var\left(\widehat{\theta }\right)}\right),$ respectively. Here, ${z}_{\alpha /2}$ represents the upper $\alpha /2$ percentage point of $N\left(\mathrm{0, 1}\right)$.

Similarly, the ACIs of $S\left(t\right)$ and $h\left(t\right)$ can also be constructed by calculating the corresponding variances, among which one of the most famous methods is Greene's delta method ^[22]. Under this method, the distribution of $\widehat{S}\left(t\right)$ ($\widehat{h}\left(t\right)$) is approximately a normal distribution with mean $S\left(t\right)$ ($h\left(t\right)$) and variance ${\sigma }_{S}^{2} = {\varDelta }_{S}{I}^{-1}(\lambda, \theta){\varDelta }_{S}^{T}{|}_{(\lambda = \widehat{\lambda }, \theta = \widehat{\theta })}({\sigma }_{h}^{2} = {\varDelta }_{h}{I}^{-1}(\lambda, \theta \left){\varDelta }_{h}^{T}{|}_{(\lambda = \widehat{\lambda }, \theta = \widehat{\theta })}\right)$ ^[23], where ${\varDelta }_{S} = \left(\begin{array}{cc}\frac{\partial S\left(x\right)}{\partial \lambda } & \frac{\partial S\left(x\right)}{\partial \theta }\end{array}\right)$, ${\varDelta }_{h} = \left(\begin{array}{cc}\frac{\partial h\left(x\right)}{\partial \lambda } & \frac{\partial h\left(x\right)}{\partial \theta }\end{array}\right)$, with the following elements

According to the above conclusions, the $100(1-\alpha)\%$ ACIs of $S\left(t\right)$ and $h\left(t\right)$ for a given $t$ are

Since the asymptotic property of MPS estimator is similar to that of MLE estimator, this paper adopts the above method to obtain the interval estimation of MPS. For more details, see Ghosh and Jammalamadaka ^[24].

4.4.   Bayesian estimation

In this section, the Bayesian estimator is used to estimate the parameters of the NGL distribution, and the corresponding highest posterior density (HPD) intervals are considered. The idea of Bayesian estimation in this paper refers to several papers ^{[25,26,27,28]}. Assume that $\lambda$ and $\theta$ are independent and obey gamma distributions. The prior PDFs of $\lambda$ and $\theta$ are:

Based on these assumptions, the joint prior density function of $\lambda$ and $\theta$ can be obtained

According to the likelihood function with the prior knowledge and Bayes' theorem, the posterior density distribution of the unknown parameters $\lambda$ and $\theta$ can be obtained:

Using Eqs (33) and (47), the posterior density distribution can be expressed:

The loss function is a key factor to make decisions in Bayesian estimation. In this paper, the SE and GE loss function are used to measure overestimation and underestimation in the investigation. The SE and GE loss functions, as symmetric loss function and asymmetric loss function, respectively, have different measures of the importance of overestimation and underestimation. The Bayesian estimator represents the posterior mean in the case of the SE loss function, and overestimation and underestimation have equal weight. The SE loss is defined as ^[29]

and the corresponding Bayesian estimator is

The GE loss, which has a different tendency to weight overestimation and underestimation, is defined as ^[30]

here, $\gamma$ is the parameter of the degree of asymmetry. Under GE loss, the Bayesian estimator is

Obviously, the integral of Eqs (47) and (49) cannot be calculated directly. Therefore, the MCMC approach, which is a very popular method for estimating parameters, is employed to calculate the corresponding HPD intervals and the Bayesian estimates (BE) of $\lambda$ and $\theta$. The full conditional posterior distribution of $\lambda$ and $\theta$ as a key factor of the MCMC method can be derived by Eq (45)

and

Because of the nonlinearity of the full conditional posterior distribution of $\lambda$ and $\theta$, the Metropolis-Hastings (M-H) algorithm is applied to obtain the unknown parameters of Bayesian estimation. We assume the normal distribution as the proposed distribution to obtain the Bayesian estimation and HPD intervals of $\lambda$, $\theta$, $S\left(t\right)$ and $h\left(t\right)$. Follow the next steps to generate the MCMC sample:

(ⅰ) Set $k = 1$.

(ⅱ) Set the initial values of $(\lambda, \theta)$ to $({\lambda }^{\left(0\right)}, {\theta }^{\left(0\right)})$.

(ⅲ) Generate ${\lambda }^{*}$ and ${\theta }^{*}$ from $N({\lambda }^{(k-1)}, {\sigma }_{\lambda }^{2})$ and $N({\theta }^{(k-1)}, {\sigma }_{\theta }^{2})$, respectively. When ${\lambda }^{*}\le 0$ or ${\theta }^{*}\notin \left(\mathrm{0, 1}\right)$, repeat step (ⅲ), where ${\lambda }^{(k-1)}$ and ${\theta }^{(k-1)}$ represent previous state, ${\sigma }_{\lambda }^{2}$ and ${\sigma }_{\theta }^{2}$ represent the variance of the previous state.

(ⅳ) Definite acceptance probability $\omega ({\lambda }_{k-1}, {\lambda }^{*}) = \mathit{min}(1, \frac{{\pi }_{1}\left({\lambda }^{*}\right|{\theta }^{(k-1)}, x)}{{\pi }_{1}\left({\lambda }^{(k-1)}\right|{\theta }^{(k-1)}, x)})$ and $\omega ({\theta }_{k-1}, {\theta }^{*}) = \mathit{min}(1, \frac{{\pi }_{2}\left({\theta }^{*}\right|{\lambda }^{\left(k\right)}, x)}{{\pi }_{2}\left({\theta }^{(k-1)}\right|{\lambda }^{\left(k\right)}, x)})$.

(ⅴ) Generate ${u}_{\left(1\right)}$ and ${u}_{\left(2\right)}$ from the uniform distribution $U\left(\mathrm{0, 1}\right)$.

(ⅵ) If ${u}_{\left(1\right)}\le \omega ({\lambda }_{k-1}, {\lambda }^{*})$, then ${\lambda }^{\left(k\right)} = {\lambda }^{*}$, otherwise ${\lambda }^{\left(k\right)} = {\lambda }^{(k-1)}$.

(ⅶ) If ${u}_{\left(2\right)}\le \omega \left({\theta }_{k-1}, {\theta }^{*}\right)$, then ${\theta }^{\left(k\right)} = {\theta }^{*}$, otherwise ${\theta }^{\left(k\right)} = {\theta }^{\left(k-1\right)}$.

(ⅷ) Calculate ${S}^{\left(k\right)}\left(t\right)$ and ${h}^{\left(k\right)}\left(t\right)$ according to the following formulas: ${S}^{\left(k\right)}\left(t\right) = (1+{\lambda }^{\left(k\right)}{\theta }^{\left(k\right)}t){e}^{-{\lambda }^{\left(k\right)}t}$ and ${h}^{\left(k\right)}\left(t\right) = {\lambda }^{\left(k\right)}-\frac{{\lambda }^{\left(k\right)}{\theta }^{\left(k\right)}}{1+{\lambda }^{\left(k\right)}{\theta }^{\left(k\right)}t},$ where $t > 0$.

(ⅸ) Set $k = k+1$.

(ⅹ) Repeat (ⅲ)–(ⅸ) $L$ times to get $\left\{{\lambda }^{\left(k\right)}\right\}, \left\{{\theta }^{\left(k\right)}\right\}, \left\{{h}^{\left(k\right)}\right(t\left)\right\}$ and $\left\{{S}^{\left(k\right)}\right(t\left)\right\}(k = \mathrm{1, 2}, \dots, L)$, and discard the first ${L}^{*}$ samples of $\left\{{\lambda }^{\left(k\right)}\right\}, \left\{{\theta }^{\left(k\right)}\right\}, \left\{{h}^{\left(k\right)}\left(t\right)\right\}$ and $\left\{{S}^{\left(k\right)}\left(t\right)\right\}$ to eliminate the influence of initial value selection.

(xi) Based on SE and GE loss functions, calculate the Bayesian estimation and HPD intervals of $\lambda$, $\theta$, $S\left(t\right)$ and $h\left(t\right)$. Take $\lambda$ for example:

● Compute the Bayesian estimate of $\lambda$ : ${\widehat{\lambda }}_{S} = \frac{1}{l}{\sum }_{i = {L}^{*}+1}^{L}{\lambda }^{\left(i\right)}$ and ${\widehat{\lambda }}_{G} = \frac{1}{l}{\left({\sum }_{i = {L}^{*}+1}^{L}[{\lambda }^{\left(i\right)}{]}^{-\gamma }\right)}^{-1/\gamma },$ where $l = L-{L}^{*}$.

● Create the HPD interval of $\lambda$ ^[31]: Let ${\lambda }_{(L*+1)}, {\lambda }_{(L*+2)}, ..., {\lambda }_{\left(L\right)}$ be the ascending values of ${\lambda }^{({L}^{*}+1)}, {\lambda }^{({L}^{*}+2)}, ..., {\lambda }^{\left(L\right)}$, the $100(1-\alpha)\%$ HPD interval of $\lambda$ can be approximated to $\left({\lambda }_{\left({k}^{*}\right)}, {\lambda }_{({k}^{*}+[(1-\alpha)l\left]\right)}\right)$, where ${k}^{*}\in \{{L}^{*}+1, {L}^{*}+2, ..., L\}$ is selected according to the formula ${\lambda }_{({k}^{*}+[(1-\alpha)l\left]\right)}-{\lambda }_{\left({k}^{*}\right)} = \underset{{L}^{*}+1\le k\le {L}^{*}+\left[\alpha l\right]}{min}\left({\lambda }_{({k}^{*}+[(1-\alpha)l\left]\right)}-{\lambda }_{\left({k}^{*}\right)}\right)$, where $\left[x\right]$ is the downward integer of $x$, that is, the largest integer less than or equal to $x$.

5.   Monte Carlo simulation

In this section, a Monte Carlo simulation will be performed to demonstrate and compare the performance of the above estimators for the NGL distribution in parameter estimation.

5.1.   Simulation plan

We simulate 1000 progressively type-Ⅱ censored samples of the NGL(0.5, 0.75) based on the parameter selection of n(Total number of samples), m(Number of valid samples), and P(Censored schemes). Meanwhile, in order to reasonably evaluate the estimates of $S\left(t\right)$ and $h\left(t\right)$, first obtain their true values at $t = 0.5$, which are 0.9248 and 0.1842, respectively. In addition, $n( = \mathrm{50, 90})$ is determined and the valid sample proportion is used to determine that the value of $m$ meets $m/n( = 60\%, 80\%)$. Also, the three progressively censored schemes used are shown in Table 2.

Next, the specific steps to generate progressively type-Ⅱ censored samples are given ^[32]:

(ⅰ) Generate $m$ observations ${\varsigma }_{i}$ (for $i = \mathrm{1, 2}, ..., m$) that follow uniform distribution $U\left(\mathrm{0, 1}\right)$.

(ⅱ) Give the corresponding value according to the specific censored schemes, set ${\zeta }_{i} = {\varsigma }_{i}^{1/(i+{\sum }_{k = m-i+1}^{m}{P}_{k})}$ (for $i = \mathrm{1, 2}, ..., m$).

(ⅲ) Set ${\xi }_{i} = 1-{\zeta }_{m}{\zeta }_{m-1}\cdots {\zeta }_{m-i+1}$ (for $i = \mathrm{1, 2}, ..., m$).

(ⅳ) Generate progressively type-Ⅱ censored samples ${x}_{i}$ of the NGL (0.5, 0.75), set

In Table 2, P = (6, 0, 0, 0, 0, 1) stands for P = (6, 0*4, 1).

After obtaining progressively type-Ⅱ censored samples, MPS estimates and ML estimates and 95% ACIs of $\lambda$, $\theta$, $S\left(t\right)$ and $h\left(t\right)$ are calculated using Matlab R2016a, as well as Bayesian estimates based on SE and GE ($\gamma ( = -2)$) loss function and HPD interval. The large 12,000 M-H samples are generated by M-H sampler, and then the first 2000 samples are deleted as fluctuation samples. Then we complete the Bayesian estimator by setting up two prior sets called $Prior-a:({\sigma }_{1}, {\sigma }_{2}, {\omega }_{1}, {\omega }_{2}) = \left(\mathrm{8, 10, 10, 5}\right)$ and $Prior-b:({\sigma }_{1}, {\sigma }_{2}, {\omega }_{1}, {\omega }_{2}) = (\mathrm{4, 5}, \mathrm{5, 2.5})$. Repeat 1000 times to ensure the accuracy and stability of the estimation results. The evaluation criteria for verifying the reliability of the estimation method are shown in Table 3.

Where ${l}_{A}\left(\right)$ represents the indicator function, and ${R}_{\left(\right)}$ and ${L}_{\left(\right)}$ denote the upper and lower bounds, respectively, for each $(1-\alpha)\%$ ACI/HPD interval.

5.2.   Simulation result

RMSE, MRAB, ACL, and CP of $\lambda$, $\theta$, $S\left(t\right)$ and $h\left(t\right)$ estimates of various estimators are shown using heat maps in Figures 5–8. The corresponding specific numerical results are shown in the appendix [Tables 4–11]. Here, SE-Pa represents the estimate of $Prior-a$ by Bayesian estimator based on SE loss, GE-Pb represents the estimate of $Prior-b$ by Bayesian estimator based on GE loss (for $\gamma ( = -2)$), and ACL-MPS corresponds to the interval estimate of the MPS.

From the heat maps in Figures 5–8 and Tables 4–11, the following conclusions can be drawn:

(1) All estimates of $\lambda$, $\theta$, $S\left(t\right)$ and $h\left(t\right)$ are good estimators because of low RMSE, MRAB, and ACL values and high CP values. Through changing the color of the heat-maps to from down to up, we can find that, in most cases, with the increase of n and m, the estimation performance of all obtained estimators will improve, corresponding to lower RMSE, MRAB, ACL, and CP values.

(2) Among all estimates, Bayesian estimates based on SE loss and GE loss are more accurate than MPS and ML estimates because of lower RMSE, MRAB, and ACL values and higher CP values.

(3) Different prior parameters will affect the effectiveness of Bayesian estimation. In Bayesian estimation, the Bayesian estimator based on SE loss and the Bayesian estimator based on GE loss have better estimation performance under the gamma prior function with $Prior-a$ as parameter than $Prior-b$ because the variance of $Prior-a$ is smaller.

(4) Different censored schemes may affect the estimated effectiveness to some extent. In most cases, estimates based on scheme-1 work better. In the face of progressively type-Ⅱ censored samples of the NGL distribution, the Bayesian estimation under the gamma prior function with $Prior-a$ can be used to estimate the corresponding unknown parameters and related functions.

(5) The Bayesian estimation under the gamma prior function with $Prior-a$ is overestimated while with $Prior-b$ is underestimated.

(6) The point (or interval) estimation of the MPS estimates and ML estimates of $S\left(t\right)$ is significantly weaker than the result obtained by Bayes.

(7) In summary, with a more complete amount of data from progressively type-Ⅱ censored, Bayesian estimation via M-H algorithm can obtain better estimates when estimating unknown parameters of the NGL distribution.

6.   Practical application

This section demonstrates the flexibility of the proposed distribution and the usefulness of the various estimation methods through a practical application. The dataset, which was previously used by Bekker et al. ^[33] and later applied by Habib et al. ^[34], is the annual survival time of 46 patients who received chemotherapy and radiation therapy. The specific data is shown in Table 12.

In order to test whether the NGL distribution is suitable for the dataset, ML estimation is first used to estimate $\lambda$ and $\theta$, and the corresponding values with standard errors (St) are obtained as 0.90788 (0.33141) and 0.20264 (0.41243). Then the Kolmogorov-Smirnov (KS) value and P-value can be obtained as 0.10158 (0.69157). Compared with the truncated Nadarajah-Haghighi Raykeigh distribution, which presents a KS value with a P-value of 0.1080 (0.6307) for this dataset ^[34], it can be considered that the NGL distribution has a better fitting effect and is suitable for this dataset. As illustrated in Figure 9, the NGL distribution to the fitting effect of the dataset is demonstrated. This includes the fitted CDF, the probability-probability (PP), the scaled total time on test (TTT) transform ^[35], and the counter of the log-likelihood function. Figure 9 indicates that the NGL distribution is very close to the real data distribution and the existence and uniqueness of the obtained MLE estimates $\widehat{\lambda }$ and $\widehat{\theta }$.

Let m = 20, and three types of progressively type-Ⅱ censored samples are obtained from this dataset, as shown in Table 13. The point estimates of $\lambda$ and $\theta$ are obtained by the above estimation method (MPS, ML, Bayesian estimation) and are shown in Table 14. At the same time, standard errors (St.es) are used to judge the accuracy of the estimated results. Here, the SE.es of $\lambda$ and $\theta$ obtained by MPS and ML estimation are the square root of the diagonal elements of ${I}^{-1}\left(\Theta \right)$, respectively.

Through the M-H algorithm, 5000 MCMC samples are generated, and the first 1000 samples are omitted. The sample sequence $\left\{{\lambda }_{\left(i\right)}\right\}, \left\{{\theta }_{\left(i\right)}\right\}$ formed can approximately obey the posterior distribution, and the variance of the posterior distribution can be estimated by the sample variance, that is:

It can be seen from Table 14 that under three different censored schemes:

● All estimation methods have a better fitting effect when fitting the dataset.

● The results estimated by ML and MPS are similar, and the results estimated by Bayes under two different losses are similar.

● The MPS and ML estimation fit the scale parameter $\lambda$ preferentially when fitting, which makes the scale parameter more accurate but tends to lead to overflow of the shape parameter $\alpha$. Bayesian estimation weighs the two parameters more than the other two estimation methods and is less prone to overflow of the range of estimates.

● In general, the results of Bayesian estimates based on GE loss functions are more consistent with the dataset than those of other methods.

7.   Conclusions

A Lindley distribution with multiple models can be applied to various fields, such as product life research and species presence distribution. This paper focused on the examination of the NGL distribution's inherent properties and significant statistical characteristics. It was determined that the proposed distribution possesses favorable properties and that such properties exhibit relatively straightforward numerical expressions. This finding is instrumental in facilitating subsequent research endeavors. Subsequently, this paper actively explored the application of maximum product spacing estimation, maximum likelihood estimation, and Bayesian estimation to the NGL distribution. The development of a comprehensive simulation plan, incorporating two distinct failure rates, three disparate censored schemes, and four evaluation criteria, was undertaken to investigate the impact of various estimation methods on the proposed distribution point estimation and interval estimation. Finally, it was concluded that the Bayesian estimator under SE loss and GE loss has relatively good estimation effect, and this estimation method can be integrated into the NGL distribution for further exploration in future research. In addition, in order to ensure that the NGL distribution has practical application significance, a real dataset, containing annual survival rates, was used to carry out research. The analysis showed that the proposed distribution has a good fitting effect on the original dataset, and we found that the proposed distribution has a good fit to the censored data of the dataset. The above estimation methods were used to finally conclude that the distribution also has a good fitting effect on the censored data of the real dataset in the Bayesian estimation. In future studies, we will continue to increase the integration of the proposed distribution with real industry scenarios.

Author contributions

Jiajie Shi: Writing-original draft, Method, Software, Writing-review & editing; Haiping Ren: Methods, Writing-original draft, Writing-review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was funded by National Natural Science Foundation of China, grant number 71661012.

Conflict of interest

The authors declare no conflict of interest.

Appendix

A.   The parameters, the reliability function, and hazard function of the NGL distribution are estimated using MPS estimation

clear

tic

a0=0.5; b0=0.75; m=30;n=50;num=1000;t=0.5;S0=0.9248;h0=0.1842;ca=0;cb=0;cS=0;ch=0;

Ra=[zeros(1,m-1),n-m];

Rb=[floor((n-m)/2),zeros(1,m-2),(n-m-floor((n-m)/2))];

Rc=[n-m,zeros(1,m-1)];

a = zeros(num,1);

b = zeros(num,1);

ACI1 = zeros(num, 2);

ACI2 = zeros(num, 2);

ACI3 = zeros(num, 2);

ACI4 = zeros(num, 2);

for i=1:num

G=rand(1,m);

for j=1:m

H(j)=G(j).^(1/(j+sum(Ra(m-j+1:m))));

end

for j=1:m

Z(j)=1-prod(H(m-j+1:m));

end

x=(-b0.*lambertw(-1,exp(-1./b0).*(Z-1)./b0)-1)./(a0.*b0);

[a(i,1),b(i,1),ACI1(i,:),ACI2(i,:),S(i,1),h(i,1),ACI3(i,:),ACI4(i,:)] =MPS(t,Ra,n,x,a0,b0);

if a0>=ACI1(i,1)&&a0<=ACI1(i,2)

ca=ca+1;

end

if b0>=ACI2(i,1)&&b0<=ACI2(i,2)

cb=cb+1;

end

if S0>=ACI3(i,1)&&S0<=ACI3(i,2)

cS=cS+1;

end

if h0>=ACI4(i,1)&&h0<=ACI4(i,2)

ch=ch+1;

end

end

aMLi=mean(a)

mse_a=sqrt(mean((a-a0).^2))

MRE_a=mean(abs((a-a0)./a0))

ACL_a=mean(ACI1(:,2)-ACI1(:,1))

ca=ca/num

bMLi=mean(b)

mse_b=sqrt(mean((b-b0).^2))

MRE_b=mean(abs((b-b0)./b0))

ACL_b=mean(ACI2(:,2)-ACI2(:,1))

cb=cb/num

SMLi=mean(S)

mse_S=sqrt(mean((S-S0).^2))

MRE_S=mean(abs((S-S0)./S0))

ACL_S=mean(ACI3(:,2)-ACI3(:,1))

cS=cS/num

hMLi=mean(h)

mse_h=sqrt(mean((h-h0).^2))

MRE_h=mean(abs((h-h0)./h0))

ACL_h=mean(ACI4(:,2)-ACI4(:,1))

ch=ch/num

toc

B.   The parameters, the reliability function, and hazard function of the NGL distribution are estimated using ML estimation

clear

tic

a0=0.5; b0=0.75; m=30;n=50;num=1000;t=0.5;S0=0.9248;h0=0.1842;ca=0;cb=0;cS=0;ch=0;

Ra=[zeros(1,m-1),n-m];

Rb=[floor((n-m)/2),zeros(1,m-2),(n-m-floor((n-m)/2))];

Rc=[n-m,zeros(1,m-1)];

a = zeros(num,1);

b = zeros(num,1);

S=zeros(num,1);

h=zeros(num,1);

ACI1 = zeros(num, 2);

ACI2 = zeros(num, 2);

ACI3 = zeros(num, 2);

ACI4 = zeros(num, 2);

for i=1:num

G=rand(1,m);

for j=1:m

H(j)=G(j).^(1/(j+sum(Ra(m-j+1:m))));

end

for j=1:m

Z(j)=1-prod(H(m-j+1:m));

end

x=(-b0.*lambertw(-1,exp(-1./b0).*(Z-1)./b0)-1)./(a0.*b0);

[a(i,1),b(i,1),ACI1(i,:),ACI2(i,:),S(i,1),h(i,1),ACI3(i,:),ACI4(i,:)] =MLi(t,Ra,n,x,a0,b0);

if a0>=ACI1(i,1)&&a0<=ACI1(i,2)

ca=ca+1;

end

if b0>=ACI2(i,1)&&b0<=ACI2(i,2)

cb=cb+1;

end

if S0>=ACI3(i,1)&&S0<=ACI3(i,2)

cS=cS+1;

end

if h0>=ACI4(i,1)&&h0<=ACI4(i,2)

ch=ch+1;

end

end

aMLi=mean(a)

mse_a=sqrt(mean((a-a0).^2))

MRE_a=mean(abs((a-a0)./a0))

ACL_a=mean(ACI1(:,2)-ACI1(:,1))

ca=ca/num

bMLi=mean(b)

mse_b=sqrt(mean((b-b0).^2))

MRE_b=mean(abs((b-b0)./b0))

ACL_b=mean(ACI2(:,2)-ACI2(:,1))

cb=cb/num

SMLi=mean(S)

mse_S=sqrt(mean((S-S0).^2))

MRE_S=mean(abs((S-S0)./S0))

ACL_S=mean(ACI3(:,2)-ACI3(:,1))

cS=cS/num

hMLi=mean(h)

mse_h=sqrt(mean((h-h0).^2))

MRE_h=mean(abs((h-h0)./h0))

ACL_h=mean(ACI4(:,2)-ACI4(:,1))

ch=ch/num

toc

C.   The parameters, the reliability function, and hazard function of the NGL distribution are estimated using Bayesian estimation

clear

a0=0.5; b0=0.75; m=30;n=50;num=1000;t=0.5;S0=0.9248;h0=0.1842;c1=0;c2=0;c3=0;c4=0;

Ra=[zeros(1,m-1),n-m];

Rb=[floor((n-m)/2),zeros(1,m-2),(n-m-floor((n-m)/2))];

Rc=[n-m,zeros(1,m-1)];

a_BMH=zeros(num,1);

b_BMH=zeros(num,1);

S_BMH=zeros(num,1);

h_BMH=zeros(num,1);

aHPD=zeros(num,1);

bHPD=zeros(num,1);

SHPD=zeros(num,1);

hHPD=zeros(num,1);

for i=1:num

G=rand(1,m);

for j=1:m

H(j)=G(j).^(1/(j+sum(Rc(m-j+1:m))));

end

for j=1:m

Z(j)=1-prod(H(m-j+1:m));

end

x=(-b0.*lambertw(-1,exp(-1./b0).*(Z-1)./b0)-1)./(a0.*b0);

[a_BMH(i,1),b_BMH(i,1),S_BMH(i,1),h_BMH(i,1),aHPD(i,1),bHPD(i,1),SHPD(i,1),hHPD(i,1),cp1,cp2,cp3,cp4]=BMH(t,Rc,n,x,0.4,0.7);

c1=c1+cp1;c2=c2+cp2;c3=c3+cp3;c4=c4+cp4;

end

aBMH=mean(a_BMH)

mse_a=sqrt(mean((a_BMH-a0).^2))

MRE_a=mean(abs((a_BMH-a0)./a0))

ACL_a=mean(aHPD)

ca=c1/num

bBMH=mean(b_BMH)

mse_b=sqrt(mean((b_BMH-b0).^2))

MRE_b=mean(abs((b_BMH-b0)./b0))

ACL_b=mean(bHPD)

cb=c2/num

SBMH=mean(S_BMH)

mse_S=sqrt(mean((S_BMH-S0).^2))

MRE_S=mean(abs((S_BMH-S0)./S0))

ACL_S=mean(SHPD)

cS=c3/num

hBMH=mean(h_BMH)

mse_h=sqrt(mean((h_BMH-h0).^2))

MRE_h=mean(abs((h_BMH-h0)./h0))

ACL_h=mean(hHPD)

ch=c4/num

Counter diagram of the parameters of the NGL distribution

clear

syms a b;

a0=0.5; b0=0.75; m=50;n=100;

Ra=[zeros(1,m-1),n-m];

Rb=[floor((n-m)/2),zeros(1,m-2),(n-m-floor((n-m)/2))];

Rc=[n-m,zeros(1,m-1)];

G=rand(1,m);

for j=1:m

H(j)=G(j).^(1/(j+sum(Ra(m-j+1:m))));

end

for j=1:m

Z(j)=1-prod(H(m-j+1:m));

end

P=Ra;

x=(-b0.*lambertw(-1,exp(-1./b0).*(Z-1)./b0)-1)./(a0.*b0);

lb=m.*log(a)-a.*sum(x.*(1+P))+sum(log(1-b+a.*b.*x))+sum(P.*log(1+a.*b.*x));

lc= n.*log(a)-a.*sum(x)+sum(log(1-b+a.*b.*x));

laa=matlabFunction(lb, 'Vars', {a, b});

a1 = 0.2:0.01:1;

b1 = 0.36:0.01:0.75;

[a1,b1] = meshgrid(a1,b1);

la=real(laa(b1,a1));

[M,c]=contour(a1,b1,la,'showText','on');

xlabel('θ');

ylabel('λ');

c.LineWidth=3;

colorbar.