PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Comput Stat Data Anal. Author manuscript; available in PMC 2010 July 1.
Published in final edited form as:
Comput Stat Data Anal. 2009 July 1; 53(9): 3314–3323.
doi:  10.1016/j.csda.2009.02.006
PMCID: PMC2678871
NIHMSID: NIHMS96244

Non-iterative sampling-based Bayesian methods for identifying changepoints in the sequence of cases of haemolytic uraemic syndrome

Summary

Diarrhoea-associated Haemolytic uraemic syndrome (HUS) is a disease that affects the kidneys and other organs. Motivated by the annual number of cases of HUS collected in Birmingham and Newcastle of England, respectively, from 1970 to 1989, we consider Bayesian changepoint analysis with specific attention to Poisson changepoint models. For changepoint models with unknown number of changepoints, we propose a new non-iterative Bayesian sampling approach (called exact IBF sampling), which completely avoids the problem of convergence and slow convergence associated with iterative Markov chain Monte Carlo (MCMC) methods. The idea is to first utilize the sampling inverse Bayes formula (IBF) to derive the conditional distribution of the latent data given the observed data, and then to draw iid samples from the complete-data posterior distribution. For the purpose of selecting the appropriate model (or determining the number of changepoints), we develop two alternative formulae to exactly calculate marginal likelihood (or Bayes factor) by using the exact IBF output and the point-wise IBF, respectively. The HUS data are re-analyzed using the proposed methods. Simulations are implemented to validate the performance of the proposed methods.

Keywords: Bayes factor, Changepoint problem, Haemolytic uraemic syndrome, IBF sampling, MCMC, Non-iterative Bayesian approach, Poisson distribution

1. Introduction

Diarrhoea-associated Haemolytic uraemic syndrome (HUS) is a disease that affects the kidneys and other organs. It poses a substantial threat to infants and young children as one of the leading causes of both acute and chronic kidney failures. HUS is most common in the warmer months of the year, following a gastrointestinal illness caused primarily by a particular strain of bacterium, Escherichia Coli O157:H7 (Milford et al., 1990). These bacteria (E. Coli O157:H7) produce extremely potent toxins which are the main cause of the symptoms related to the gastrointestinal illness. Table 1 displays the annual number of cases of HUS collected in Birmingham and Newcastle of England, respectively, from 1970 to 1989 (Tarr et al., 1989; Henderson and Matthews, 1993). The primary concern is the incidence of HUS and when the frequency of cases increases sharply. In the mean-corrected cumulative sum plot (Figure 1), the annual totals appear to increase abruptly at about 1980 for the Birmingham series and 1976, 1984 for the Newcastle series. Therefore, a changepoint analysis of the data with Poisson models seems to be appropriate.

Figure 1
Mean-corrected cumulative sum plot for the number of cases at Birmingham and Newcastle.
Table 1
Counts of cases of HUS at Birmingham and Newcastle (Tarr et al., 1989)

Changepoint problems (CPPs) are often encountered in medicine and other fields, e.g., economics, finance, psychology, signal processing, industrial system control and geology. Typically, a sequence of data is collected over a period of time, we wish to make inference about the location of one or more points of the sequence at which there is a change in the model. The literature on CPPs is extensive. For Poisson process CPPs, a well-known example concerns British coal-mining disasters from 1851–1962 (originally gathered by Maguire et al. (1952) and corrected by Jarrett (1979)). Frequentist investigations appear in Worsley (1986) and Siegmund (1988), while traditional Bayesian analysis and Markov chain Monte Carlo (MCMC) hierarchical Bayesian analysis are presented in Raftery and Akman (1986) and Carlin, Gelfand and Smith (1992), respectively. Arnold (1993) considered the application of the Gibbs sampler to a Poisson distribution with a changepoint. For binomial CPPs, Smith (1975) presented the conventional Bayesian approach for a finite sequence of independent observations with details on binomial single-changepoint model. Smith (1980) studied binomial multiple-changepoint model, which were investigated by Stephens (1994) using the Gibbs sampler. For binary CPPs, Halpern (1999) applied a novel changepoint statistic based on the minimum value, over possible changepoint locations of Fisher’s Exact Test to assessing recombination in genetic sequences of HIV. For multiple change-point models, Chib (1998) provided a comparison study and Fearnhead and Liu (2007) proposed an on-line algorithm. Three comprehensive reviews on CPPs are provided by Brodsky and Darkhovsky (1993), Chen and Gupta (2000) and more recently by Wu (2005).

The primary objective in the analysis of CPPs is to make inferences about unknown changepoints and the associated parameters. Although the MCMC methods can be employed in such Bayesian analyses, in our viewpoint, the difficulties lie in monitoring the convergence of the Markov chains. In addition, they could suffer from slow convergence. These issues have prompted some researchers to take the view that the MCMC methods are to be used only when there is no better alternative (see, e.g., discussions in Evans and Swartz (1995, 2000) and Hobert and Casella (1996)). In this article, we first propose a new non-iterative Bayesian sampling approach (called exact IBF sampling), which completely avoids the problem of convergence and slow convergence. The idea is to first utilize the sampling-wise inverse Bayes formulae (IBF, Tan et al., 2003) to derive the conditional distribution of the missing data given the observed data, and then to draw iid samples from the complete-data posterior distribution.

In practice, we are generally uncertain about the number of changepoints. Hence, model determination is the first task in changepoint analysis. Let Ms represent a model with s changepoints. A classical approach of selecting the most appropriate model is the likelihood ratio test by comparing Ms with Ms+1 (e.g., Henderson and Matthews, 1993). Gelfand and Dey (1994) reviewed the behavior of the likelihood ratio statistic and well-known adjustments to it. In the context of Bayesian analysis, Bayes factor is a useful tool for model choice. However, the calculation of Bayes factor itself has proved extremely challenging (Kass and Raftery, 1995). Approximate computation of Bayes factor (equivalently, marginal likelihood) can be implemented by using the Gibbs output (Chib, 1995) or the more general MCMC output (Chen, 2005). In this paper, we will we develop two alternative formulae to exactly calculate marginal likelihood (or Bayes factor) by using the exact IBF output and the point-wise IBF (Tan et al., 2003), respectively.

The rest of this paper is organized as follows. In Section 2, we formulate the general Bayesian changepoint models. In Section 3, we propose a non-iterative Bayesian sampling approach (called the exact IBF sampling) and derive two simple formulae to exactly calculate the marginal likelihood. Section 4 considers Poisson models with single and multiple changepoints and the corresponding Bayesian model selection. In Section 5, we re-analyze the HUS data using the proposed methods. Simulations are conducted to validate the performance of the proposed methods in Section 6. Conclusion and comment are presented in Section 7.

2. Bayesian formulation for changepoint problems

Let Yobs={yi}i=1n denote a realization of the sequence of independent random variables {Yi}i=1n of length n. The random variables {Yi}i=1n are said to have a changepoint at r (1 ≤ rn) if Yi ~ f1(y|θ1) (i = 1,…, r) and Yi ~ f2(y|θ2) (i = r + 1,…, n), where f1(y|θ1) ≠ f2(y|θ2), θ1 and θ2 could be vector-valued. In particular, the point r = n represents ‘no change’. Thus, the likelihood function becomes

L(Yobsr,θ1,θ2)=Πi=1rf1(yiθ1)·Πi=r+1nf2(yiθ2).
(2.1)

Using π(r, θ1, θ2) as a joint prior distribution for r, θ1, and θ2, the joint posterior distribution is given by

f(r,θ1,θ2Yobs)L(Yobsr,θ1,θ2)·π(r,θ1,θ2).
(2.2)

This single-changepoint problem can be easily generalized to incorporate multiple changes in the sequence. The Bayesian formulation for the multiple-changepoint problem is almost identical with that for the single-changepoint problem. Let Ms represent a model with s changepoints denoted by r = (r1,…, rs)[top top]. Similar to (2.2), under Ms (s is given), we have

f(r,θ1,,θs+1Yobs)L(Yobsr,θ1,,θs+1)·π(r,θ1,,θs+1)={Πj=1s+1Πi=rj1+1rjfj(yiθj)}·π(r,θ1,,θs+1),
(2.3)

where r0 [equivalent] 0, rs+1 [equivalent] n, and the changepoints r take values in the domain

S(rYobs)={r:1r1<<rsn,rjisaninteger,j=1,,s}.
(2.4)

The primary objective is to make inferences about the unknown changepoints r and the unknown parameters (θ1,…, θs+1).

3. Exact IBF sampling and marginal likelihood calculation

For a given model, let Yobs denote the observed data, Z the missing data or latent data (e.g., changepoints) and θ the model-specific parameter vector. We further denote the likelihood function by L(Yobs|θ) and the marginal density of Yobs (equivalently, the marginal likelihood) by m(Yobs). Within a Bayesian framework, we assume that the prior density is π(θ). Two basic tasks are (i) for the purpose of Bayesian inferences, how to obtain iid samples from the observed posterior f(θ|Yobs) or equivalently from the joint posterior f(Z, θ|Yobs), and (ii) for the purpose of Bayesian model choice, how to exactly calculate the marginal likelihood m(Yobs) for the given model.

3.1 Exact IBF sampling

In general, we can obtain explicit expressions for both the complete-data posterior distribution f(θ|Yobs, Z) and the conditional predictive distribution f(Z|Yobs, θ), that is, the sampling from the two distributions and the evaluation of the two densities can routinely be implemented. The fundamental conditional sampling principle implies

f(Z,θYobs)=f(ZYobs)·f(θYobs,Z),

which states that if we could draw Z()iidf(ZYobs) and simulate θ([ell]) ~ f(θ|Yobs, Z([ell])), then {(Z(),θ())}1L are iid samples from the joint posterior f(Z, θ|Yobs). Therefore, the key is to be able to generate iid samples from f(Z|Yobs).

Let An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(θ|Yobs) and An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(Z|Yobs) denote the conditional supports of θ|Yobs and Z|Yobs, respectively. The sampling IBF shows that (Tan et al., 2003)

f(ZYobs)f(ZYobs,θ0)f(θ0Yobs,Z),foranarbitraryθ0S(θYobs)andallZS(ZYobs).
(3.1)

Consider the case where Z is a discrete random variable/vector taking finite values on the conditional support An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(Z|Yobs). For example, in (2.2), the changepoint r takes values in {1,…, n}; and in (2.3), the s changepoints (r1,…, rs) take values in An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(r|Yobs) defined by (2.4). Without loss of generality, we denote the conditional support of Z|(Yobs, θ) by An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg (Z|Yobs, θ) = {z1,…, zK}. Since f(Z|Yobs, θ) is available, firstly, we can directly identify {zk}1K from the model specification and thus all {zk}1K are known. Secondly, we assume that {zk}1K do not depend on the value of θ, therefore, we have

S(ZYobs)=S(ZYobs,θ)={z1,,zK}.

Because of the discreteness of Z, the notation f(zk|Yobs) will used to denote the pmf, i.e., f(zk|Yobs) = Pr{Z = zk|Yobs}. Thus, the key is to find pk = f(zk|Yobs) for k = 1, …, K. For some θ0 [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(θ|Yobs), let

qk=qk(θ0)=Pr{Z=zkYobs,θ0}/f(θ0Yobs,zk),k=1,,K.
(3.2)

As both f(Z|Yobs, θ) and f(θ|Yobs, Z) are available, the computation of (3.2) is straight-forward. Observing that all {qk}1K depend on θ0, we denote qk by qk(θ0) to emphasize its dependency on θ0. From the sampling IBF (3.1), we obtain

pk=qk(θ0)/k=1Kqk(θ0),k=1,,K,
(3.3)

where {pk}1K do not depend on θ0 since {pk}1K are normalizing probabilities of {qk}1K. Thus, it is easy to sample from f(Z|Yobs) since it is a discrete distribution with probabilities {pk}1K on {zk}1K (e.g., the built-in S-plus function “ sample” is especially designed for this purpose). We summarize the algorithm as follows.

The exact ibf sampling

Given both the complete-data posterior distribution f(θ|Yobs, Z) and the conditional predictive distribution f(Z|Yobs, θ),

  1. Identify An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(Z|Yobs) = {z1,…, zK} from f(Z|Yobs, θ) and calculate {pk}1K according to (3.3) and (3.2);
  2. Generate iid samples {Z()}=1L of Z from the pmf f(Z|Yobs) with probabilities {pk}1K on {zk}1K;
  3. Generate θ([ell]) ~ f(θ|Yobs, Z([ell])) for [ell]= 1,…,L, then {θ()}1L are iid samples from the observed posterior distribution f(θ|Yobs).

3.2 Exact calculation of marginal likelihood

In this subsection, we provide two alternative formulae to calculate the marginal likelihood m(Yobs). Let {(Z(),θ())}1L denote the output from the exact IBF sampling. From Bayes formula: m(Yobs) = L(Yobs|θ)π(θ)/f(θ|Yobs), which holds for any θ, we have

logm(Yobs)=logL(Yobsθ0)+logπ(θ0)logf(θ0Yobs),θ0S(θYobs).
(3.4)

For estimation efficiency, θ0 is generally taken to be a high-density point in the support of the posterior (e.g, the posterior mode/mean as suggested by Chib (1995)). Since the observed posterior density can be written as

f(θYobs)=f(θYobs,Z)f(ZYobs)dZ,

we obtain a Monte Carlo estimate of f(θ|Yobs) at θ0,

f^(θ0Yobs)=(1/L)=1Lf(θ0Yobs,Z()),
(3.5)

where {Z([ell])} are iid samples from f(Z|Yobs). Note that this estimate is simulation consistent, i.e., f(θ0|Yobs) → f(θ0|Yobs) as L → ∞. Combining (3.4) with (3.5), we have an approximate formula to calculate m(Yobs).

On the other hand, note that Z is a discrete random variable taking values on {zk}1K, using the point-wise IBF (Tan et al., 2003): f(θ|Yobs) = {∫ f(Z|Yobs, θ)/f(θ|Yobs, Z) dZ}−1, we explicitly have

f(θ0Yobs)={k=1Kqk(θ0)}1=p1/q1(θ0),
(3.6)

where pk and qk(θ0) are defined in (3.3) and (3.2), respectively. Substituting (3.6) into (3.4) gives an exact formula to calculate m(Yobs).

4. Poisson models with changepoints

4.1 A single changepoint

Let Ms represent a model with s changepoints, Poisson(θ) a Poisson distribution with mean θ and Poisson(·|θ) the corresponding probability mass function. We first consider the single-changepoint model M1. In (2.1), we let fj(y|θj) = Poisson(y|θj) for j = 1, 2. As a joint prior distribution for (r, θ1, θ2), we assume that r, θ1 and θ2 are independent, r has a discrete uniform distribution on {1,…,n},

θ1Ga(a1,b1)andθ2Ga(a2,b2),
(4.1)

where Ga(a, b) is a gamma distribution with density Ga(x|a, b) = baxa−1ebx/Γ(a), x ≥ 0. Thus, the joint posterior distribution (2.2) becomes

f(r,θ1,θ2Yobs)θ1a1+Sr1e(b1+r)θ1·θ2a2+SnSr1e(b2+nr)θ2,

where Sri=1ryi. Direct calculation yields

f(θ1,θ2Yobs,r)=Ga(θ1a1+Sr,b1+r)·Ga(θ2a2+SnSr,b2+nr),
(4.2)
f(rYobs,θ1,θ2)=(θ1/θ2)Srexp{(θ2θ1)r}i=1n(θ1/θ2)Siexp{(θ2θ1)i},r=1,,n.
(4.3)

We can treat the changepoint r as latent variable Z and (θ1, θ2) as parameter vector θ. By using (3.1)–(3.3), for any given (θ1,0, θ2,0) [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(θ1, θ2|Yobs), we immediately obtain

f(rYobs)=Γ(a1+Sr)Γ(a2+SnSr)/[(b1+r)a1+Sr(b2+nr)a2+SnSr]i=1nΓ(a1+Si)Γ(a2+SnSi)/[(b1+i)a1+Si(b2+ni)a2+SnSi],
(4.4)

where r = 1,…,n. It again confirms that the right-hand side of (4.4) does not depend on (θ1,0, θ2,0). Based on (4.4) and (4.2), we can obtain iid posterior samples for the changepoint r and the parameters (θ1, θ2) by using the exact IBF sampling.

4.2 Multiple changepoints

Now we consider the multiple-changepoints model Ms. In (2.3), let fj(y|θj) = Poisson(θj) for j = 1,…, s + 1, where θ = (θ1,…, θs+1)[top top] is the mean vector and r = (r1,…, rs)[top top] denote the s changepoints taking integer values on the domain An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(r|Yobs) defined in (2.4). We use independent priors for r, θ, and r has a discrete uniform prior on An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(r|Yobs),

θjGa(aj,bj),j=1,,s+1.
(4.5)

Thus, the joint posterior distribution (2.3) becomes

f(r,θYobs)Πj=1s+1θjaj+SrjSrj11e(bj+rjrj1)θj,
(4.6)

where Sri=1ryi, r0 [equivalent] 0 and rs+1 [equivalent] n. From (4.6), we have

f(θYobs,r)=Πj=1s+1Ga(θjaj+SrjSrj1,bj+rjrj1),
(4.7)
f(rYobs,θ)Πj=1s(θj/θj+1)Srje(θj+1θj)rj,rS(rYobs).
(4.8)

We treat r as latent variables and θ as parameter vector. By using (3.1)–(3.3), for any given θ0 [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms96244ig1.jpg(θ|Yobs), we immediately obtain

f(rYobs)j=1s+1Γ(aj+SrjSrj1)(bj+rjrj1)aj+SrjSrj1,rS(rYobs).
(4.9)

Based on (4.9) and (4.7), we can obtain iid posterior samples for the changepoints r and the parameter vector θ by using the exact IBF sampling.

4.3 Determining the number of changepoints via Bayes factor

In practice, we are generally uncertain about the number of changepoints. Hence, model determination is the first task in changepoint analysis. Let Ms represent the Poisson model with s changepoints r = (r1,…, rs)[top top] and θ = (θ1,…, θs+1)[top top] the mean vector. Further let Θ = (r, θ) and [Theta with circumflex] = ([r with circumflex], [theta w/ hat]) denote the posterior means obtained via the exact IBF output. Under model Ms, from (3.4), the log-marginal likelihood is given by

logm(YobsMs)=logL(YobsΘ^,Ms)+logπ(Θ^Ms)logf(Θ^Yobs,Ms),
(4.10)

where f([Theta with circumflex]|Yobs,Ms) = f([r with circumflex]|Yobs,Ms) · f([theta w/ hat]|Yobs, [r with circumflex],Ms). We choose the model with the largest log-marginal likelihood. Essentially, the marginal likelihood approach is the same as the Bayes factor approach. A Bayes factor is defined as the ratio of posterior odds versus prior odds, which is simply a ratio of two marginal likelihoods. For comparing models Ms and Ms+1, the Bayes factor for model Ms vs. model Ms+1 is

Bs,s+1=m(YobsMs)m(YobsMs+1).
(4.11)

Jeffreys (1961, Appendix B) suggested interpreting Bs,s+1 in half-units on the log10 scale, i.e., when Bs,s+1 falls in (1, 3.2), (3.2, 10), (10, 100) and (100, +∞), the evidence against Ms+1 is not worth more than a bare mention, substantial, strong and decisive, respectively.

5. Analysis of the HUS data

We first analyze the Birmingham data in Table 1. Denote the number of cases of HUS in Birmingham in year i by yi (i = 1,…,n with n = 20, and i = 1 denotes the year 1970). To determine the number of changepoints via Bayes factor, we can not use non-informative prior distributions because they are improper. We investigate models M0, M1 and M2 and choose standard exponential prior distributions, specified by setting all aj = bj = 1 in (4.5). Based on (4.10), we calculate log-marginal likelihoods for the three models, and we obtain logm(Yobs|M0) = −86.14, logm(Yobs|M1) = −57.56 and logm(Yobs|M2) = −57.00. Therefore, M2 seems to be the most appropriate choice. From (4.11), the Bayes factor for M1 versus M0 is 2.583 × 1012, while the Bayes factor for M2 versus M1 is 1.751. That is, the difference between M2 and M1 is not worth to mention. Therefore, we select M1, which is consistent with the pattern indicated in Figure 1.

Under M1, we assume that y1,,yriidPoisson(θ1) and yr+1,,yniidPoisson(θ2), where r is the unknown changepoint and θ1θ2. Table 2 contains the exact posterior probabilities for the changepoint r using (4.4). The changepoint occurs at r = 11 (i.e., year 1980) with posterior probability 0.9795. Based on (4.4) and (4.2), we generate 20, 000 iid posterior samples by using the exact IBF sampling, and the Bayes estimates of r, θ1 and θ2 are given by 11.013, 1.593 and 9.609. The corresponding Bayes standard errors are 0.143, 0.370 and 0.985. The 95% Bayes credible intervals for θ1 and θ2 are [0.952, 2.393] and [7.800, 11.621], respectively. Figures 2(a) and 2(b) show the histogram of the changepoint r and the posterior densities of θ1 and θ2, which are estimated by a kernel density smoother based on iid posterior samples. Figure 2(c) depicts the annual numbers of HUS for Birmingham series, the identified changepoint position, and the average number of cases before and after the changepoint.

Figure 2
Birmingham data set. (a) Histogram of the changepoint r. (b) The posterior densities of θ1 and θ2 estimated by a kernel density smoother based on 20, 000 iid samples generated via the exact IBF sampling. (c) The annual numbers of cases ...
Table 2
Exact posterior probabilities for the changepoint r for Birmingham series

Now we analyze the Newcastle data in Table 1. Similarly, three log-marginal likelihoods are given by logm(Yobs|M0) = −85.24, logm(Yobs|M1) = −64.13 and logm(Yobs|M2) = −64.10. From (4.11), the Bayes factor for M2 versus M0 is 1.5169 × 109, and the Bayes factor for M2 versus M1 is 1.03. Therefore, we select M2, which is consistent with the pattern as indicated in Figure 1. In addition, the selection of M2 is also identical to that obtained by Henderson and Matthews (1993).

Under M2, we assume that y1,,yr1iidPoisson(θ1),yr1+1,,yr2iidPoisson(θ2), and yr2+1,,yniidPoisson(θ3), where (r1, r2) are the unknown changepoints and θ1θ2θ3. Using the standard exponential prior distributions, specified by letting aj = bj = 1 (j = 1, 2, 3) in (4.5), we obtained exact joint posterior probabilities for the changepoint pair (r1, r2) from (4.9). Two changepoints occur at r1 = 7 and r2 = 15 (i.e., year 1976 and year 1984) with the joint posterior probability being 0.3589. Based on (4.9) and (4.7), we generated 20, 000 iid posterior samples. The resulting Bayes estimates of r1, r2, θ1, θ2 and θ3 are given by 7.638, 15.47, 1.805, 3.591 and 9.643. The 95% Bayes credible intervals for θ1, θ2 and θ3 are [0.7461, 3.620], [1.5085, 11.50] and [0.2806, 13.32], respectively. Figures 3(a) and 3(b) display the histograms of r1 and r2. Figures 3(c) shows the posterior densities of θj (j = 1, 2, 3). Figure 3(d) depicts the annual numbers of HUS, two identified changepoints, and the average number of cases before and after the two changepoints.

Figure 3
Newcastle data set. (a) Histogram of the changepoint r1. (b) Histogram of the changepoint r2. (c) The posterior densities of θ1, θ2 and θ3 estimated by a kernel density smoother based on 20, 000 iid samples generated via the exact ...

6. Simulation Studies

The first simulated dataset consists of 100 observations with y1,,y50iidPoisson(3) and y51,,y100iidPoisson(0.5). The simulated observations are shown in Figure 4(c). We again use standard exponential distributions as priors of θ. From (4.10), log-marginal likelihoods for three models M0, M1 and M2 are given by logm(Yobs|M0) = −187.1, logm(Yobs|M1) = −148.8 and logm(Yobs|M2) = −149.3. From (4.11), we have B10 = 4.3 × 1016 and B12 = 1.649, which suggest that M1 is appropriate. Computations according to (4.4) show that the changepoint occurs at r = 50 with posterior probability 0.807. Based on (4.4) and (4.2), we generate 30, 000 iid posterior samples by using the exact IBF sampling. The Bayes means, standard errors, and 95% credible intervals for (r, θ1, θ2) are given by (50.6, 2.8226, 0.5249), (1.642, 0.240, 0.103) and [50, 56], [2.373, 3.315], [0.341, 0.742], respectively. Figures 4(a) and 4(b) show the histogram of r and the posterior densities of θ1 and θ2. Figure 4(c) displays the 100 simulated observations, the identified changepoint position, and the Bayes estimates of θ1 and θ2.

Figure 4
Simulated dataset with one changepoint. (a) Histogram of r. (b) The posterior densities of θ1 and θ2. (c) The 100 simulated observations. The dotted vertical line denotes the identified changepoint position (r = 50), the left horizontal ...

The second simulated dataset consists of 100 observations: y1,,y20iidPoisson(5.5),y21,,y70iidPoisson(0.8), and y71,,y100iidPoisson(3.5). The simulated observations are shown in Figure 5(d). Similarly, we have logm(Yobs|M0) = −249.7, logm(Yobs|M1) = −222.2 and logm(Yobs|M2) = −185.6, B20 = 6.891 × 1027, and B21 = 7.856 ×, which suggest that M2 is appropriate. Computations according to (4.9) show that the changepoints occur at r1 = 20 and r2 = 70 with the joint posterior probability 0.7912. Based on (4.9) and (4.7), we generate 30, 000 iid posterior samples by using the exact IBF sampling. The Bayes means, standard errors, and 95% credible intervals for (r1, r2, θ1, θ2, θ3) are given by (20.0091, 69.7871, 5.7120, 0.8427, 3.8190), (0.1019, 0.4457, 0.5230, 0.1296, 0.3517) and [20, 20], [69, 70], [4.735, 6.789], [0.610, 1.116], [3.161, 4.537], respectively. Figures 5(a) and 4(b) show the histogram of r1 and r2. Figures 5(c) shows the posterior densities of θ1, θ2 and θ3. Figure 5(d) displays the 100 simulated observations, the two identified changepoint positions, and the Bayes estimates of θ1, θ2 and θ3.

Figure 5
Simulated dataset with two changepoints. (a) Histogram of r1. (b) Histogram of r2. (c) The posterior densities of θ1, θ2 and θ3. (d) The 100 simulated observations, two identified changepoints (20, 70), and three Bayes estimates ...

7. Discussion

It is noted that Barry and Hartigan (1992, 1993) and Fearnhead (2006) describe methods for calculating marginal likelihoods for multiple change-point problems. The latter also discusses methods for simulating from the change-point positions. Barry and Hartigan (1992, 1993) assumed a specific prior structure on the number and position of changepoints; while Fearnhead (2006) considered models with a fixed number of change-points. Although it was claimed that these methods can deal with arbitrarily large numbers of change-points, these methods are quite complicated in implementation. For example, the parameter values may be estimated exactly in O(n3) calculations, or to an adequate approximation by MCMC methods that are O(n) in the number of observations.

In this paper, we considered Poisson changepoint analysis by using an exact IBF sampling approach. The advantages of the proposed exact IBF sampling method over MCMC methods are that: (i) there is no requirement to diagnose whether the MCMC algorithms has converged, i.e., the former entirely avoids the problems of convergence and slow convergence associated with MCMC methods; (ii) because the samples generated from the observed posterior distribution are independent it is straightforward to quantify uncertainty in estimates of features of the posterior distributions based on them. To determine the number of changepoints, we developed simple methods to exactly calculate marginal likelihood (or Bayes factor). Two simulations are conducted to validate the performance of the proposed methods.

We should point out that the proposed approach is limited. For example, let Ms represent a model with s changepoints and the s-changepoint is denoted by r = (r1,…, rs)[top top]. For a large number of observations or s ≥ 4, the calculation of (4.9) becomes prohibitive. In this cases, the general IBF sampler (Tian et al., 2003) is a feasible alternative.

In the re-analysis of the HUS data using the proposed methods, we have focused on the annual numbers of cases rather than the incidence of the HUS because accurate population are difficult to obtain for the catchment areas. Possibilities for further analysis might consider extra-Poisson variation, treads in mean, the influence of some covariates (e.g., age, race) and so on.

Acknowledgments

We are grateful to the Editor, an Associate Editor and three referees for their constructive comments and suggestions. GL Tian and M Tan’s research was supported in part by U.S. National Cancer Institute grants CA106767 and CA119758. The research of KW Ng was partially supported by a research grant of the University of Hong Kong. Special thanks should go to one referee for drawing our attention to several recent papers on multiple change-point problems.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errorsmaybe discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  • Arnold SF. Gibbs sampler. In: Rao CR, editor. Handbook of Statistics. Vol. 9. Elsevier Science Publishers B. V; 1993. pp. 599–625.
  • Barry D, Hartigan JA. Product partition models for change point problems. Ann Statist. 1992;20:260–279.
  • Barry D, Hartigan JA. A Bayesian analysis for change point problems. Journal of the American Statistical Association. 1993;88:309–319.
  • Brodsky BE, Darkhovsky BS. Series: Mathematics and Its Applications. Vol. 243. Springer; New York: 1993. Nonparametric Methods in Change Point Problems.
  • Carlin BP, Gelfand AE, Smith AFM. Hierarchical Bayesian analysis of changepoint problems. Appl Statist. 1992;41:389–405.
  • Chen MH. Computing marginal likelihoods from a single MCMC output. Statistica Neerlandica. 2005;59:16–29.
  • Chen J, Gupta AK. Parametric Statistical Change Point Analysis. Springer; New York: 2000.
  • Chib S. Marginal likelihood from the Gibbs output. Journal of the American Statistical Association. 1995;90:1313–1321.
  • Chib S. Estimation and comparison of multiple change-point models. Journal of Econometrics. 1998;86:221–241.
  • Evans M, Swartz T. Methods for approximating integrals in statistics with special emphasis on Bayesian integration problems (with discussions) Statist Sci. 1995;10:254–272.
  • Evans M, Swartz T. Approximating Integrals via Monte Carlo and Deterministic Methods. Oxford University Press; Oxford: 2000.
  • Fearnhead P. Exact and efficient inference for multiple changepoint problems. Statist Comput. 2006;16:203–213.
  • Fearnhead P, Liu Z. On-line inference for multiple changepoint problems. Journal of the Royal Statistical Society, Series B. 2007;64:589–605.
  • Gelfand AE, Dey DK. Bayesian model choice: asymptotics and exact calculations. Journal of the Royal Statistical Society, Series B. 1994;56:501–514.
  • Halpern AL. Minimally selected p and other tests for a single abrupt changepoint in a binary sequence. Biometrics. 1999;55:1044–1050. [PubMed]
  • Henderson R, Matthews JNS. An investigation of changepoints in the annual number of cases of haemolytic uraemic syndrome. Appl Statist. 1993;42:461– 471.
  • Hobert JP, Casella G. The effect of improper priors on Gibbs sampling in hierarchical linear models. J Am Statist Assoc. 1996;91:1461–1473.
  • Jarrett RG. A note on the intervals between coal-mining disasters. Biometrika. 1979;66:191–193.
  • Jeffreys H. Theory of Probability. 3. Oxford: Oxford University Press; 1961.
  • Kass RE, Raftery AE. Bayes factors. Journal of the American Statistical Association. 1995;90:773–795.
  • Maguire BA, Pearson ES, Wynn AHA. The time intervals between industrial accidents. Biometrika. 1952;38:168–180.
  • Milford DV, Taylor CM, Gutteridge B, Hall SM, Rowe B, Kleanthous H. Haemolytic uraemic syndromes in the British Isles 1985–8: association with verocytotoxin producing Escherichia coli. Part 1: Clinical and epidemiological aspects. Arch Dis Child. 1990;65(7):716–721. [PMC free article] [PubMed]
  • Raftery AE, Akman VE. Bayesian analysis of a Poisson process with a change-point. Biometrika. 1986;73:85–89.
  • Siegmund D. Confidence sets in change point problems. Int Statist Rev. 1988;56:31–48.
  • Smith AFM. A Bayesian approach to inference about a change-point in a sequence of random variables. Biometrika. 1975;62:407–416.
  • Smith AFM. Change-point problems: Approaches and applications. In: Bernardo JM, DeGroot MH, Lindley DV, Smith AFM, editors. Bayesian Statistics. Vol. 1. Valencia: Valencia University Press; 1980. pp. 83–89.
  • Smith AFM, Cook DG. Straight lines with a change-point: A Bayesian analysis of some renal transplant data. Appl Statist. 1980;29:180–189.
  • Stephens DA. Bayesian retrospective multiple-changepoint identification. Appl Statist. 1994;43:159–178.
  • Tan M, Tian GL, Ng KW. A noniterative sampling method for computing posteriors in the structure of EM-type algorithms. Statistica Sinica. 2003;13:625–639.
  • Tarr PI, Neill MA, Allen J, Siccardi CJ, Watkins SL, Hickman RO. The increasing incidence of the hemolytic-uremic syndrome in King County, Washington: lack of evidence for ascertainment bias. Am J Epidemiol. 1989;129(3):582–586. [PubMed]
  • Worsley KJ. Confidence regions and tests for a change-point in a sequence of exponential family random variables. Biometrika. 1986;73:91–104.
  • Wu YH. Series: Lecture Notes in Statistics. Vol. 180. Springer; New York: 2005. Inference for Change Point and Post Change Means After a CUSUM Test.