Bayesian Detection of Changes of a Poisson Process Monitored at Discrete Time Points Where the Arrival Rates are Unknown

Abstract

We look at a Poisson process where the arrival rate changes at some unknown integer. At each integer, we count the number of arrivals that happened in that time interval. We assume that the arrival rates before and after the change are unknown. For a loss function consisting of the cost of late detection and a penalty for early stopping, we develop, using dynamic programming, the one- and two-step look-ahead Bayesian stopping rules. We provide some numerical results to illustrate the effectiveness of the detection procedures.

Keywords:

• Bayesian stopping rules
• Change-point detection
• Dynamic programming
• Poisson processes
• Risk
• Unknown arrival rates

Subject Classifications:

• 62L15
• 60G40

1. INTRODUCTION

Detection of changes in the distribution of random variables has become very important in many aspects of life today. When there is an increase in the arrival rates of patients coming to a hospital, it is important to detect this change as soon as possible. This could be due to environmental factors or other issues. This is also important in industry, where quality control depends upon being able to detect changes in the process mean as soon as possible.

There are many published papers dealing with the topic of detection of changes in the distribution of data. Shewhan's (1926) control charts for means and standard deviations, Page's (1954) continuous inspection schemes, and Page's (1962) cumulative sum schemes are all process control procedures designed to detect changes in the distributions of sequences of observed data. Chernoff and Zacks (1964), Hinkley and Hinkley (1970), Kander and Zacks (1966), and Smith (1975) investigated the Bayesian detection procedures for fixed a sample sizes. See Zacks (1983) for a survey of papers on change-point problems. Shiryaev (1978) studied optimal Bayesian procedures for detecting changes when the process is stopped as soon as a change is detected. Zacks and Barzily (1981) extended the work of Shiryaev for the case of Bernoulli sequences. Optimal stopping rules were described based on dynamic programming procedures. In general, finding optimal stopping rules is quite difficult. Zacks and Barzily gave explicit formulas for finding one- and two-step-ahead suboptimal procedures.

Several studies were published recently on detecting changes in the intensity of a homogeneous ordinary Poisson process. Among these studies, we mention Peskir and Shiryaev (2002), Herberts and Jensen (2004), and Brown and Zacks (2006a). These papers dealt with a Poisson process that is monitored continuously.

Brown and Zacks (2006b) also studied a Poisson process that is monitored only at discrete time points. In that paper, we look at the sequence of random variables X _i , where X _i is the number of arrivals that occur in the time interval (i − 1, i]. Thus, we see only the number of arrivals that happen in each time interval, not exactly where the arrivals occurred within that time interval. In all of the above papers, it was assumed that the arrival rates both before and after the change are known. In the present paper, we assume that the arrival rates before and after the change are unknown and the change is positive. We assume that, the change point τ is an unknown integer. We use a Bayesian approach. We put Shiryaev's (1978) geometric prior distribution on the change point τ and a uniform 0 < λ₁ < λ₂ < M prior distribution on the two arrival rates λ₁ and λ₂.

In Section 2, we discuss the posterior process of calculating the probability that a change has already occurred by time n. In Section 3, we discuss the optimal stopping rule based on dynamic programming procedures. We assume there is a cost associated with stopping early and a cost associated with late detection. We develop optimal stopping rules that will minimize the expected cost. In Section 4, we give an explicit formulation of the two-step-ahead stopping rule. We conclude with a numerical example showing how our method works.

2. THE BAYESIAN FRAMEWORK

We have a Poisson process, but we are only able to observe this process at fixed discrete time points. Let X _i be the number of arrivals in the interval (i − 1, i]. We denote ℑ _n as the σ-field generated by {X ₁,…, X _n }. At the end of some unknown interval indexed by an integer τ, the rate of arrivals increases from λ₁ to λ₂. We assume the arrival rates before and after the change, namely, λ₁ and λ₂, are unknown, as is the change point τ. Our goal is to find a detection procedure that will detect the change point τ as soon as possible. The event {τ = 0} represents the case where the change happens prior to the first observed interval. In this case all the observations are taken from a process with the arrival rate at λ₂. The event {τ = r} means that the change happens at the end of the rth interval. Define as the total number of arrivals before time r and as the total number of arrivals after time r. T ₀ ≡ 0 and T ₀∗ = 0. Thus, when τ = r, T _r is a Poisson (rλ₁) and T∗ _n−r is Poisson ((n − r)λ₂). If (λ₁, λ₂) are both known, the likelihood of τ is given by

When (λ₁, λ₂) are unknown, the profile likelihood function of τ is the expected value of (2.1) with respect to the prior distribution of (λ₁, λ₂). We assume a priori that (λ₁, λ₂) is uniformly distributed in some interval 0 < λ₁ < λ₂ < M. That is, we know that when the change happens, the arrival rates have increased, and we also know that the arrival rates are bounded and uniform in a certain interval. Thus the profile likelihood function of τ, given ℱ _n , is

For r = 0, let

where is the survival function of the Poisson distribution with parameter λ. For 1 ≤ r ≤ n − 1, let

Finally, let

We put Shiryaev's geometric prior distribution h on the change point τ,

Thus, we obtain that the posterior probability that there has been a change by time n is

Let  = 1 − ψ _n be the probability that there has not been a change by time n.

3. OPTIMAL STOPPING RULE

In this section, we find an optimal stopping rule. Suppose the cost of stopping early is c ₁ loss units and the cost per time unit of stopping late is c ₂ loss units. We want to find a stopping rule that will minimize the expected cost. We assume without loss of generality that c ₁ = 1 and c ₂ = c is the relative cost of stopping late. Let R _n be the risk (or minimal expected loss) at time n. Therefore,

The stopping rule, then, would be to stop sampling when the risk of stopping is smaller than the risk of continuing. According to (3.1), this rule is equivalent to saying we stop sampling at the first n such that

In order to calculate R _n , we use the principle of dynamic programming. Suppose we must stop after n∗ observations if we haven't stopped before. Let be the risk at time n when at most j more observations are allowed. When j = 0,n∗ = n. Therefore at time n we must stop and our risk is

In general, for j ≥ 1, we have a choice: we can stop and say that a change has occurred, or we can continue and take another observation. If we decide to continue, our risk is . Therefore, our risk at time n is . We can further write this as

Theorem 3.1

The sequence {(,ℑ _n ), n ≥ 1} is a super martingale.

Proof

Taking the conditional expectation of ψ _n+1 given the σ-field generated by the first n observations, we obtain

Therefore,

Accordingly, the sequence {,n ≥ 0} is a super martingale.

Corollary 3.1

Proof

From equation (3.6),

Thus, the lim sup _n→∞ E() = 0. Since  ≥ 0 with probability 1 for all n, we obtain (3.7).

Thus, the one-step look-ahead procedure states that we stop sampling and declare that a change has happened at the first n such that

Define for j ≥ 1, and . Thus, the functions satisfy the recursive relationship

Hence, the risk can be written as a function of , namely,

We then make the decision to stop when , which is equivalent to stopping when . We can define the jth-step-ahead boundary to be .

Lemma 3.1

exists. That is, there exists an optimal stopping boundary.

Proof

By induction one can show that are decreasing in j. Thus, are decreasing functions of j. Also, for all j. By the monotone convergence theorem, the limit exists.

Lemma 3.2

.

Proof

Let

Therefore, if i ∊ I, then  < q∗. Since by Corollary 2.1,  → 0,

Therefore, . Thus, . By induction, we can show this for all j.

4. EXPLICIT FORMULATION OF THE TWO-STEP-AHEAD BOUNDARY

In this section, we evaluate the two-step-ahead boundary. Substituting j = 2 into equation (3.10), we obtain

To evaluate the two-step-ahead risk, we find an expression for . First, we express recursively as a function of :

where, according to (2.5),

and is defined as L _r when the likelihood is being calculated with respect to ℑ _n , the σ-field generated by the first n observations. According to (3.10),

We only need to take the expectation over those X _n+1 for which  > q∗.

Lemma 4.1

If X _n+1 ≥ M, then  < .

Proof

Let X _n+1 = m ≥ M. Let . For r = 0,

If we divide equation (4.5) by given in (4.3), we obtain, , where F and G are defined as

We show that . This is equivalent to showing that the first term in F is greater than the first term in G. (See Lemma A.1 Appendix for details). Therefore, . Now, let 0 < r < n:

To show that , it suffices to show that each term in the numerator times is greater than the corresponding term in the denominator of (4.7) for each j, r, and n. This is equivalent to showing F _j  > G _j , where F _j and G _j axe given by

The inequality F _j  > G _j can be shown by induction on j, making use of Lemma A.1 in the appendix. Therefore, for all r < n. Thus,

From (4.8) and (4.9), we see that . Therefore,

Therefore, if the (n + 1)th observation is greater than M, the posterior probability that a change has not occurred will decrease.

We use Lemma 4.1, to evaluate the two-step-ahead stopping rule. According to equation (4.1), the two-step-ahead risk function is where is defined as . We only need to take the expectation over those X _n+1 for which  > q∗. Note that if  < q∗, then  < q∗ for k ≥ M. Thus for  < q∗,

Thus the two-step-ahead boundary is

Notice when  < q∗. Also when  > q∗, both the one- and two-step look ahead procedure will tell us to continue sampling. Thus the two-step-ahead procedure is equivalent to the procedure which tells us to stop when  < B _{n, 2}. In Table 1, we do an example. We look at simulated values for λ₁ = 1, when n ≤ 10 and λ₂ = 2 when n > 10. We use a prior distribution of p = 0.01, π = 0.01. We assume we know that the upper bound on the arrival rate M = 3, and the cost c of stopping late is c = 0.06. In this example, using formula 3.9, we calculate the one-step-ahead boundary q∗ = 0.8571. We show at each n, the total number of arrivals up to time n, the posterior probability that a change has not taken place and the two-step-ahead boundary B _{n, 2}. One can observe that the B _{n, 2} approaches q∗ for large n. This illustrates Lemma 3.2.

Table 1. Simulated values to show detection using the 2-step-ahead boundary

n	T n		B n, 2
1	0	0.9770	0.7788
2	1	0.9828	0.7907
3	1	0.9877	0.7786
4	4	0.9055	0.8236
5	4	0.9520	0.8100
6	6	0.9186	0.8246
7	6	0.9484	0.8144
8	7	0.9450	0.8164
9	8	0.9420	0.8181
10	9	0.9392	0.8196
11	14	0.8018	0.8552
12	15	0.8231	0.8526
13	17	0.7960	0.8571
14	19	0.7704	0.8571
15	23	0.6370	0.8571
16	24	0.6744	0.8571
17	28	0.5437	0.8571
18	29	0.5859	0.8571
19	30	0.6217	0.8571
20	31	0.6518	0.8571

From looking at Table 1 one observes that we would stop at time 11. Since the change took place at time 10, this shows that the detection rule is quick. We also see the convergence between the one- and two-step-ahead stopping rules.

ACKNOWLEDGMENTS

I would like to thank Dr. Shelemyahu Zacks for his helpful advice in preparing this article. I would also like to thank the Editor, Associate Editor, and referees for their helpful suggestions.

Notes

Recommended by N. Mukhopadhyay

Appendix A

Lemma A.1

For all n, k, and m positive integers, and m ≥ M > 0,

Proof

Taking differences and utilizing the relationship between the cumulative distribution function of the Poisson distribution and the gamma distribution, we obtain

Simplifying, we obtain that equation (rm A.2) is equal to

We will divide the integral into where y < x and a second integral where y > x. Thus, equation (rm A.3) is equivalent to M ^2k+m+1 times

Note that the second integral in (rm A.4) can be rewritten as

Thus, equation (rm A.2) is

g (x) = x ^m e ^−Mx is an increasing function over (0, 1) when m > M. Hence,

for all 0 < y < x < 1. Therefore,

for all n, k, m, and M such that m ≥ M.