Survival and Reliability Analysis with an Epsilon-Positive Family of Distributions with Applications

We introduce a new class of distributions called the epsilon–positive family, which can be viewed as generalization of the distributions with positive support. The construction of the epsilon– positive family is motivated by the ideas behind the generation of skew distributions using symmetric kernels. This new class of distributions has as special cases the exponential, Weibull, log–normal, log–logistic and gamma distributions, and it provides an alternative for analyzing reliability and survival data. An interesting feature of the epsilon–positive family is that it can viewed as a finite scale mixture of positive distributions, facilitating the derivation and implementation of EM–type algorithms to obtain maximum likelihood estimates (MLE) with (un)censored data. We illustrate the flexibility of this family to analyze censored and uncensored data using two real examples. One of them was previously discussed in the literature; the second one consists of a new application to model recidivism data of a group of inmates released from the Chilean prisons during 2007. The results show that this new family of distributions has a better performance fitting the data than some common alternatives such as the exponential distribution.


Introduction
The statistical analysis of reliability and survival data is an important topic in several areas, including medicine, epidemiology, biology, economics, engineering, and environmental sciences, to name a few. When using a parametric approach, one of the first steps for modeling the data is to choose a suitable distribution that can capture relevant features of the observations of interest. In this context, the gamma and Weibull distributions have become popular choices due to their flexibility that allows for a non-constant hazard rate function and to model skewed data. Although several alternatives were considered to accommodate different cases, researchers have continued to develop extensions and modifications of the standard distributions to increase the flexibility of the models see [1][2][3] for a few examples.
In this paper, we consider a generalization of the distributions with positive support and propose a new family of distributions, called the epsilon-positive family, whose construction is motivated by the ideas behind the generation of skew distributions using symmetric kernels. Specifically, we build upon the ideas from [4], where the authors start with a symmetric around zero distribution f , and define a family of distributions indexed by a parameter γ > 0 as the set of densities of the form h(x; γ) = 2 where 1 A denotes the indicator function of the set A. Some extensions of this family include epsilon-skew-normal family introduced by [5] and the epsilon-skew-symmetric family introduced by [6], both discussed in some details in [7].
Here, starting with a probability density function g with positive support, we obtain a general class that extends the family of distributions with positive support, and that contains the Weibull, gamma and exponential distributions as special cases, depending on the choice of g. Furthermore, we discuss a stochastic representation and how to obtain maximum likelihood estimators for the members of this family. We also derive the corresponding survival and hazard functions and note that one interesting feature of this new class is that the hazard function is not necessarily constant.
The rest of the paper is organized as follows: in Section 2 we define the epsilonpositive family and obtain the hazard and survival functions, mean residual life and stressstrength parameters for this family. In addition, we discuss maximum likelihood estimation and how to obtain such estimates using an EM-type algorithm for the general case. In Section 3, we focus on one specific member of the family introduced in Section 2, namely epsilon-exponential distribution, and discuss its applicability in the analysis of survival data. In Section 4 we discuss two real data examples and we finish with a brief discussion in Section 5. We include Appendixes A and B with some of the technical details.

The Epsilon-Positive Family
Let g(·) = g Y (·; Ψ) be a probability density function (pdf) with positive support and parameters Ψ ∈ p . Then, for 0 < ε < 1, the corresponding epsilon-positive (EP) family of distributions is defined as If a random variable X has the density given in (2), we say that X has an epsilonpositive distribution and write X ∼ EP(Ψ, ε).
Observe that as ε ↓ 0, f X (x; Ψ, ε) → g(x)1 {x>0} and therefore the distribution g Y (·; Ψ) can be seen as a particular member of the family.
To draw observations from an epsilon-positive distribution, we first notice that for From this stochastic representation, it follows that we can generate EP random variables according to the Algorithm 1: Algorithm 1 Algorithm to generate observations from an epsilon-positive distribution.

Reliability Properties
From the definition, we obtain that the survival function S X (x; Ψ, ε) = P(X > x) for this family is given by where S Y (·) is the survival function associated with the density g Y (·; Ψ). Similarly, the hazard function λ X (x; Ψ, ) = f X (x; Ψ, ε)/S X (x; Ψ, ε) is given by where , and R(·) = S Y (·)/g Y (·) is the Mills ratio. Table 1 shows some examples of the densities that can be extended using the definition of the epsilon-positive family, with the corresponding densities, survival and hazard functions. Figures 1 and 2 show the pdf, survival and hazard functions of the epsilonexponential, epsilon-Weibull, epsilon-log-logistic and epsilon-gamma distributions. We can see that in the case of the epsilon-Weibull, epsilon-log-logistic and epsilon-gamma distributions a bimodal shape is obtained when the value of the parameter ε is 0.9.
Examples of the probability density f (x), survival S(x) and hazard λ(x) functions of epsilon-log-logistic distribution, ELL(σ, ε), and epsilon-gamma distribution, EG(α, σ, ε). Please note that the log-logistic and gamma distributions correspond to the case ε = 0. Table 1. Hazard rate, λ(·), Survival, S(·), and density, f (·), functions of some probability models that can be generalized using the definition in (2) In the table

Mean Residual Life
The mean residual life or life expectancy is an important characteristic of the model. It gives the expected additional lifetime given that a component has survived until time t. For a non-negative continuous random variable X ∼ EP(Ψ, ε) the mean residual life (mlr) function is defined as where t > 0. The above conditional expectation is given by Calculation of the integral in (8) is done in the same way as the calculation of the mean. Thus, Making the changes of variables z = x 1+ε and u = x 1−ε we have (7) can be written as

Stress-Strength Parameter
An important concept in reliability theory is the stress-strength parameter. Let X 1 denote the strength of a system or component with a stress X 2 . Then, the stress-strength parameter is defined as R = P(X 2 < X 1 ), which can be viewed as a measure of the system performance. In the next theorem, we look at this quantity when X 1 and X 2 are independent random variables with epsilon-positive distributions. Theorem 1. Suppose X 1 and X 2 are random variables independently distributed as X 1 ∼ EP(Ψ 1 , ε 1 ) and X 2 ∼ EP(Ψ 2 , ε 2 ), the reliability of the system with stress variable (X 2 ) and strength variable (X 1 ) is given by where a = 1+ε 1 Proof of Theorem 1. Making the changes of variables z = x 1 1+ε and u = x 1 1−ε we have Observe that the same concept can be used to make comparisons between two systems. For example, if X 1 and X 2 denote instead the lifetimes of systems S 1 and S 2 respectively, then, a probability P(X 1 < X 2 ) > 0.5 would indicate that the system S 2 is better than the system S 1 in a stochastic sense.

Maximum Likelihood Estimation
LetX n = (X 1 , . . . , X n ) be a random sample from an EP(Ψ, ε) distribution. Then, the maximum likelihood estimator (MLE) of θ = (Ψ, ε) is given by is the log-likelihood. Although the MLE for the EP family is conceptually straightforward, typically closed form solutions are not available and the MLE need to be obtained numerically. One possibility is the Newton-Raphson algorithm, with iteration equation where θ (k) be the current estimate of θ, u(θ) denote the vector of first derivatives of (θ;X n ), and H(θ). A disadvantage of this approach is that it requires redthe calculation of the second derivatives of the likelihood function and repeated inversion of potentially large matrices, which can be computationally intensive. Instead, we can consider an expectationmaximization (EM) approach see [8] as a general iterative method for data sets with missing (or incomplete) data.
The mixture representation proposed in (4) is particularly useful in order to use an EMtype algorithm to estimate the model parameters, since it provides a hierarchical scheme for the EP family. Next, we show how to implement maximum likelihood estimation using an EM-type algorithm for the EP family.
The E and M steps are alternated repeatedly until a convergence criteria is satisfied. For the variance estimation of the MLEs we consider the bootstrapping method suggested in [9].

The Epsilon-Exponential Distribution
If we take g Y (y; Ψ = σ) = 1 σ e −y/σ 1 {y>0} , the pdf of an exponential distribution, the expression in (2) becomes where σ > 0 and 0 < ε < 1. We say that a random variable X has epsilon-exponential (EE) distribution with scale parameter σ and shape parameter ε if its density has the form in (13), and we write X ∼ EE(σ, ε).
Recall that the rth moment, r = 1, 2, . . . , of Y ∼ Exp(σ) is E(Y r ) = r!σ r . From (3), when X ∼ EE(σ, ε), we obtain that the mean, variance, skewness (CS) and kurtosis (CK) coefficients are, respectively, Please note that for any value of ε ∈ (0, 1), CS > 0 and CK > 0. It can be seen that 2 < CS < 2.3 and 9 < CK < 11.023. Figure 3 depicts the behavior of the skewness (CS) and kurtosis (CK) coefficients as a function of ε. In the figures we observe that the maximum skewness is attained ar ε = 0.4, while the maximum kurtosis coefficient is obtained when ε = 0.37. Recall that the survival function of the exponential distribution is of the form S Y (y) = e −y/σ 1 {y>0} and the mean residual life is of the form mrl Y (t) = σ. Then, it follows from (5), (6) and (9) that if X ∼ EE(σ, ε), then the survival and hazard functions and mean residual life are (respectively) Interestingly, in contrast to the exponential distribution, it can be shown that the hazard function λ(x; σ, ε) of the EE distribution is not constant, but decreasing in x. This feature can be easily observed in Figure 1 (top panel), where we show the pdf, survival and hazard functions of the EE distribution for different values of the parameter ε when σ = 1. Please note that λ(x; σ, ε) −→ λ Y (x; σ) = 1/σ as ε −→ 0. Additionally, mrl(t) −→ mrl Y (t) = σ as ε −→ 0. Suppose X 1 and X 2 are random variables independently distributed as X 1 ∼ EE(σ 1 , ε 1 ) and X 2 ∼ EE(σ 2 , ε 2 ), using Theorem 1 the reliability of the system with stress variable (X 2 ) and strength variable (X 1 ) is given by Please note that when (ε 1 , ε 2 ) −→ (0, 0), R −→ σ 1 σ 1 +σ 2 which corresponds to the stress-strength reliability model of the exponential distributions with parameter σ 1 for X 1 (strength), and with parameter σ 2 for X 2 (stress), respectively.

Numerical Experiments
To illustrate the properties of the estimators we performed a small simulation study considering 5000 simulated datasets for different pair of values of σ and ε using Algorithm 1. The goal of the study is to observe the behavior of the MLEs for the model parameters using our proposed EM algorithm considering different sample sizes. Table 2 summarizes the simulation results. In the table, the columns labeled as "estimate" show the average of the estimators obtained in the simulations, and the columns labeled "SD" show the sample standard deviation of the corresponding estimators. To obtain the standard errors we used the bootstrap method with B = 150 samples. We observe that the estimates are quite stable and fairly accurate, reporting (on average) numbers close to the true values of the parameters in all cases. Please note that as expected, the precision of the estimates improves as the sample size increase.

Survival and Reliability Analysis
Let T ≥ 0 represent the survival time until the occurrence of a "death" event. In this context, suppose we have n subjects with lifetimes determined by a survival function S(t), and that the ith subject is observed for a time t i . If the individual dies at time t i , its contribution to the likelihood function is given by is the corresponding hazard function. On the other hand, if the ith individual is still alive at time t i and therefore, under non-informative censoring, all we can say is that their lifetime exceeds t i . It follows that the contribution of a censored observation to the likelihood is simply given by L i = S(t i ). Notice that in either case we evaluate the survival function at time t i , because in both cases the ith subject was alive until (at least) time t i . A death will multiply this contribution by the hazard λ(t i ), but a censored observation will not.
We can combine these contributions into a single expression using a death indicator d i , taking the value one if individual i died and the value zero otherwise. The resulting likelihood function L is of the form In the next section we will assume that the random variable T follows an epsilonpositive family, and show how estimate the model parameters using the EM algorithm.

Estimation Using the EM Algorithm
Let T 1 , . . . , T n denote the survival times, T i ∼ EP(Ψ, ε). Using the notation introduced in the previous section, the observed data are a collection of pairsX n = { (T 1 , d 1 ), . . . , (T n , d n )}, where the d i , i = 1, . . . , n, are the censoring indicators.

EM Algorithm for the Epsilon-Exponential Distribution
Suppose that the survival times T i ∼ EE(σ, ε). Then the EM algorithm takes the form: • E-step: For i = 1, . . . , n, compute • M-step: Given ε (s) and σ (s) , compute

Real Data Examples
In this section, we use two examples to illustrate the proposed distributions using (un)censored data sets.

Example 1: Maintenance Data
First, we consider a real data set originally analyzed by [10]. The complete data set correspond active repair times (in hours) for an airborne communication transceiver, and can be found in Table 3. Using the EM algorithm described in Section 4.1.1 we fit an epsilon-exponential (EE) distribution to the active repair times. We obtain that the maximized log-likelihood value −103.806. Alternatively, we also fit an exponential (Exp), exponentiated-exponential (EExp) and Weibull (Wei) distribution, obtaining the maximized log-likelihood values of −105.006, −104.983 and −104.470, respectively. For model comparison, we use the Akaike information criterion (AIC) introduced in [11], given by AIC = −2l + 2k wherê l is the maximized log-likelihood and k is the number of parameters of the distribution under consideration.
The best model is deemed to be the one with the smallest AIC. For the data set in the example, we obtain AIC EE = 211.611, AIC Exp = 212.012, AIC EExp = 213.966, and AIC Wei = 212.939. It follows that in terms of the AIC criteria, the epsilon-exponential shows the best performance when fitting these data. Figure 4 shows the fit of the different models used in the example. Figure 5 displays the three estimated survival functions for this data set.

Example 2: Recidivism Data
For the second example we use real data obtained from the official records of Gendarmerie of Chile on all inmates released from the Chilean prisons during 2007 after serving a sentence of imprisonment by robbery.
The data set contains records of 9477 inmates and the follow-up period from release until 30 April 2012. In this study, recidivism is defined to occur when a released prisoner goes back to prison for the original or any other offense.
Overall, 52.2% of the inmates in the cohort were convicted of one or more offenses and returned to prison within 64 months of release. About 11.8% of the cohort returned to prison within three months of release, and 30% returned within a year of release. Table 4 shows the observed proportion of the cohort returning to prison within 1, 3, 6, 12, 18, 24, 36, 48 and 64 months of release. We observe that the cumulative proportion of recidivism grew quickly over the first 12 months after release, increasing by more than 7% every 3 months. After that, the percent increases were smaller and over longer periods. To analyze the time to recidivism, we determined the number of days between an inmate's release and his return to prison. Because some inmates did not reoffend, we have censored data and we used the EM algorithm described in Section 4.1.1 to fit an epsilon-exponential distribution to the time to recidivism. The maximized log-likelihood value for an assumed epsilon-exponential distribution is easily calculated to be −42,067.84. In comparison, we also fit an exponential distribution yielding a maximized log-likelihood value of −42,632.81.
Looking at the AIC values, we obtain AIC EE = 84,139.68 and AIC Exp = 85,267.62 for the epsilon-exponential and exponential model respectively, and therefore we conclude, the epsilon-exponential is a better model for these data, based on this criteria.
Finally, we also analyzed the survival time using the Kaplan-Meier estimator. Figure 5 displays the three estimated survival functions for this data set. We observe a close agreement between the Kaplan-Meier survival curve and the epsilon-exponential distribution.

Discussion
We introduced a new class of distributions with positive support called epsilonpositive which are generated from any distributions with positive support. This new class of distributions includes the exponential, Weibull, log-normal, etc. ones as special cases. We discussed a stochastic representation for this family, as well as parameter estimation, using the maximum likelihood approach via the Newton-Raphson. In addition, we showed that the elements of this new family can be expressed as a finite scale mixture of positive distributions, which facilitates the implementation of EM-type of algorithms.
We then focused on particular member of this family, called epsilon-exponential distribution, and discuss its applicability in the analysis of survival and reliability data. In this context, we considered the censored data case, and we show how we can use this new family to analyze this type of data sets. For the new class of distributions and, in particular, for the epsilon-exponential distribution we estimate the model parameters using the EM algorithm. An interesting feature of the hazard function of the epsilon-exponential distribution is that is not constant at difference of the exponential one. This feature increases the flexibility of the models allowing its use in a broader range of scenarios. This greater flexibility is corroborated in the two examples considered in this paper where the AIC criteria shows a better performance our proposed epsilon-exponential model when compared to commonly used alternatives such as the exponential one. The results suggests that the epsilon-exponential distribution should be considered to be a legitimate alternative for the analysis of survival and reliability data in both the censored and uncensored case.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.