Exact Statistical Inference for Nonhomogeneous Poisson Processes

Bo Henry Lindqvist

bo@math.ntnu.no

Department of Mathematical Sciences Norwegian University of Science and Technology N-7491 Trondheim

1. Introduction and Summary

Nonhomogenous Poisson processes (NHPP) are widely used as models for events occuring in time, for example failures of a repairable system. In practice it may be of interest to check the NHPP property by statistical tests given failure data. The main purpose of this paper is to develop exact tests for null hypotheses that failure times follow NHPPs of particular parametric forms.

One way of testing NHPP models is to embed the models in more general parametric failure models and use likelihood ratio tests based on asymptotic chi-square distributions. A problem is here that the observed number of failures is often too small to justify the use of asymptotic distributions. This motivates the use of exact statistical inference whenever possible (Baker 1996, Gaudoin 1999).

In the present paper we consider exact (Monte Carlo) tests obtained by conditioning on a sufficient statistic under the NHPP model of the null hypothesis. A new method for simulating conditionally on sufficient statistics (Engen and Lillegård 1997, Lindqvist and Taraldsen 2000) will be used. A closely related way of obtaining exact confidence intervals in parametric models (Lillegård and Engen 1999) is also briefly considered.

2. Observations, Likelihood Functions and Sufficient Statistics

We assume that a repairable system is observed from time t = 0 and until n failures have occurred. This is known as failure truncation. Most of our results may be modified to the case of time truncation, where the process is observed until a given time t_{0}.

Suppose we model the failure process by an NHPP with intensity function λ(t) and failure times T = (T_{1}, …, T_{n}). The log likelihood function resulting from observed failure times {t_{j}} is then given as ∑_{j=1}^{n} log λ(t_{j}) − ∫_{0}^{t_{n}} λ(u)du (Crowder et al. 1991).

Two popular parametrizations of NHPPs are the power law, with intensity function given by λ_{pow}(t) = abt^{b−1} and the log linear law, with intensity function λ_{log}(t) = exp(a + bt). Substituting each of these functions into the above log likelihood and using the factorization criterion for sufficiency, gives that the statistic S = (log T_{n}, ∑_{j=1}^{n−1} log T_{j}) is sufficient in the power law case, while S = (T_{n}, ∑_{j=1}^{n−1} T_{j}) is sufficient in the log linear case.

3. Conditional Testing Given a Sufficient Statistic

Let T be the vector of failure times as described in the previous section, and consider testing of the null hypothesis H_{0} that these data come from an NHPP of a particular parametric type. Let W ≡ W (T) be any test statistic for revealing departure from the null model. Suppose that S ≡ S(T) is a sufficient statistic for the unknown parameters under the null hypothesis. Then if w_{obs} is the observed value of W we can obtain a conditional p-value by computing the conditional probability P_{H_{0}}(W ≤ w_{obs}|S = s). For this we need to know the conditional distribution of W given S = s under H_{0}, which by sufficiency is independent of the unknown parameters and thus in principle can be found. This may, however, be difficult in many practical cases and we shall thus rely on simulations. Our tests will thus be Monte Carlo tests.

4. Monte Carlo Conditioning on a Sufficient Statistic

In general, let the model for the observation T (vector) under the null hypothesis be specified in terms of a parameter (vector) θ. Suppose that for a given value of θ we can simulate realizations of T by T = χ(U, θ) for some function χ and a random vector U with known distribution. Further, suppose S(T) ≡ S is sufficient for θ. Then S can be simulated by the function τ(U, θ) ≡ S(χ(U, θ)).

Consider now computation of conditional expectations of the form E(ϕ(T)|S = s) for given functions ϕ, where s is the observed value of S. By sufficiency, this conditional expectation is independent of θ and the clue of the approach is that it can be expressed in terms of ordinary expectations of functions of U (Theorems 1 and 2 below). The basic idea comes from Engen and Lillegård (1997) who, however, erroneously claimed that the conclusion of Theorem 2 will also hold under the more general assumptions of Theorem 1.

Theorem 1 (Lindqvist and Taraldsen, 2000) Suppose θ and S(T) take values in R^{k} and suppose the equation τ(u, θ) = s has the unique solution θ = θ̂(u, s) for each fixed u, s. Let f(θ) be a nonnegative function defined on the parameter space, and let det ∂_{θ}τ(u, θ) be the determinant of the matrix of partial derivatives of τ(u, θ) for fixed u. Then

E(ϕ(T)|S = s) = E_{U} [ϕ(χ(U, θ̂(U, s)))|\frac{f(θ)}{det ∂_{θ}τ(U,θ)}|_{θ=θ̂(U,s)}] ; E_{U} [|\frac{f(θ)}{det ∂_{θ}τ(U,θ)}|_{θ=θ̂(U,s)}]

It is tacitly assumed above that f is such that the expectations exist and the denominator is positive, but f may otherwise be arbitrarily chosen. The idea is that the expectations can be computed by simulation by drawing an i.i.d. sample of U and averaging the expressions inside the expectations. It follows that if a function f(·) can be chosen so that |f(θ)/ det ∂_{θ}τ(U, θ)|_{θ=θ̂(U,s)} does not depend on U,

then we will have E(ϕ(T)|S = s) = E_{U} [ϕ(χ(U, θ̂(U, s)))] which means that the function χ(U, θ̂(U, s)) can be used to sample directly from the conditional distribution of T given S = s. Unfortunately, it is not always possible to find such an f, but the following sufficient condition can be given:

Theorem 2 (Lindqvist and Taraldsen, 2000) Let the situation be as in the previous theorem. Assume that there exist functions r and τ̃ with τ(u, θ) = τ̃(r(u), θ), such that the equation τ̃(v, θ) = s has a unique solution v = ṽ(θ, s) for all (θ, s). Then χ(U, θ̂(U, s)) is distributed as the conditional distribution of T given S = s.

The new assumption of Theorem 2 means that τ(u, θ) depends on u only through r(u), which usually has a much lower dimension than u, and has the property that for given θ, r(u) is uniquely determined by s. Note that ṽ(θ, S) is a pivotal quantity in the classical meaning.

5. Conditional Simulation for Parametric NHPP Models

The first n events of an NHPP with intensity function λ(·) can be simulated by letting U = (U_{1}, U_{2}, …, U_{n}) be the first n events of a homogeneous Poisson process with unit intensity, and then letting

T_{j} = Λ^{−1}(U_{j}) ; j = 1, …, n  (1)

where Λ^{−1} is the inverse function of the cumulative intensity function Λ(t) ≡ ∫_{0}^{t} λ(u)du. 5.1 Power Law Intensity From Section 2 follows that T_{j} = (U_{j}/a)^{1/b}. With notation from the previous section we simulate T by

χ(u; a, b) = ((u_{1}/a)^{1/b}, …, (u_{n}/a)^{1/b})

while the sufficient statistic S = (log T_{n}, ∑_{j=1}^{n−1} log T_{j}) is simulated by

τ(u; a, b) = ((log u_{n} − log a)/b, (∑_{j=1}^{n−1} log u_{j} − (n − 1) log a)/b)

Thus the pivotal condition in Theorem 2 holds with r(u) = (log u_{n}, ∑_{j=1}^{n−1} log u_{j}). Letting the observed times be t = (t_{1}, …, t_{n}) and letting s = (t_{n}, ∑_{j=1}^{n−1} t_{j}), it is straightforward to obtain θ̂(u, s) ≡ (â(u, s), b̂(u, s)) by solving τ(u; a, b) = s for a and b.

Samples t̃ = (t̃_{1}, …, t̃_{n}) from the conditional distribution of T given S = s can now be obtained by first sampling u = (u_{1}, …, u_{n}) and then computing t̃= χ(u; θ̂(u, s)). We get

t̃_{j} = (u_{j}/u_{n})^{1/b̂(u,s)}t_{n} ; j = 1, …, n

Note that t̃ can be easily simulated by noting that the u_{j}/u_{n} are distributed as the order statistic of n − 1 independent uniforms on [0, 1]. 5.2 Log Linear Intensity Now (1) becomes T_{j} = log(1 + be^{−a}U_{j})/b, so T is simulated by

χ(u; a, b) = (log(1 + be^{−a}u_{1})/b, …, log(1 + be^{−a}u_{n})/b)

and S = (T_{n}, ∑_{j=1}^{n−1} T_{j}) by

τ(u; a, b) = (log(1 + be^{−a}u_{n})/b, ∑_{j=1}^{n−1} log(1 + be^{−a}u_{j})/b)

It is clear that the pivotal condition of Theorem 2 is not satisfied here. This should be no surprise since it is well known that the log linear NHPP model has no interesting pivotal statistic. In order to compute p-values of conditional tests, we thus have to use Theorem 1 with some arbitrarily chosen function f(a, b).

Now we can first find θ̂(u, s) ≡ (â(u, s), b̂(u, s)) by solving for a, b the equations

t_{n} = \frac{1}{b} log(1 + be^{−a}u_{n})  (2)
∑_{j=1}^{n−1} t_{j} = \frac{1}{b} ∑_{j=1}^{n−1} log(1 + be^{−a}u_{j})  (3)

Solving (2) for a and substituting into (3) gives

b ∑_{j=1}^{n−1} t_{j} = ∑_{j=1}^{n−1} log [1 + \frac{u_{j}}{u_{n}} (e^{t_{n}b} − 1)]  (4)

which is an equation in b only. By differentiating twice it is seen that the right hand side of (4) is convex in b. Further consideration leads to the conclusion that the equation (4) in addition to the trivial solution b = 0 has a unique additional solution, which is the one that solves our problem, and which is easily found by numerical methods. The solution for a is then finally found from (2).

Now t̃ = χ(u; â(u, s), b̂(u, s)) is given by t̃_{j} = log (1 + (e^{t_{n}b̂(u,s)} − 1)(u_{j}/u_{n})) /b̂(u, s) and a straight-

forward computation shows that the determinant det∂_{a,b}τ(u; a, b) with (â(u, s), b̂(u, s)) substituted for (a, b) is given by

(t_{n} ∑_{j=1}^{n−1} h_{j} − h_{n} ∑_{j=1}^{n−1} t_{j}) /b̂(u, s)^{2}

where h_{j} = ((e^{t_{n}b̂} − 1)(u_{j}/u_{n}))/(1 + (e^{t_{n}b̂} − 1)(u_{j}/u_{n})) and b̂ ≡ b̂(u, s). As for the power law case we simulate the u_{j}/u_{n} for j = 1, …, n − 1 as the order statistic from a set of n − 1 i.i.d. uniforms on [0, 1].

6. Statistical Inference in NHPP Models

6.1 Goodness-of-Fit Testing

The identity (1) is equivalent to U_{j} = Λ(T_{j}). It follows that if Λ(·) is the cumulative intensity function of the NHPP T_{1}, T_{2}, …, then Λ(T_{1}), Λ(T_{2}), … is a homogeneous Poisson process with unit intensity. Thus the transformed times V_{j} = Λ(T_{j})/Λ(T_{n}) for j = 1, …, n − 1 are distributed as the order statistic of n − 1 i.i.d. uniform variables on [0, 1]. If Λ^{∗}(·) is an estimate of Λ(·) based on data t = (t_{1}, …, t_{n}), then we shall define estimated transformed times v_{1}^{∗}, …, v_{n−1}^{∗} by v_{j}^{∗} = Λ^{∗}(t_{j})/Λ^{∗}(t_{n}). One then anticipates these to behave much similar to uniform variables, and goodness-of-fit testing may thus be based on comparing the behaviour of the estimated transformed times to that of uniform variates.

Baker (1996) shows for the power law process that when Λ^{∗} is based on the maximum likelihood estimates for the parameters, then the estimated transformed times are pivots, i.e. have distributions which do not depend on the unknown parameters. This follows in fact from the representation t_{j} = (u_{j}/a)^{1/b} in Section 5.1, noting that b^{∗} = −n/ ∑_{j=1}^{n−1} log(t_{j}/t_{n}) is the maximum likelihood estimate of

b based on the data (Crowder et al. 1991). Then we have v_{j}^{∗} = (t_{j}/t_{n})^{b^{∗}} = (u_{j}/u_{n})^{−n/ ∑_{j=1}^{n−1} log(u_{j}/u_{n})} which is independent of the parameters.

Baker (1996) derives a class of score tests based on the estimated transformed times. A special case, which we shall use for illustration, is W = ∑_{j=1}^{n}(v_{j}^{∗} − v_{j−1}^{∗})^{2} where v_{0}^{∗} = 0, v_{n}^{∗} = 1. The null hypothesis of NHPP is rejected for either too small or too large values of this statistic. Note that since the v_{j}^{∗} are pivots, we can in the power law case compute (by simulation) the unconditional probabilities P_{H_{0}}(W ≤ w_{obs}).

In the case of log linear intensity (Section 5.2) we get v_{j}^{∗} = (e^{b^{∗}t_{j}} −1)/(e^{b^{∗}t_{n}} −1) where the maximum likelihood estimate b^{∗} is given as the solution to an equation given in Crowder et al. (1991). Conditional p-values for the test based on W can then be computed by using Theorem 1 with ϕ(T) = I(W ≤ w_{obs}). 6.2 Exact Confidence Intervals

Lillegård and Engen (1999) show how the b̂(u, s) of Section 5 can be used to obtain exact confidence intervals for the parameter b, both for the power law and the log linear law cases: Draw a (large) number of independent realizations u_{1}, …, u_{m} of U. Let s be the observed value of the sufficient statistic and let b̃_{(1)} < ⋯ < b̃_{(m)} be the ordered values of the b̂(u_{j}, s). Then (b̃_{(k)}, b̃_{(m−k+1)}) is an exact 1 − 2k/(m + 1) confidence interval for b.

For the power law case it can be seen that the above interval (for m → ∞) is the same as the classical one based on the pivotal statistic 2nb/b^{∗}, which is known to be chi-square distributed with 2(n − 1) degrees of freedom. The interval obtained for the log linear case, however, appears to be new. 6.3 Application to Data Set

We apply the above results to data from a reliability growth program, taken from Leitch (1995, p. 98). There are n = 10 failures, at times 103, 315, 801, 1183, 1345, 2957, 3909, 5702, 7261, 8245. Suppose first that it is of interest to know if the data are consistent with a power law or a log linear NHPP.

In the power law case we get w_{obs} = 0.1263. Simulating the distribution of W we find P_{H_{0}}(W ≤ 0.1263) ≈ 0.024, which implies some evidence against the power law NHPP. In the log linear case we get w_{obs} = 0.1466, s = (8245, 23576) and P_{H_{0}}(W ≤ w_{obs}|S = s) ≈ 0.217, which would not lead to rejection of the log linear NHPP assumption.

An exact 90% confidence interval for b in the power law model, using the chi-square distribtuion as explained above, follows by first computing the maximum likelihood estimate b^{∗} = 0.6249. The resulting interval is (0.2933, 0.9020). Since the interval does not contain 1, we can conclude that we have reliability growth. For the log linear case the maximum likelihood estimate is b^{∗} = −1.75 · 10^{−4} and an exact 90% confidence interval is (−5.60 · 10^{−4}, −3.83 · 10^{−5}), again indicating reliability growth.

References

1. Baker, R. D. (1996). Some new tests of the power law process. Technometrics, 38, 256-265.

2. Crowder, M. J., Kimber, A. C., Smith, R. L. and Sweeting, T. J. (1991). Statistical Analysis of Reliability Data. Chapman & Hall, Great Britain.

3. Engen, S. and Lillegård, M. (1997). Stochastic simulations conditioned on sufficient statistics. Biometrika, 84, 235-240.

4. Gaudoin, O. (1999). CPIT goodness-of-fit tests for reliability growth models. In: Ionescu, V. and Limnios, N. (editors): Statistical and Probabilistic Methods in Reliability, Birkhauser Boston, pp. 27-37.

5. Leitch, R. D. (1995). Reliability Analysis for Engineers. Oxford University Press, New York.

6. Lillegård, M. and Engen, S. (1999). Exact confidence intervals generated by conditional parametric bootstrapping. Journal of Applied Statistics, 26, 447-459.

7. Lindqvist, B. H. and Taraldsen, G. (2000). Monte Carlo conditioning on a sufficient statistic. (In preparation.)