APP下载

Robust and Efficient Reliability Estimation for Exponential Distribution

2021-12-15MuhammadAslamMohdSafariNurulkamalMasseranandMuhammadHilmiAbdulMajid

Computers Materials&Continua 2021年11期

Muhammad Aslam Mohd Safari,Nurulkamal Masseran and Muhammad Hilmi Abdul Majid

1Department of Mathematics,Faculty of Science,Universiti Putra Malaysia,UPM Serdang,43400,Selangor,Malaysia

2Department of Mathematical Sciences,Faculty of Science and Technology,Universiti Kebangsaan Malaysia,UKM Bangi,43600,Selangor,Malaysia

Abstract:In modeling reliability data,the exponential distribution is commonly used due to its simplicity.For estimating the parameter of the exponential distribution,classical estimators including maximum likelihood estimator represent the most commonly used method and are well known to be efficient.However,the maximum likelihood estimator is highly sensitive in the presence of contamination or outliers.In this study,a robust and efficient estimator of the exponential distribution parameter was proposed based on the probability integral transform statistic.To examine the robustness of this new estimator,asymptotic variance,breakdown point,and gross error sensitivity were derived.This new estimator offers reasonable protection against outliers besides being simple to compute.Furthermore,a simulation study was conducted to compare the performance of this new estimator with the maximum likelihood estimator,weighted likelihood estimator,and M-scale estimator in the presence of outliers.Finally,a statistical analysis of three reliability data sets was conducted to demonstrate the performance of the proposed estimator.

Keywords:Exponential distribution;M-estimator;probability integral transform statistic;robust estimation;reliability

1 Introduction

Exponential distribution is the most widely used parametric distribution for modeling reliability and failure time data due to its mathematical simplicity and ability to create a realistic failure time model [1-5].The advantages of applying a parametric model like exponential distribution in modeling failure time data are the following:this distribution can be concisely described with only one parameter instead of having to report an entire curve apart from providing smooth estimates of failure time distributions [1].Besides the reliability analysis,the exponential distribution is also applied for modeling income distribution [6,7].In the previous studies,research to propose a new model related to the family of exponential distributions has attracted considerable interest from many researchers [8-20].The main objective of introducing an extension and modification to the exponential distribution is to offer more flexible distribution structures for fitting data.

AssumeXas a random variable that follows the exponential distribution.Thus,the respective probability density function (PDF),cumulative distribution function (CDF),and survival function of the exponential distribution can be defined as follows:

and

whereγis the rate or inverse scale parameter.A high value ofγindicates high risk and short survival while a low value ofγindicates low risk and long survival [3].The parameterγalso measures the heaviness of exponential tail whereby the tail becomes heavier asγdecreases.In the analysis of reliability data,reliability properties are often defined using the mean time to failure(MTTF),survival (or reliability) function and failure rate function.Based on the exponential distribution,the survival function is given in Eq.(3),whereas the MTTF and failure rate function are expressed as the following:

and

The reliability data can be sometimes contaminated with outliers.Outliers are observations much deviated from the bulk of the data [21]and are also referred to as abnormalities,discordants,deviants,or anomalies in statistics literature [22].Outliers may arise due to errors or simply by natural deviations in a data set.In statistical modeling using parametric distribution,the classical estimator known as maximum likelihood estimator (MLE) is the most commonly used method for estimating the parameters of any parametric model.In fact,for any typical parametric distribution,MLE is well known to be efficient especially for a large sample size.However,in the presence of data contamination in which the outliers are present in the data set,the MLE is not robust and severely biased [23].Consequently,this condition may affect the application of exponential model in reliability analysis or any assessment.Therefore,when the outliers are present,a robust estimator should be utilized for estimating the parameter of any parametric distribution as an alternative to the MLE.

Several robust methods have been proposed in the literature for estimating the parameter of exponential distribution.Thall [24]proposed a robust estimator for the parameter of exponential distribution called Huber-sense robust M-estimator.However,this estimator is only suitable for the case when the proportion of outliers is small.Trimmed-mean and Winsorized-mean type estimators have also been proposed and studied by several researchers [25-27].The advantage of these types of estimator is that they are simple to compute.Further,Gather et al.[28]proposed another simple robust estimator called standardized median estimator (SME).However,the SME has low asymptotic relative efficiency (ARE) of only 48%.Gather et al.[28]also compared the performance of SME with two other estimators called RCS-estimator and RCQ-estimator in which the respective AREs for both estimators were 55% and 74%.Meanwhile,Ahmed et al.[29]introduced another robust estimator called the weighted likelihood estimator (WLE).When estimating the parameter of exponential distribution using WLE,observations with small likelihood are assigned with zero weights.This procedure causes the outliers to be eliminated from the data set and then,the ordinary MLE is used to estimate the parameter of interest based on the remaining observations.In recent years,Shahriari et al.[23]proposed a robust M-scale estimator for the parameter of exponential distribution in which this estimator has a maximum ARE of 71%.This estimator uses families of bisquare function as the weight for estimating the parameter of interest.

This study aims to develop a new robust estimator for the parameter of exponential distribution based on the probability integral transform statistic,which offers reasonable protection against outliers.Generally,probability integral transform statistic is used to transform the random variable of any continuous distribution to a standard uniform distribution [30-32].It should be noted that the new robust estimator proposed in this study is a class of M-estimator.This estimator has a simple form and easy to compute.In this study,the robustness of this new estimator was examined based on the measures of asymptotic variance,breakdown point,and gross error sensitivity.

The rest of this paper is structured as follows.Section 2 provides a brief explanation of M-estimators and then presents the new robust estimation method for the rate parameter of exponential distribution.Section 3 compares the performance of the proposed estimator with several estimators in the presence of outliers through a simulation study.Section 4 applies the proposed method for estimating the rate parameter of exponential distribution to real data sets.Finally,Section 5 concludes the paper.

2 Proposed Estimator

In this section,a brief explanation of the M-estimators is given followed by the introduction of the new estimator developed based on the probability integral transform statistical approach.Furthermore,the robustness of this new estimator was compared with the ordinary MLE and discussed in this section.

2.1 M-Estimators

M-estimators are generalized ML estimators that provide tools for measuring the robustness of the maximum likelihood-type estimates.As described in Huber [33],an estimatorβndefined by

or

is known as an M-estimator.ρis a measurable function onX×Θandψ(x,β)=(∂/∂β)ρ(x,vβ)denotes the derivative of the functionρwith respect toβ(when it exists).Note that ifρ(x,β)=-logf(x;β),thenis the ordinary MLE.Letx1,x2,...,xnbe a random sample from an exponential distribution as described in (2).The MLE for parameterγis given by

2.2 Probability Integral Transform Statistic Estimator

Assume thatx1,x2,...,xnis a random sample from the exponential distribution.Since the CDF of exponential distribution in Eq.(2) is continuous and strictly increasing,the random variablesF(x1),F(x2),...,F(xn)follow a standard uniform distribution,which isF(X)~Unif(0,1).Define

whereκ>0 denotes a tuning parameter to be used in regulating the balance between efficiency and robustness.Whenβ=γandκ=1,e-γxi=S(xi)=1-F(xi),is a random variable that follows the standard uniform distribution.Suppose thatu1,u2,...,unis a random sample from a standard uniform distribution.According to the strong law of large numbers,converges toasn→∞with probability 1.Thus,the probability integral transform statistic estimator (PITSE) of parameteris defined as the solution of the following equation:

Note that Eq.(10) can be solved using any numerical method such as bisection,secant,or Newton-Raphson method.

Lemma 1.For any fixedκ>0,the equationHn,κ(β)=1/(κ+1)in (10) has exactly one solution.

Proof of Lemma 1.Note thatHn,κ(β)is continuous on [0,∞) and thatHn,κ(0)=1>1/(κ+1)while limβ→∞Hn,κ(β)=0<1/(κ+1).Based on the intermediate value theorem,Hn,κ(β)=1/(κ+1)for someβin [0,∞).That is,Eq.(10) has at least one solution.Then,it can be shown thatfor allβ>0.Thus,Hn,κ(β)is strictly monotonic with respect toβand it can be concluded the Eq.(10) has a unique solution.

The PITSE for the rate parameter of exponential distribution is a class of M-estimator with

To investigate the properties of PITSE,the function introduced by Huber [33]can be utilized as follows:

Thus,based on Eq.(12),the following is obtained:

Theorem 1.For any fixedκ>0,the PITSEconverges toγasn→∞with probability 1.

Proof of Theorem 1.Based on Proposition 2.1 and Corollary 2.2 in Chapter 3 of Huber [33],λ(β)>0 forβ<γandλ(β)<0 forβ>γ.As the PITSE is uniquely defined,is a consistent estimator ofγ.

2.3 Efficiency of PITSE

The MLE is well known to be efficient in the sense that it has minimum asymptotic variance.For this reason,MLE is useful in providing a quantitative benchmark on the measure of efficiency.Note that the MLE for the rate parameter given in Eq.(8) is asymptotically normal with a meanγand varianceγ2/n,i.e.,MLE~N(γ,γ2/n).The asymptotic distributionPITSEis given in Theorem 2.

Theorem 2.For any fixedκ>0,the estimatoris asymptotically normal with mean 0 and varianceγ2(κ+1)2/(2κ+1),i.e.,~N(γ,γ2(κ+1)2/[n(2κ+1)]).

Proof of Theorem2.Based on Corollary 2.5 in Chapter 3 of Huber [33],λ′(γ)=-κ/[γ(κ+1)2]andTherefore,the asymptotic variance ofisσ20/[λ′(γ)]2=γ2(κ+1)2/(2κ+1).Note that for the finite sample,the variance ofisγ2(κ+1)2/[n(2κ+1)].

The ARE of PITSE is defined as the ratio of the asymptotic variance of the MLE to the asymptotic variance of the PITSE.In other words,the ARE measures the relative efficiency of the estimatorPITSEcompared withMLE.The ARE of PITSE for parameterγis given by

As the value ofκincreases,the value of ARE decreases.In other words,when the value ofκincreases,the PITSE gains robustness but with its relative efficiency reduced.By takingκclose to 0,ARE of PITSE can be made arbitrarily close to 1.The guidelines for choosing suitable ARE of PITSE in practical application are given in the simulation study section.

2.4 Breakdown Point of PITSE

The breakdown point (BP) is a useful measure of the robustness of a statistical approach in which the degree of sensitivity of an estimator to data contamination is measured.According to Huber [34],BP is defined as the largest contamination proportion that can be tolerated by an estimator before breaking down.A higher value of BP is an indication that an estimator is more robust against data contamination.In general,two types of BP exist;lower breakdown point(LBP) and upper breakdown point (UBP).In the present context ofγestimation,LBP is the largest proportion of lower contamination that can be tolerated by an estimator before forcing→∞.On the other hand,the UBP is the largest proportion of upper contamination that can be tolerated by an estimator before forcing→0.Since the contamination of the upper-end of the distribution is crucial in typical applications,only the UBP was emphasized in this study.

Theorem 3.The respective UBP and LBP of PITSE areκ/(κ+1)and 1/(κ+1).For the finite sample,the UBP and LBP of PITSE areand,respectively.

Proof of Theorem3.For any integer 1≤m≤n,the estimatoris defined asSuppose thatx1,x2,...,xmtake on values that approach ∞.Certainly,=1/(κ+1)<(n-m)/nbecause the terme-κxiPITSEis less than 1.This solution is valid if and only if 1/(κ+1)<(n-m)/n,that is,m<nκ/(κ+1).Thus,the finite sample UBP is,whereis the ceiling function.By takingn→∞in the finite sample UBP,UBP equal toκ/(κ+1)will be obtained.Similarly,letx1,x2,...,xmtake on values approaching 0.it can be shown that=1/(κ+1).This solution is finite if and only if 1/(κ+1)>m/n,that is,m<n/(κ+1).Therefore,the finite sample LBP isand by lettingn→∞,the LBP will be equal to 1/(κ+1).

2.5 Gross Error Sensitivity of PITSE

The breakdown point (BP) is a useful measure of the robustness of a statistical approach in which the degree of sensitivity of Gross error sensitivity (GES) is also an important measure of the robustness of an estimator.According to Hampel et al.[35],GES is the supremum of the absolute value of the influence function (IF).The IF of an estimatorβnthat satisfies Eq.(7) can be defined by the following:

whereβ(F) denotes the solutionβnof Eq.(7) with samples generated from the CDFF.Then,the GES is given by

An estimator with small GES should be more robust than that with larger GES.The MLE,which is,has an infinite GES becauseψ(x;)is an unbounded function.This result proves that MLE is a non-robust estimator.Then,the GES of PITSE is provided in Theorem 4.

Theorem 4.The GES of the PITSE is max{γ(κ+1)/κ,γ(κ+1)}.

Proof of Theorem 4.By applying the formula of IF in Eq.(15),the IF of PITSE can be demonstrated asIF(x;ψ,F)=-ψ(x;γ)γ(κ+1)2/κ.Asψ(x;γ) is monotone inx,the maximum value occurs when eitherx→∞orx→0.Note thatψ(x;γ)=-1/(κ+1)whenx→∞andψ(x;γ)=κ/(κ+1)whenx→0.Then,using the formula in Eq.(16),the GES of PITSE is obtained.

Tab.1 shows the properties of the PITSE based on ARE,BP,and GES for different values of tuning parameterκ.The properties of the PITSE can be summarized as follow:

· As the value ofκincreases,the PITSE loses efficiency (ARE decreases) but also becomes more resistant to upper contamination.

· As the value ofκincreases,the UBP increases,and the LBP,which suggests that the robustness of the PITSE,increases against upper contamination.On the other hand,the LBP decreases as the value ofκincreases,which indicates that the robustness of the PITSE decreases against lower contamination.

· As the value ofκincreases up to 1,the GES of the PITSE decreases.But then,the GES of the PITSE started to increase when the value ofκis greater than 1.

Table 1:The ARE,BP,and GES of PITSE for different values of tuning parameter κ

3 Simulation Study

In this section,a simulation study to compare the performance of the MLE,PITSE,WLE,and M-scale estimator in the presence of outliers is conducted.The design and results of the simulation study are given in the next two following subsections.Then,the guidelines for selecting the suitable ARE of the PITSE for practical application purposes are provided.

3.1 Design of Simulation

In the simulation study,the methods considered for comparison were MLE,PITSE (90% and 70% ARE),WLE (90% and 70% ARE),and M-scale estimator (70% ARE).The data sets were simulated from an exponential distribution,Exp(γ),with different values of rate parameter,γ=0.5,1,1.5,and 2.The sample size used in this simulation study was divided into two cases:small sample sizes (n=30,50,and 70) and large sample sizes (n=100,300,and 500).Then,some observations in the data sets were randomly selected and replaced with outliers.According to Lin et al.[36],the observations are from exponential distribution with parameterωγ,which is Exp(ωγ) with 0<ω<1 considered as outliers.Thus,in this simulation study,the outliers were generated from an exponential distribution with parameter 0.05γ,which was Exp(0.05γ).For the case of small sample sizes,the outliers were generated from Exp(0.05γ) for several fixed numbers,k=1,3,and 5.For the case of large sample sizes,the outliers were simulated from Exp(0.05γ)for several fixed proportions,ε=1%,5%,and 10%.This simulation was repeated for 10,000 data sets found based on 10,000 simulation runs.

The performance of the estimators was assessed in terms of percentage relative root mean square error (RRMSE).For a given true value of the rate parameter of exponential distribution,γ,theRRMSEis given by

3.2 Results of Simulation

Results based on the simulation study are presented in Tabs.2-7.For the case of small sample sizes,whenn=30 andk=1,for all values of the parameterγ,the value ofRRMSEfound for PITSE and WLE with 90% ARE were nearly similar and smaller than those for other estimators as shown in Tab.2,indicating that these two estimators performed better than other estimators for this particular case.Furthermore,from Tab.2,when the number of outliers increased tok=3,the PITSE,WLE,and M-scale estimator with 70% ARE performed almost equally and outperformed other estimators due to smaller values ofRRMSE.Whenk=5,it was observed that the PITSE and M-scale estimator with 70% ARE performed almost equally and outperformed other estimators as presented in Tab.2.As the sample size increases,the performance of all estimators improves for all cases considered because the values ofRRMSEbecame smaller as shown in Tabs.2-4.For all cases of small sample sizes,the MLE produced a high value ofRRMSE,indicating a poor performance of this method when the outliers were present in the data set.

Table 2:Results of RRMSE for estimation of parameter γ with n=30 and k=1,3,5

The simulation results for the case of large sample sizes are shown in Tabs.5-7.Based on Tab.5,forn=100,the PISTE and WLE with 90% and 70% ARE performed almost equally and outperformed other estimators in the case of a small degree of contamination,where∈=1%.When the degree of contamination increased to moderate,which was∈=5% the PITSE,WLE,and M-scale estimator with 70% ARE produced much better performance than other estimators.At a high degree of contamination where∈=10%,the PITSE and M-scale estimator with 70%ARE performed almost equally and outperformed other estimators.Although WLE is a robust estimator,it performed substantially worse than the PITSE and the M-scale estimator in the case of a high degree of contamination.The reason is that a large proportion of observations (outliers)in the data set was not being considered,which reduced the sample size and,consequently causing the efficiency of the WLE to be decreased.As shown in Tabs.5-7,when the sample size increases,the performance of all estimators improves.However,similar to the case of small sample size,the performance of MLE was also poor when outliers were present in the data set for the case of large sample size.Thus,the use of the MLE for estimating the parameterγshould be avoided since this estimator does not protect outliers even in the small degree of contamination.

Table 3:Results of RRMSE for estimation of parameter γ with n=50 and k=1,3,5

Table 4:Results of RRMSE for estimation of parameter γ with n=70 and k=1,3,5

Table 5:Results of RRMSE for estimation of parameter γ with n=100 and ∈=1%,5%,10%

Table 6:Results of RRMSE for estimation of parameter γ with n=300 and ε=1%,5%,10%

3.3 Guidelines for Selecting Suitable ARE of PITSE in Practical Application

In practice,the suitable ARE of PITSE can be selected based on sample size and number or proportion of outliers.Based on a comprehensive simulation study,the guidelines for selecting the suitable ARE of PITSE in practical application are provided in Tab.8.

It should be noted that when there are a large number or proportion of outliers in a data set,PITSE with ARE not less than 60% should be used to make sure that the PITSE has a reasonable efficiency in estimating the rate parameter of exponential distribution.

Table 7:Results of RRMSE for estimation of parameter γ with n=500 and ∈=1%,5%,10%

Table 8:Guidelines for selecting suitable ARE of PITSE in practical application

4 Application on Reliability Data Sets

In this section,three applications of exponential distribution are proposed using three real data sets to compare the performance of MLE,PITSE,WLE,and M-scale estimator.This comparative study considered several ARE levels for PITSE and WLE,while a fixed ARE level of 70% was used for the M-scale estimator.The first data set (Set 1) was obtained from Linhart et al.[37],which represents the failure times of the air-conditioning system of an airplane.The second data set (Set 2) consisted of the reliability data of a 180-tonne rear dump truck obtained from Badr et al.[38].Lastly,the third data set (Set3) was obtained from Proschan [39]consisting of the number of successive failures for the air conditioning system of each member in a fleet of 13 Boeing 720 jet airplanes.The pooled data,yielding a total of 213 observations,have also been analyzed by other authors such as Ku¸s [40]and Tahmasbi et al.[8].All the data sets used are given in Appendix A.

Figure 1:Generalized boxplot for data (a) Set 1,(b) Set 2 and (c) Set 3

Table 9:Descriptive statistics for data sets 1,2 and 3

Table 10:Parameter estimates and goodness of fit of the exponential distributions of data sets 1,2 and 3

To identify the presence of outliers in all data sets,the generalized boxplot method [41],which is suitable for skewed and/or heavy-tailed distributions,was utilized.According to Bruffaerts et al.[41],this boxplot relies on a simple rank-preserving transformation that allows the data to fit a Tukey g-and-h distribution.It was found that the generalized boxplot was robust against outliers and has clear advantages over the standard boxplot particularly for the case of skewed and/or heavy-tailed distributions.Fig.1 shows the generalized boxplot for all considered data sets,while Tab.9 presents the descriptive statistics of these data sets.

Figure 2:Fitted exponential densities on the histograms of data (a) Set 1,(b) Set 2 and (c) Set 3 based on several different estimators

To compare the performance of all considered methods in estimating the parameter of exponential distribution,the Kolmogorov-Smirnov (KS) tests were employed as a goodness-of-fit assessment.The best method was determined by choosing the smallest values of KS statistics as well as the highestp-values of the KS test.Tab.10 reports the estimated parameter and goodnessof-fit of the exponential model for all data sets.For data Set 1,it can be observed that the PITSE (70% ARE) and M-scale estimator (70% ARE) performed almost equally and provided a better estimation of exponential parameter compared to other methods based on their smallest KS statistic and highestp-value of the KS test.For data Set 2,the PITSE has outperformed other methods in estimating the rate parameter.It can be also seen for data Set 2 that the MLE and M-scale estimators failed to provide a good estimation of rate parameters based on the small value of the KS test found for these two estimators.For data Set 3,the PITSE was found to be the best method to estimate the rate parameter.Nevertheless,other methods were also found quite reliable for estimating the rate parameter based on the reasonablep-value of KS test (p-value>0.05).Fig.2 demonstrates the fitted exponential density on the histogram of each data set based on several different estimators.

Since the measures of reliability based on the exponential distribution depend on the parameterγ,it is crucial to employ the appropriate method for estimating parameterγ.Based on the simulation study and its application to real data sets,it was demonstrated that the proposed PITSE can represent a viable alternative for estimating the rate parameter of the exponential distribution,particularly for the case when outliers are present in the data.

5 Conclusion and Discussion

In this study,a robust and efficient estimator for the parameter of exponential distribution called PITSE has been introduced based on probability integral transform statistic.The asymptotic variance,BP,and GES were derived to study the PITSE properties.The advantage of PITSE is that it is conceptually simple and easy to compute.According to the simulation study,PITSE performed better than MLE and was comparable with WLE and M-scale estimator when outliers are present in the data set.However,in the case of a high degree of contamination,the performance of WLE was worse than PITSE and M-scale estimator.On the other hand,the M-scale estimator only has the maximum ARE of about 71%,which makes it unsuitable for estimating the parameter of exponential distribution in the case of a small degree of contamination.Therefore,the PITSE introduced in this study was considered a viable alternative for estimating the parameter of exponential distribution in the presence of outliers.Finally,the application on three real data sets showed that the PITSE provided desirable protection against outliers.The R commands for PITSE are available in Appendix B.

Although the PITSE proposed in this study was considered viable for estimating the exponential parameter,there existed a limitation regarding the scope of the current work.It was the reliability modeling under exponential model assumptions that was limited to certain cases(constant failure rate) since the exponential distribution is a simple model consisting of one parameter.There are two-parameter distributions such as Weibull that could provide a better fit to the reliability data,hence providing a better reliability estimation.Therefore,for future work,the robust and efficient estimator for the Weibull parameters can be developed based on the probability integral transform statistical approach.From there,it is believed that a better reliability estimation can be obtained particularly for the case when the outliers are present in the data set.

Acknowledgement:The authors would like to thank the editor for the time spent reviewing this manuscript.

Funding Statement:This work is supported by the Universiti Kebangsaan Malaysia [Grant Number DIP-2018-038].

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A.Data Sets 1,2 and 3

Set 1:

23,261,87,7,120,14,62,47,225,71,246,21,42,20,5,12,120,11,3,14,71,11,14,11,16,90,1,16,52,95

Set 2:

0.01,0.01,0.01,0.01,0.01,0.02,0.02,0.02,0.02,0.03,0.04,0.06,0.08,0.1,0.1,0.12,0.12,0.12,0.13,0.14,0.15,0.15,0.15,0.16,0.16,0.17,0.18,0.18,0.19,0.2,0.21,0.22,0.23,0.25,0.26,0.28,0.28,0.3,0.32,0.34,0.36,0.38,0.39,0.41,0.41,0.42,0.43,0.44,0.44,0.45,0.45,0.5,0.53,0.56,0.58,0.58,0.61,0.62,0.62,0.62,0.64,0.66,0.7,0.7,0.7,0.72,0.77,0.78,0.78,0.8,0.82,0.83,0.85,0.86,0.96,0.97,0.98,0.99,1.05,1.06,1.07,1.18,1.35,1.36,1.42,1.55,1.59,1.65,1.73,1.77,1.79,1.8,1.91,2.09,2.14,2.15,2.15,2.31,2.33,2.36,2.43,2.45,2.5,2.51,2.58,2.64,2.68,3.08,3.94,4.12,4.33,4.42,4.53,4.88,4.97,5.11,5.32,5.55,6.63,6.89,7.62,11.41,11.76,11.85,12.36,13.22

Set 3:

194,15,41,29,33,181,413,14,58,37,100,65,9,169,447,184,36,201,118,34,31,18,18,67,57,62,7,22,34,90,10,60,186,61,49,14,24,56,20,79,84,44,59,29,118,25,156,310,76,26,44,23,62,130,208,70,101,208,74,57,48,29,502,12,70,21,29,386,59,27,153,26,326,55,320,56,104,220,239,47,246,176,182,33,15,104,35,23,261,87,7,120,14,62,47,225,71,246,21,42,20,5,12,120,11,3,14,71,11,14,11,16,90,1,16,52,95,97,51,11,4,141,18,142,68,77,80,1,16,106,206,82,54,31,216,46,111,39,63,18,191,18,163,24,50,44,102,72,22,39,3,15,197,188,79,88,46,5,5,36,22,139,210,97,30,23,13,14,359,9,12,270,603,3,104,2,438,50,254,5,283,35,12,130,493,100,7,98,5,85,91,43,230,3,130,230,66,61,34,487,18,14,57,54,32,67,59,134,152,27,14,102,209

Appendix B.R Commands for PITSE

### PITSE for the rate parameter ###

f<-function(data,l,k){

n<-length(data)

fx<-(sum(exp(-l * k * data))/n)-(1/(k+1))

return(fx)

}

#solve using secant method

# k-tuning parameter

# l1-1stinitial value;l2-2ndinitial value

pitse<-function(data,k,l1,l2,num=1000,eps=1e-05,eps1=1e-05)

{

i=0

while ((abs(l1-l2)>eps) &&(i<num)) {

c=l2-f(data,l2,k) * (l2-l1)/(f(data,l2,k)-f(data,l1,k))

l1=l2

l2=c

i=i+1

}

rate<-l2

return(rate)

}