Smooth Tests of Fit for the Lindley Distribution

We consider the little-known one parameter Lindley distribution. This distribution may be of interest as it appears to be more flexible than the exponential distribution, the Lindley fitting more data than the exponential. We give smooth tests of fit for this distribution. The smooth test for the Lindley has power comparable with the Anderson-Darling test. Advantages of the smooth test are discussed. Examples that illustrate the flexibility of this distributions is given.


Introduction
The Lindley distribution was introduced by Lindley [1] to analyze failure time data.The Lindley distribution belongs to the exponential family and it can be written as a mixture of the exponential and gamma distributions.Its motivation arises from its ability to model failure data with increasing, decreasing, unimodal, and bathtub-shaped hazard rates.See [2].Ghitany et al. [3] give many properties of the Lindley distribution.They suggest it is often a better model than the traditional exponential distribution that is commonly used to model lifetime or waiting time data.For example, the exponential distribution is limited to lifetime data with coefficients of variation 1, skewness 2, and kurtosis 6.The Lindley distribution allows a greater range for these coefficients: namely 1/ √ 2 to 1, √ 2 to 2 and 6 to 9 respectively.Similarly, the hazard rate function of the exponential is more limited than that of the Lindley.
Ghitany et al. [3] examine the fit of the Lindley distribution to some waiting time data by looking at plots and by showing the Lindley likelihood is better than the exponential likelihood.However, this does not demonstrate that the Lindley distribution fits the data well, only that it fits better than the exponential.Assessment of the plots is subjective and here we derive a smooth test of fit to give a more objective assessment of goodness of fit of the Lindley model.We also consider the widely used Anderson-Darling test.
The Lindley distribution has probability density function (1 + x) e −θx for x > and zero otherwise, in which θ > 0 and cumulative distribution function x > and zero otherwise.
Smooth tests of fit can be found using the second and third order smooth test components.Section 2 discusses smooth testing for goodness of fit.Section 3 looks at the approach to the asymptotic chi-squared distributions of the smooth test statistics and finds them to be quite slow.It is suggested that p-values be found using the parametric bootstrap.A slightly expanded version of an algorithm in Ghitany et al. [3] generating random Lindley variates is given in Section 3 so that these p-values can be calculated.Some powers of the smooth test and the Anderson-Darling test are compared in Section 4 while Section 5 gives two examples.

The Smooth Test Statistics
Tests of fit are an important element in determining the suitability of statistical models.The genesis of the smooth tests is Neyman [4], who nested the null probability density function in a family of distributions indexed by the elements of a vector parameter and tested if that parameter was consistent with zero.Neyman only considered testing for distributions with no nuisance parameters and hence, after transforming by the probability integral transformation, testing for the simple uniform (0, 1) distribution.Best and Rayner [5] give an early discussion of smooth tests for distributions with nuisance parameters by considering a smooth test for normality.The subject has undergone considerable development, especially in recent years.See, for example, [6,7].
Of particular interest here is [6] (Section 9.2) in which the following is shown.Suppose we have a random sample X 1 , . . ., X n from a distribution hypothesized to have probability density function f (x; β) in which β = (β 1 , . . ., β q ) T is a q × 1 vector of nuisance parameters.Moreover, assume that the densities f (x; β) form a multi-parameter exponential family.An alternative probability density function is taken to be in which θ = (θ 1 , . . ., θ k ) T , it is assumed a normalizing constant C(θ, β) exists and {h i (x; θ, β)} is a set of orthonormal functions with weight function f (x; β).
When testing H: θ = 0 against K: θ = 0 the smooth test statistic is a sum of squares of the components V q+1 , . . ., V k , that is, √ n, β being the maximum likelihood estimator of β when θ = 0.The model g k (x; θ, β) is over-parametrized, with the θ 1 , . . ., θ q playing the same role as β 1 , . . ., β q .One way of dealing with this technically is to drop θ 1 , . . ., θ q from g k (x; θ, β).The over-parametrization is apparent when it is realized that the likelihood equations are equivalent to ∑ n j=1 h r (X j ; β)/ √ n = 0 for r = 1, . . ., q and so V 1 ≡ . . .≡ V q ≡ 0. The non-degenerate components are asymptotically independent and asymptotically χ 2  1 distributed.Since the components V r involve the orthonormal polynomials, we now give the orthonormal polynomials of a random variable X up to order three.These results are true for any distribution for which the moments exist.The notation suppresses the dependence on the nuisance parameter.We take h 0 (x) = 1 for all x.Write µ for the mean of X and µ r , r = 2, 3, . . .for the central moments of X.Then and Should further orthonormal polynomials be required the recurrence relation in [8] can be used.See [9] (Appendix D) for details of the numerical construction of orthogonal polynomials in R.
For the exponential distribution the moments are relatively simple and lead to The early literature suggested that the components V r were diagnostic.That is, if the test statistic V 2 q+1 + . . .+ V 2 k was found to be significant at some prescribed level then the V r could be used to diagnose the model failure.For example, a significant V 3 would indicate the third moment of the hypothesized distribution was not consistent with the data.In fact, a significant V r could be attributed to moments up to the 2r th .However, in most applications it is reasonable to say that a significant V r suggests model failure in moments up to the r th .This could be confirmed (or not) by plotting the data.
For distributions from a one parameter exponential family the maximum likelihood and method of moments estimators coincide.It follows that θ can be readily found by solving The choice of the order k of the smooth test is typically a trade-off between the test power and the alternatives detected.A larger k means the test is more omnibus, seeking to detect alternatives to the null in a richer family of distributions.A smaller k gives a more focused test with greater power for the alternatives for which the test is sensitive, but no power for other alternatives.Based solely on intuition Neyman [4] took k to be four.Rayner et al. [6] discuss modeling approaches for the choice of k.A larger k requires more orthonormal polynomials and hence more moments to calculate them.Here we make the pragmatic option of considering only two non-degenerate components through V 2 2 , V 2  3 and S = V 2 2 + V 2 3 .In the next section we briefly consider the approach of these test statistics to their asymptotic chi-squared distributions.Use of S is similar to the test in [10] which is also based on the first two non-zero smooth test components.

The Approach to the Asymptotic Distribution
In Table 1 we look, for θ = 0.5 and 1.5, at the approach to the asymptotic χ 2 1 distribution of V 2 2 and V  and the approach to the asymptotic χ 2 2 distribution of S = V 2 2 + V 2 3 .The results in Table 1 are 5% critical values found using 100,000 simulations of Monte Carlo samples of size n.A random variate generator, given below, is needed for these results, the powers of the next section and bootstrap p-values.
Table 1.5% critical values based on 100,000 simulations of samples of size n for V 2 2 , V 2 3 and S when θ = 0.5 and 1.5.The convergence to the asymptotic values is quite slow and so we suggest in applications that for smaller sample sizes the parametric bootstrap will be needed to find p-values.The Table 1 results are similar for θ = 0.5 and θ = 1.5.
To find parametric bootstrap p-values generate Lindley samples of size n many (10,000 subsequently) times and take the p-values as the proportion of the samples with test statistics greater than or equal to the values of the test statistics for the original dataset.To find parametric bootstrap powers for each alternative distribution generate many (10,000 subsequently) samples of size n.Then the power is the proportion of these samples with p-value less than the significance level.

Power Comparisons
For a significance level of α = 0.05 and a sample size of n = 20, Tables 2-4 give some powers for tests based on V 2 2 , V 2 3 , S and AD where AD is the Anderson-Darling test statistic in which {z (i) } are ordered values of {z i }, where z i = F(x i ; θ).The Anderson-Darling test is included as it is well-known and usually performs well in power studies for other distributions.The tests based on S and AD have similar powers.Notice that we used critical values such that all four of the tests we examined had equal test size.This is necessary so that power advantages are not due to poor approximation to the finite sample null distribution by, say, use of the asymptotic critical values.Random values from the alternative distributions were generated using the IMSL software package.Alternatives were chosen to give powers not all near 0.05 or 1.00, to demonstrate differences in the effectiveness of the tests.
The test based on V 2 2 has slightly less power than the tests based on S and AD.It has little power if the alternative has similar variance to the Lindley variance, which is quite reasonable as it is particularly sensitive to distributions with the Lindley variance.The test based on V 2 3 is, roughly, testing for distributions with the Lindley skewness when V 2 2 is not significant.This is why it is useful to apply V 2 2 and V 2 3 together, either separately as in exploratory data analysis, or more formally together, via S.

Examples
In both the examples following the bootstrap p-values are based on 1000 samples.Waiting Time Data.Ghitany [3] gives the waiting times (in minutes) before service of 100 bank customers.On the basis of a superior log likelihood they conclude that the Lindley distribution gives a better fit than the exponential.
In testing for the Lindley distribution, we find bootstrap p-values for the tests based on V 2 2 , V 2 3 , S and AD to be 0.61, 0.49, 0.70 and 0.50 respectively.Using the asymptotic chi-squared distributions the p-values for V 2 2 , V 2 3 and S are 0.62, 0.60 and 0.77 respectively.We can conclude that for a sample size this large there is acceptable agreement between the p-values based on the bootstrap and the asymptotic chi-squared distributions.Moreover, the Lindley distribution fits the data well.
Operational Lifetime Data.Angus [11]  In testing for the Lindley distribution bootstrap p-values for the tests based on V 2 2 , V 2 3 , S and AD are 0.14, 0.27, 0.17 and 0.30 respectively.Using the asymptotic chi-squared distributions the p-values for V 2 2 , V 2 3 and S are 0.21, 0.44 and 0.34 respectively.For this sample size there is not good agreement between the p-values based on the bootstrap and the asymptotic chi-squared distributions.Moreover, it appears that, unlike the exponential distribution, the Lindley fits the data well.

Conclusions
Tests of fit are an important element in determining the suitability of statistical models.They give objective model assessments, which may be complemented by subjective graphical methods such as Q-Q plots.Comparison of statistics such as likelihoods can show that one model is a better fit than another, but not whether either model is a good fit or not.We have given a smooth test of fit statistic S for the Lindley distribution.Two datasets illustrate the applicability of the Lindley distribution, which provides a good model for both datasets, while each of the exponential models fails.For small sample sizes, we suggest that p-values be given using the parametric bootstrap.For larger sample sizes, the asymptotic χ 2 distribution may be useful.

Table 3 .
Powers of tests based on V 2 2 , V 2 3 , S and AD for α = 0.05 and n = 50.
[6]se data are analysed in[6](Example 6.3.1 and Example 11.7.1)where it is found that a test for exponentiality is significant at the 0.05 level while several tests for the generalised Pareto distribution are not significant at the 0.05 level.