Introduction
Detecting structural changes in time-series is important across many areas such as climate change [16 ] , environmental science [10 ] , finance [2 ] , bioinformatics [11 ] , medical imaging [27 ] , software reliability (Ke and Huang, 2020) and physics [29 ] , amongst many others. See Li et al. (2026) for a recent review of the area and applications. Whilst traditionally the focus of methods to detect these changepoints has been on offline analysis, there has been increasing interest in online methods for detecting changepoints due to the importance of detecting changes in real-time, from data such as medical monitoring signals [15 ] , telecommunications network data [21 ] , and industrial sensor measurements [32 ] .
Frequentist approaches to online changepoint detection are well established and include approaches based on sequential likelihood ratio tests [18 , 17 , 22 ] (Yu et al. 2023) and moving window approaches [13 ] . However, for some applications there are advantages of using a Bayesian approach. Such methods can incorporate prior information about the type and size of change, and they directly quantify uncertainty in both the presence and location of a change. For offline analysis, Bayesian changepoint methods allow inference for multiple changes and tend to involve Markov Chain Monte Carlo [26 , 9 , 5 , 4 ] . For a certain class of models, exact calculation and simulation from the posterior is possible [3 , 6 , 8 ] through solving a set of recursions.
The latter class of recursive methods can be implemented in an online manner. These apply to models where the likelihood of the data after a change does not depend on the parameters associated with segments before that change. Exact inference has a computational cost that increases linearly with time, so processing the observation at time T T T O ( T ) O(T) O ( T ) O ( 1 ) O(1) O ( 1 ) [7 , 1 ] .
For online applications, we often only need to detect whether or not there has been a changepoint prior to the current time. This then involves only considering models for the data with no or just one change. In this paper we consider such a setting. We specifically consider models where the data can be written as a unknown signal plus noise, and where a change corresponds to a change in the structure of the signal. It is straightforward to construct a general recursion for calculating the posterior in such a model. Importantly, as we are only allowing at most one change, this recursion can be applied to a wider range of models than those of Fearnhead and Liu [7 ] and Adams and MacKay [1 ] as the data after the change can depend on the parameters associated with the pre-change data. This allows us to, for example, detect changes in slope and changes in mean in the presence of seasonal effects.
Solving our recursion exactly leads to a computational cost per iteration that increases linearly with the amount of data. To avoid this we could use the pruning (or resampling) ideas from Fearnhead and Liu [7 ] and Adams and MacKay [1 ] : the increasing computational cost comes from having to consider all possible changepoint locations prior to the current time, and pruning involves ignoring some of these possibilities which have lowest probability. We instead propose an alternative, novel approach to reducing the computational cost. This involves merging rather than pruning: we partition the set of possible changepoint locations into a set of consecutive intervals, and then jointly update the posterior probabilities associated with values within an interval. This has the advantage over pruning in that it retains an approximation to the posterior distribution of a change at all time-points, and we show that this leads to a more accurate approximation of the posterior of the change location.
Merging has previously been suggested, for a different class of changepoint models, by Shamp et al. [24 ] , though without the idea of approximating the posterior distribution on the time of a change that motivates our use. Our implementation of how we merge is also different. Merging involves two steps: deciding which two current intervals of changepoint times to merge, and then choosing an approximate posterior for the parameters of the model given a change occurring within the merged interval. The approximation of merging comes from assuming the same posterior for the parameters for each changepoint time in the interval – as the true posteriors will differ for different changepoints. Shamp et al. [24 ] choose to merge the two intervals with the lowest current posterior probability, and then approximate the posterior for this merged interval based on minimising some distance to the posteriors that were associated with the two intervals that have been merged. We instead use a criterion that aims to minimise the difference between the posteriors on the post-change parameters prior to and after merging. Importantly, we do not just merge intervals based on their current posterior probability, but also on how similar the posteriors associated with each interval are.
The intuition behind our merging criteria is that the future updates of the posterior probabilities associated with each changepoint location depend on the posterior for the post-change parameter given a change at that location. By merging intervals with similar posteriors for the post-change parameter, we ensure that future approximations to the posterior on the changepoint location are more accurate. We demonstrate empirically that this is the case, and that often we can get substantially more accurate results with our merging criteria than if we just merged based on the posterior probabilities associated with a change in each interval. We name our approach Bagel, from BAyesian online chanGEpoint detection for Linear models.
The outline of the paper is as follows. In Section 2, we define the Bayesian real-time changepoint detection problem and introduce our proposed algorithm. Sections 3, and 4 show the performance of our algorithm on simulated data and real datasets. More examples, proofs, and detailed simulation results can be found in Supplement A. Data and reproducible code can be found in Supplement B.
Throughout this paper, all vectors are assumed to be column vectors unless explicitly stated otherwise. We denote by 1 1 : t \mathbf{1}_{1:t} 1 1 : t ( 1 , ⋯ , t ) T (1,\cdots,t)^{\mathsf{T}} ( 1 , ⋯ , t ) T 1 1 : t \mathbf{1}_{1:t} 1 1 : t 0 1 : t \mathbf{0}_{1:t} 0 1 : t t t t y 1 : t \mathbf{y}_{1:t} y 1 : t ( y 1 , ⋯ , y t ) T (y_1,\cdots,y_t)^{\mathsf{T}} ( y 1 , ⋯ , y t ) T μ \mu μ σ 2 \sigma^2 σ 2 N ( x ; μ , σ 2 ) N(x;\mu,\sigma^2) N ( x ; μ , σ 2 ) ν \nu ν l l l e e e t ν ( x ; l , e ) t_\nu(x;l,e) t ν ( x ; l , e ) ν \nu ν ι \iota ι I G ( x ; ν , ι ) IG(x;\nu,\iota) I G ( x ; ν , ι ) A A A A i : j × k : l A_{i:j\times k:l} A i : j × k : l i i i j j j k k k l l l
Univariate real-time Bayesian changepoint detection
The changepoint problem
The model
Assume an independent data stream y 1 , y 2 , … , y t y_1,y_2,\ldots,y_t y 1 , y 2 , … , y t t t t t t t
Let τ \tau τ τ = ∞ \tau= \infty τ = ∞ y t y_t y t t t t
y t = a t T β + b t , τ T γ + σ ϵ t . y_t = a_t^{\mathsf{T}}\beta+ b_{t,\tau}^{\mathsf{T}}\gamma+ \sigma\epsilon_t. y t = a t T β + b t , τ T γ + σ ϵ t . Here a t a_t a t b t , τ b_{t,\tau} b t , τ d 1 d_1 d 1 d 2 d_2 d 2 β \beta β γ \gamma γ d 1 d_1 d 1 d 2 d_2 d 2 a t T β a_t^{\mathsf{T}}\beta a t T β t t t b t , τ T γ b_{t,\tau}^{\mathsf{T}}\gamma b t , τ T γ b t , τ = 0 1 : d 2 b_{t,\tau} = 0_{1:d_2} b t , τ = 0 1 : d 2 τ ≥ t \tau\geq t τ ≥ t ϵ t \epsilon_t ϵ t σ > 0 \sigma> 0 σ > 0
Consequently, at any time t t t t t t y 1 : t y_{1:t} y 1 : t
y 1 : t = A 1 : t β + B 1 : t , τ γ + σ ϵ 1 : t , y_{1:t} = A_{1:t}\beta+ B_{1:t,\tau}\gamma+ \sigma\epsilon_{1:t}, y 1 : t = A 1 : t β + B 1 : t , τ γ + σ ϵ 1 : t , where ϵ 1 : t \epsilon_{1:t} ϵ 1 : t N ( 0 , 1 ) N(0,1) N ( 0 , 1 ) t t t A 1 : t A_{1:t} A 1 : t B 1 : t , τ B_{1:t,\tau} B 1 : t , τ a t T a_t^{\mathsf{T}} a t T b t , τ T b_{t,\tau}^{\mathsf{T}} b t , τ T τ ≥ t \tau\geq t τ ≥ t B 1 : t , τ B_{1:t,\tau} B 1 : t , τ τ \tau τ t t t t t t τ = 1 , 2 , … , t − 1 \tau= 1,2,\ldots,t-1 τ = 1 , 2 , … , t − 1
The prior for the model parameters
The likelihood of the observation depends on the parameters ( β , γ ) (\beta,\gamma) ( β , γ ) ( β , γ , σ ) (\beta,\gamma,\sigma) ( β , γ , σ )
For our model, the conjugate priors are Gaussian for ( β , γ ) (\beta,\gamma) ( β , γ ) σ \sigma σ σ \sigma σ β \beta β σ \sigma σ γ \gamma γ τ \tau τ β \beta β ( β , γ ) (\beta,\gamma) ( β , γ ) σ \sigma σ τ \tau τ
β , γ ∣ σ 2 , τ ∼ N ( μ τ , σ 2 Σ τ ) \beta,\gamma\mid\sigma^2,\tau\sim N(\mu^\tau,\sigma^2\Sigma^\tau) β , γ ∣ σ 2 , τ ∼ N ( μ τ , σ 2 Σ τ ) with prior mean μ τ ∈ R d 1 + d 2 \mu^\tau\in\mathbb{R}^{d_1+d_2} μ τ ∈ R d 1 + d 2 Σ τ ∈ R ( d 1 + d 2 ) × ( d 1 + d 2 ) \Sigma^\tau\in\mathbb{R}^{(d_1+d_2)\times(d_1+d_2)} Σ τ ∈ R ( d 1 + d 2 ) × ( d 1 + d 2 ) β \beta β σ \sigma σ μ 1 : d 1 τ \mu^\tau_{1:d_1} μ 1 : d 1 τ Σ 1 : d 1 × 1 : d 1 τ \Sigma^\tau_{1:d_1\times1:d_1} Σ 1 : d 1 × 1 : d 1 τ β \beta β τ \tau τ σ \sigma σ
σ 2 ∼ I G ( ν , ι ) . \sigma^2 \sim IG(\nu,\iota). σ 2 ∼ I G ( ν , ι ) . While it is often possible to use improper priors for β \beta β σ \sigma σ γ \gamma γ t t t γ \gamma γ
Here we give three examples to illustrate our model framework (with derivations of the stated priors given in Appendix § 6 ). In each case we present the priors for the case where the noise variance is known, but extending them to the unknown variance case only requires specifying the inverse Gamma prior for the variance. Through we denote user-chosen constants in the prior mean by μ i \mu_i μ i i = 1 , 2 , … i=1,2,\ldots i = 1 , 2 , … δ i > 0 \delta_i > 0 δ i > 0 i = 1 , 2 , … i=1,2,\ldots i = 1 , 2 , …
Example 1 (Change-in-mean model): We first consider the change-in-mean model where d 1 = d 2 = 1 d_1=d_2=1 d 1 = d 2 = 1
A 1 : t = 1 1 : t , B 1 : t , τ = [ 0 1 : τ 1 τ + 1 : t ] A_{1:t}=1_{1:t}, B_{1:t,\tau} =
\begin{bmatrix}
0_{1:\tau}\\
1_{\tau+1:t}
\end{bmatrix} A 1 : t = 1 1 : t , B 1 : t , τ = [ 0 1 : τ 1 τ + 1 : t ] Here β \beta β γ \gamma γ
distribution in this setting: one where the pre-change mean and the change in mean are independent, and another where the pre-change and post-change means are independent and identically distributed.
For the former case we have μ = [ μ 1 0 ] \mu= \begin{bmatrix} \mu_1 \\ 0 \end{bmatrix} μ = [ μ 1 0 ] Σ = [ δ 1 0 0 δ 2 ] \Sigma= \begin{bmatrix} \delta_1 & 0 \\ 0 & \delta_2 \end{bmatrix} Σ = [ δ 1 0 0 δ 2 ] Σ = [ δ 1 − δ 1 − δ 1 2 δ 1 ] \Sigma= \begin{bmatrix} \delta_1 & -\delta_1 \\ -\delta_1 & 2\delta_1 \end{bmatrix} Σ = [ δ 1 − δ 1 − δ 1 2 δ 1 ]
Example 2 (Change-in-slope models): Next, we consider the change-in-slope where d 1 = d 2 = 2 d_1 = d_2 = 2 d 1 = d 2 = 2
A 1 : t = [ 1 1 : t 1 : t ] , B 1 : t , τ = [ 0 1 : τ 0 1 : τ 1 ( τ + 1 ) : t ( τ + 1 ) : t ] A_{1:t} = [1_{1:t}\ \ 1:t], \quad B_{1:t,\tau} = \begin{bmatrix} 0_{1:\tau} & 0_{1:\tau} \\ 1_{(\tau+1):t} & (\tau+1):t \end{bmatrix} A 1 : t = [ 1 1 : t 1 : t ] , B 1 : t , τ = [ 0 1 : τ 1 ( τ + 1 ) : t 0 1 : τ ( τ + 1 ) : t ] Here β \beta β γ \gamma γ
There are two natural types of change we can allow in this model, depending on whether or not we require continuity of the trend line at the changepoint, as illustrated in Figure 1 . If we wish to detect a continuous change, also known as a change-in-slope ([8 ] ), then a natural model is that the change-in-slope, γ 2 \gamma_2 γ 2 γ 1 \gamma_1 γ 1 γ 2 \gamma_2 γ 2 τ \tau τ τ \tau τ γ 1 + τ γ 2 = 0 \gamma_1 + \tau\gamma_2 = 0 γ 1 + τ γ 2 = 0 ( β , γ ) (\beta,\gamma) ( β , γ )
Figure 1. An example of continuous change (left) and discontinuous change (right) in a linear trend model.
μ τ = ( μ 1 , μ 2 , 0 , 0 ) T , and Σ τ = [ δ 1 0 0 0 0 δ 2 0 0 0 0 τ 2 δ 3 − τ δ 3 0 0 − τ δ 3 δ 3 ] . \mu^\tau= (\mu_1,\mu_2,0,0)^\mathsf{T}, \quad\text{and} \quad\Sigma^\tau= \begin{bmatrix} \delta_1 & 0 & 0 & 0 \\ 0 & \delta_2 & 0 & 0 \\ 0 & 0 & \tau^2\delta_3 & -\tau\delta_3 \\ 0 & 0 & -\tau\delta_3 & \delta_3 \end{bmatrix}. μ τ = ( μ 1 , μ 2 , 0 , 0 ) T , and Σ τ = δ 1 0 0 0 0 δ 2 0 0 0 0 τ 2 δ 3 − τ δ 3 0 0 − τ δ 3 δ 3 . For a discontinuous change a natural choice of prior would assume that the distribution of the trend line starting at time 0 is independent of and has the same distribution of the trend line starting at time τ \tau τ τ \tau τ β 1 + γ 1 + τ ( β 2 + γ 2 ) \beta_1+\gamma_1+\tau(\beta_2+\gamma_2) β 1 + γ 1 + τ ( β 2 + γ 2 ) β 1 \beta_1 β 1 ( β 2 + γ 2 ) (\beta_2+\gamma_2) ( β 2 + γ 2 ) β 2 \beta_2 β 2
μ τ = ( μ 1 , μ 2 , − μ 2 τ , 0 ) T , and Σ τ = [ δ 1 0 − δ 1 0 0 δ 2 0 − δ 2 − δ 1 0 2 δ 1 + τ 2 δ 2 − τ δ 2 0 − δ 2 − τ δ 2 2 δ 2 ] . \mu^\tau= (\mu_1,\mu_2,-\mu_2\tau,0)^\mathsf{T}, \quad\text{and} \quad\Sigma^\tau= \begin{bmatrix} \delta_1 & 0 & -\delta_1 & 0 \\ 0 & \delta_2 & 0 & -\delta_2 \\ -\delta_1 & 0 & 2\delta_1+\tau^2\delta_2 & -\tau\delta_2 \\ 0 & -\delta_2 & -\tau\delta_2 & 2\delta_2 \end{bmatrix}. μ τ = ( μ 1 , μ 2 , − μ 2 τ , 0 ) T , and Σ τ = δ 1 0 − δ 1 0 0 δ 2 0 − δ 2 − δ 1 0 2 δ 1 + τ 2 δ 2 − τ δ 2 0 − δ 2 − τ δ 2 2 δ 2 . Example 3 (Change-in-seasonality model): The framework can include seasonality. For simplicity, we consider the seasonal model with one Fourier-harmonic seasonality period P P P d 1 = d 2 = 2 d_1=d_2=2 d 1 = d 2 = 2
A 1 : t = [ sin 2 π ( 1 : t ) P , cos 2 π ( 1 : t ) P ] , B 1 : t , τ = [ 0 1 : τ 0 1 : τ sin 2 π ( ( τ + 1 ) : t ) P cos 2 π ( ( τ + 1 ) : t ) P ] . A_{1:t}=\left[\sin\frac{2\pi(1:t)}{P},\ \cos\frac{2\pi(1:t)}{P}\right],\quad B_{1:t,\tau}=\left[\begin{array}{cc}0_{1:\tau}&0_{1:\tau}\\ \sin\frac{2\pi((\tau+1):t)}{P}&\cos\frac{2\pi((\tau+1):t)}{P}\end{array}\right]. A 1 : t = [ sin P 2 π ( 1 : t ) , cos P 2 π ( 1 : t ) ] , B 1 : t , τ = [ 0 1 : τ sin P 2 π (( τ + 1 ) : t ) 0 1 : τ cos P 2 π (( τ + 1 ) : t ) ] . Here β \beta β γ \gamma γ
A natural choice of prior would be for the distribution of the Fourier coefficients before and after the change to be independent. In this case the mean will be μ T = [ μ sin , μ cos , 0 , 0 ] \mu^{T}=[\mu_{\sin},\mu_{\cos},0,0] μ T = [ μ s i n , μ c o s , 0 , 0 ] ( β 1 , β 2 ) T (\beta_1,\beta_2)^{T} ( β 1 , β 2 ) T
Σ s = [ δ sin 0 0 δ cos ] , \Sigma_{\mathrm{s}}=\begin{bmatrix}\delta_{\sin}&0\\0&\delta_{\cos}\end{bmatrix}, Σ s = [ δ s i n 0 0 δ c o s ] , then the prior variance will be
Σ τ = [ Σ s − Σ s − Σ s 2 Σ s ] . \Sigma^{\tau}=\begin{bmatrix}\Sigma_{\mathrm{s}}&-\Sigma_{\mathrm{s}}\\-\Sigma_{\mathrm{s}}&2\Sigma_{\mathrm{s}}\end{bmatrix}. Σ τ = [ Σ s − Σ s − Σ s 2 Σ s ] . Prior for the changepoint
To complete our model specification we need a prior for τ \tau τ t t t Pr t ( τ = i ) \Pr_t(\tau=i) Pr t ( τ = i ) i = 1 , … , t − 1 i=1,\ldots,t-1 i = 1 , … , t − 1 Pr t ( τ ≥ t ) \Pr_t(\tau\ge t) Pr t ( τ ≥ t ) t t t
One approach to defining the prior at time t t t τ \tau τ 1 , … 1,\ldots 1 , … τ = ∞ \tau=\infty τ = ∞ Pr t ( τ ) \Pr_t(\tau) Pr t ( τ ) t t t τ \tau τ τ \tau τ t t t
Thus we suggest the following approach. First let us define the Bayes factor for a change at τ = i \tau= i τ = i τ ≥ t \tau\ge t τ ≥ t
B F t ( i ) = p ( y 1 : t ∣ τ = i ) p ( y 1 : t ∣ τ ≥ t ) , BF_t(i) = \frac{p(y_{1:t} \mid\tau= i)}{p(y_{1:t} \mid\tau\ge t)}, B F t ( i ) = p ( y 1 : t ∣ τ ≥ t ) p ( y 1 : t ∣ τ = i ) , where the terms on the right-hand side are marginal likelihoods, defined as
p ( y 1 : t ∣ τ ) = ∬ p ( y 1 : t ∣ τ , β , γ , σ ) p ( β , γ ∣ τ , σ ) d β d γ , p(y_{1:t} \mid\tau) = \iint p(y_{1:t} \mid\tau, \beta, \gamma, \sigma)p(\beta, \gamma\mid\tau, \sigma)\,\mathrm{d}\beta\mathrm{d}\gamma, p ( y 1 : t ∣ τ ) = ∬ p ( y 1 : t ∣ τ , β , γ , σ ) p ( β , γ ∣ τ , σ ) d β d γ , in the case where σ \sigma σ σ \sigma σ t t t τ \tau τ t t t
B F t = ∑ i = 1 t − 1 B F t ( i ) Pr t ( τ = i ∣ τ < t ) . BF_t = \sum_{i=1}^{t-1} BF_t(i)\Pr_t(\tau= i \mid\tau< t). B F t = i = 1 ∑ t − 1 B F t ( i ) t Pr ( τ = i ∣ τ < t ) . So the posterior probability of a change before time t t t
p t ( τ < t ∣ y 1 : t ) = ( 1 − Pr t ( τ ≥ t ) ) B F t Pr t ( τ ≥ t ) + ( 1 − Pr t ( τ ≥ t ) ) B F t . p_t(\tau< t \mid y_{1:t}) = \frac{(1-\Pr_t(\tau\ge t))BF_t}{\Pr_t(\tau\ge t) + (1-\Pr_t(\tau\ge t))BF_t}. p t ( τ < t ∣ y 1 : t ) = Pr t ( τ ≥ t ) + ( 1 − Pr t ( τ ≥ t )) B F t ( 1 − Pr t ( τ ≥ t )) B F t . If we use a monitoring scheme which detects a change if p t ( τ < t ∣ y 1 : t ) ≥ c t p_t(\tau< t \mid y_{1:t}) \ge c_t p t ( τ < t ∣ y 1 : t ) ≥ c t c t c_t c t
B F t ≥ c t Pr t ( τ ≥ t ) ( 1 − c t ) ( 1 − Pr t ( τ ≥ t ) ) . BF_t \ge\frac{c_t\Pr_t(\tau\ge t)}{(1-c_t)(1-\Pr_t(\tau\ge t))}. B F t ≥ ( 1 − c t ) ( 1 − Pr t ( τ ≥ t )) c t Pr t ( τ ≥ t ) . In the absence of specific prior information, it would be natural to detect a change based on the Bayes Factor being above a threshold that is constant over time. This is easiest to implement with c t c_t c t Pr t ( τ ≥ t ) \Pr_t(\tau\ge t) Pr t ( τ ≥ t ) Pr t ( τ ≥ t ) \Pr_t(\tau\ge t) Pr t ( τ ≥ t )
The test does depend on our prior for τ \tau τ τ < t \tau< t τ < t
Pr t ( τ = t − k ∣ τ < t ) ∝ p k , for k = 1 , … , t − 1. \Pr_t(\tau= t-k \mid\tau< t) \propto p^k,\quad\text{for } k=1,\ldots,t-1. t Pr ( τ = t − k ∣ τ < t ) ∝ p k , for k = 1 , … , t − 1. This is a (truncated geometric) prior for the time since the change. This prior is time dependent, and is updated sequentially as new observations arrive, ensuring that Pr t ( τ < t ) \Pr_t(\tau< t) Pr t ( τ < t ) p = 1 p=1 p = 1 [20 ] .
Sequential Updating
Recall that we wish to calculate the posterior probability of a change, its location and the parameters of the model at each time t t t t − 1 t-1 t − 1 t t t y t y_t y t
The posterior probability for τ \tau τ t t t
Pr t ( τ ∣ y 1 : t ) ∝ Pr t ( τ ) p ( y 1 : t ∣ τ ) . \Pr_t(\tau\mid y_{1:t}) \propto\Pr_t(\tau)p(y_{1:t}\mid\tau). t Pr ( τ ∣ y 1 : t ) ∝ t Pr ( τ ) p ( y 1 : t ∣ τ ) . We will derive a recursion for the right-hand side. To do this at time t t t w i , t = Pr t ( τ = i ) p ( y 1 : t ∣ τ = i ) w_{i,t}=\Pr_t(\tau=i)p(y_{1:t}\mid\tau=i) w i , t = Pr t ( τ = i ) p ( y 1 : t ∣ τ = i ) i = 1 , … , t − 1 i=1,\ldots,t-1 i = 1 , … , t − 1 w 0 , t = Pr t ( τ ≥ t ) p ( y 1 : t ∣ τ ≥ t ) w_{0,t}=\Pr_t(\tau\ge t)p(y_{1:t}\mid\tau\ge t) w 0 , t = Pr t ( τ ≥ t ) p ( y 1 : t ∣ τ ≥ t ) i , t i,t i , t t t t i i i τ \tau τ i = 0 i=0 i = 0 t t t
The following result gives the update for these weights in terms of the predictive density for y t y_t y t y 1 : t − 1 y_{1:t-1} y 1 : t − 1
Theorem 2.1. For t > 1 t > 1 t > 1
w i , t = w i , t − 1 Pr t ( τ = i ) Pr t − 1 ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) for i = 1 , … , t − 2 , and t > 2 w_{i,t}=w_{i,t-1}\frac{\Pr_t(\tau=i)}{\Pr_{t-1}(\tau=i)}p(y_t\mid y_{1:t-1},\tau=i)\quad\text{for }i=1,\ldots,t-2,\text{ and }t>2 w i , t = w i , t − 1 Pr t − 1 ( τ = i ) Pr t ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) for i = 1 , … , t − 2 , and t > 2 w 0 , t = w 0 , t − 1 p ( y t ∣ y 1 : t − 1 , τ ≥ t ) , and w_{0,t}=w_{0,t-1}p(y_t\mid y_{1:t-1},\tau\geq t),\quad\text{and} w 0 , t = w 0 , t − 1 p ( y t ∣ y 1 : t − 1 , τ ≥ t ) , and w t − 1 , t = w 0 , t − 1 Pr t ( τ = t − 1 ) Pr t − 1 ( τ ≥ t − 1 ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) . w_{t-1,t}=w_{0,t-1}\frac{\Pr_t(\tau=t-1)}{\Pr_{t-1}(\tau\geq t-1)}p(y_t\mid y_{1:t-1},\tau=t-1). w t − 1 , t = w 0 , t − 1 Pr t − 1 ( τ ≥ t − 1 ) Pr t ( τ = t − 1 ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) . With the weights initialised at t = 1 t=1 t = 1
w 0 , 1 = Pr t ( τ ≥ 1 ) p ( y 1 ∣ τ ≥ 1 ) . w_{0,1}=\Pr_t(\tau\geq1)p(y_1\mid\tau\geq1). w 0 , 1 = t Pr ( τ ≥ 1 ) p ( y 1 ∣ τ ≥ 1 ) . Proof. See Appendix B
We can implement these recursions by recursively calculating the posterior of the model parameters, and hence the predictive distribution, for each possible value of τ \tau τ σ \sigma σ σ \sigma σ σ \sigma σ
As the prior for the post-change parameter can depend on the time of change, in the case of no changepoint we will update only the marginal posterior for β \beta β t − 1 t-1 t − 1 β \beta β σ \sigma σ
p ( β ∣ y 1 : t − 1 , σ , τ ≥ t − 1 ) ∼ N ( μ t − 1 , 0 , σ 2 Σ t − 1 , 0 ) , p(\beta\mid y_{1:t-1},\sigma,\tau\geq t-1)\sim N\left(\mu^{t-1,0},\sigma^2\Sigma^{t-1,0}\right), p ( β ∣ y 1 : t − 1 , σ , τ ≥ t − 1 ) ∼ N ( μ t − 1 , 0 , σ 2 Σ t − 1 , 0 ) , and
p ( σ 2 ∣ τ ≥ t − 1 ) ∼ I G ( ν t − 1 , 0 , ι t − 1 , 0 ) , p(\sigma^2\mid\tau\geq t-1)\sim IG(\nu^{t-1,0},\iota^{t-1,0}), p ( σ 2 ∣ τ ≥ t − 1 ) ∼ I G ( ν t − 1 , 0 , ι t − 1 , 0 ) , where μ t − 1 , 0 \mu^{t-1,0} μ t − 1 , 0 d 1 d_1 d 1 Σ t − 1 , 0 \Sigma^{t-1,0} Σ t − 1 , 0 d 1 × d 1 d_1\times d_1 d 1 × d 1 i = 1 , … , t − 2 i=1,\ldots,t-2 i = 1 , … , t − 2
p ( β , γ ∣ y 1 : t − 1 , σ , τ = i ) ∼ N ( μ t − 1 , i , σ 2 Σ t − 1 , i ) , p(\beta,\gamma\mid y_{1:t-1},\sigma,\tau=i)\sim N\left(\mu^{t-1,i},\sigma^2\Sigma^{t-1,i}\right), p ( β , γ ∣ y 1 : t − 1 , σ , τ = i ) ∼ N ( μ t − 1 , i , σ 2 Σ t − 1 , i ) , and
p ( σ 2 ∣ τ = i ) ∼ I G ( ν t − 1 , i , ι t − 1 , i ) , p(\sigma^2\mid\tau=i)\sim IG(\nu^{t-1,i},\iota^{t-1,i}), p ( σ 2 ∣ τ = i ) ∼ I G ( ν t − 1 , i , ι t − 1 , i ) , where now μ t − 1 , i \mu^{t-1,i} μ t − 1 , i d 1 + d 2 d_1+d_2 d 1 + d 2 Σ t − 1 , i \Sigma^{t-1,i} Σ t − 1 , i ( d 1 + d 2 ) × ( d 1 + d 2 ) (d_1+d_2)\times(d_1+d_2) ( d 1 + d 2 ) × ( d 1 + d 2 ) ( β , γ ) (\beta,\gamma) ( β , γ ) σ \sigma σ t − 1 t-1 t − 1 β \beta β σ \sigma σ t − 1 t-1 t − 1
Theorem 2.2. Let the prior mean for a change at t − 1 t-1 t − 1 μ t − 1 = [ μ β T , μ γ T ] T \boldsymbol{\mu}^{t-1} = [\mu_{\beta}^{\mathsf{T}},\mu_{\gamma}^{\mathsf{T}}]^{\mathsf{T}} μ t − 1 = [ μ β T , μ γ T ] T
σ 2 Σ t − 1 = σ 2 [ Σ β , β Σ β , γ Σ γ , β Σ γ , γ ] \sigma^{2}\boldsymbol{\Sigma}^{t-1} = \sigma^{2}
\begin{bmatrix}
\Sigma_{\beta,\beta} & \Sigma_{\beta,\gamma}\\
\Sigma_{\gamma,\beta} & \Sigma_{\gamma,\gamma}
\end{bmatrix} σ 2 Σ t − 1 = σ 2 [ Σ β , β Σ γ , β Σ β , γ Σ γ , γ ] so that μ β \mu_{\beta} μ β β \beta β Σ β , β \Sigma_{\beta,\beta} Σ β , β β \beta β Σ β , γ \Sigma_{\beta,\gamma} Σ β , γ β \beta β γ \gamma γ ( β , γ , σ ) (\beta,\gamma,\sigma) ( β , γ , σ ) y 1 : t − 1 y_{1:t-1} y 1 : t − 1 τ = t − 1 \tau=t-1 τ = t − 1
μ t − 1 , t − 1 = [ μ t − 1 , 0 μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) ] \boldsymbol{\mu}^{t-1,t-1} =
\begin{bmatrix}
\boldsymbol{\mu}^{t-1,0}\\
\mu_{\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}(\mu^{t-1,0}-\mu_{\beta})
\end{bmatrix} μ t − 1 , t − 1 = [ μ t − 1 , 0 μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) ] Σ t − 1 , t − 1 = [ Σ t − 1 , 0 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ] \boldsymbol{\Sigma}^{t-1,t-1} =
\begin{bmatrix}
\Sigma^{t-1,0} & \Sigma^{t-1,0}\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}\\
\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma^{t-1,0} & \Sigma_{\gamma,\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}(\Sigma^{t-1,0}-\Sigma_{\beta,\beta})\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}
\end{bmatrix} Σ t − 1 , t − 1 = [ Σ t − 1 , 0 Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ] with ν t − 1 , t − 1 = ν t − 1 , 0 \nu^{t-1,t-1}=\nu^{t-1,0} ν t − 1 , t − 1 = ν t − 1 , 0 ι t − 1 , t − 1 = ι t − 1 , 0 \iota^{t-1,t-1}=\iota^{t-1,0} ι t − 1 , t − 1 = ι t − 1 , 0
Proof. See Appendix B
We can update these parameters as follows.
Theorem 2.3. Fix iteration t − 1 t-1 t − 1 h i h_i h i
h i = { [ a t − 1 T , b t − 1 , i T ] T 0 < i ≤ t − 2 , a t − 1 i = 0. h_i =
\begin{cases}
[a_{t-1}^{\mathsf{T}},b_{t-1,i}^{\mathsf{T}}]^{\mathsf{T}} & 0<i\le t-2,\\
a_{t-1} & i=0.
\end{cases} h i = { [ a t − 1 T , b t − 1 , i T ] T a t − 1 0 < i ≤ t − 2 , i = 0. For i = 0 , 1 , … , t − 2 i=0,1,\ldots,t-2 i = 0 , 1 , … , t − 2 e i = y t − 1 − h i μ t − 1 , i e_i=y_{t-1}-h_i\boldsymbol{\mu}^{t-1,i} e i = y t − 1 − h i μ t − 1 , i Q = h i T Σ t − 1 , i h i + 1 Q=h_i^{\mathsf{T}}\boldsymbol{\Sigma}^{t-1,i}h_i+1 Q = h i T Σ t − 1 , i h i + 1 A = Σ t − 1 , i h i / Q A=\boldsymbol{\Sigma}^{t-1,i}h_i/Q A = Σ t − 1 , i h i / Q y t − 1 y_{t-1} y t − 1
Σ t , i = Σ t − 1 , i − A A T Q , μ t , i = μ t − 1 , i + A e , \boldsymbol{\Sigma}^{t,i}=\boldsymbol{\Sigma}^{t-1,i}-AA^{\mathsf{T}}Q,\qquad\boldsymbol{\mu}^{t,i}=\boldsymbol{\mu}^{t-1,i}+Ae, Σ t , i = Σ t − 1 , i − A A T Q , μ t , i = μ t − 1 , i + A e , ν t , i = ν t − 1 , i + 1 2 , ι t , i = ι t − 1 , i + 1 2 e i 2 / Q . \nu^{t,i}=\nu^{t-1,i}+\frac{1}{2},\qquad\iota^{t,i}=\iota^{t-1,i}+\frac{1}{2}e_i^2/Q. ν t , i = ν t − 1 , i + 2 1 , ι t , i = ι t − 1 , i + 2 1 e i 2 / Q . Finally, given the posterior for the parameters at time t − 1 t-1 t − 1 t t t
Theorem 2.4. Using the same notation for the parameters of the posteriors as in Theorem 2.3 , the predictive density can be calculated as
P ( y t ∣ y 1 : t − 1 , τ = i ) = N ( y t ; h i μ t , i , σ 2 ( 1 + h i T Σ t , i h i ) ) P(y_t\mid y_{1:t-1},\tau=i)=N\left(y_t;h_i\boldsymbol{\mu}^{t,i},\sigma^2(1+h_i^{\mathsf{T}}\boldsymbol{\Sigma}^{t,i}h_i)\right) P ( y t ∣ y 1 : t − 1 , τ = i ) = N ( y t ; h i μ t , i , σ 2 ( 1 + h i T Σ t , i h i ) ) if σ \sigma σ
P ( y t ∣ y 1 : t − 1 , τ = i ) = t 2 ν t , i ( y t ; h i μ t , i , t t , i ν t , i ( 1 + h i T Σ t , i h i ) ) P(y_t|y_{1:t-1}, \tau=i)=t_{2\nu_{t,i}}\left(y_t;h_i\mu^{t,i},\frac{t^{t,i}}{\nu_{t,i}}\left(1+h_i^{\mathsf{T}}\Sigma^{t,i}h_i\right)\right) P ( y t ∣ y 1 : t − 1 , τ = i ) = t 2 ν t , i ( y t ; h i μ t , i , ν t , i t t , i ( 1 + h i T Σ t , i h i ) ) if σ \sigma σ t ν ( y t ; l , e ) t_{\nu}(y_t;l,e) t ν ( y t ; l , e ) t t t ν \nu ν l l l e e e
Reducing the computational complexity by merging
Implementing the univariate Bayesian changepoint detection methods gives a cost per iteration that increases with iteration t t t T T T T T T M M M [1 , 23 ] ), in which the computational complexity can be reduced to O ( T M ) O(TM) O ( T M ) Figure 2 .
Figure 2. A failure case example for LP-P. Figure (a) presents the simulated data with a change at 1500. Figures (b) and (c) show the posterior distribution obtained from the exact approach without pruning (grey line). The black dots represent the candidates not pruned by the pruning approach with particle budget M = 50 M = 50 M = 50
Merging criteria
To overcome these issues, we first use the idea of merging the posterior distribution for the parameters for different change locations as a way to reduce the computational cost whilst keeping the full support of the posterior. Second, we take account of the similarity of the posterior distributions of the parameters when deciding which distributions to merge. We will look at the similarity of the posteriors in terms of the parameters, or the function of parameters, that affect the distribution of the data after the change. Denote such parameters as θ = g ( β , γ ) \theta= g(\beta,\gamma) θ = g ( β , γ ) g ( ⋅ ) g(\cdot) g ( ⋅ ) ( β , γ ) (\beta,\gamma) ( β , γ ) θ \theta θ
To explain our approach, it is simplest to first consider what happens when we first merge posteriors. Fix the time, t t t f i ( θ ) f_i(\theta) f i ( θ ) i i i t − 1 t - 1 t − 1 { w i , t , f i ( θ ) } \{w_{i,t}, f_i(\theta)\} { w i , t , f i ( θ )} 1 ≤ i ≤ t − 1 1 \leq i \leq t - 1 1 ≤ i ≤ t − 1
Our merge procedure will replace f i ( θ ) f_i(\theta) f i ( θ ) f i + 1 ( θ ) f_{i+1}(\theta) f i + 1 ( θ ) i i i i ∈ { 1 , … , t − 2 } i \in\{1, \ldots, t - 2\} i ∈ { 1 , … , t − 2 } θ \theta θ T V ( h 1 ( θ ) , h 2 ( θ ) ) = 1 2 ∫ ∣ h 1 ( θ ) − h 2 ( θ ) ∣ d θ TV(h_1(\theta), h_2(\theta)) = \frac{1}{2} \int|h_1(\theta) - h_2(\theta)|\mathrm{d}\theta T V ( h 1 ( θ ) , h 2 ( θ )) = 2 1 ∫ ∣ h 1 ( θ ) − h 2 ( θ ) ∣ d θ h 1 ( θ ) h_1(\theta) h 1 ( θ ) h 2 ( θ ) h_2(\theta) h 2 ( θ ) h 1 ( θ ) = ∑ i = 1 t − 1 w i , t f i ( θ ) h_1(\theta) = \sum_{i=1}^{t-1} w_{i,t}f_i(\theta) h 1 ( θ ) = ∑ i = 1 t − 1 w i , t f i ( θ ) f i ( θ ) f_i(\theta) f i ( θ ) f i + 1 ( θ ) f_{i+1}(\theta) f i + 1 ( θ ) h 2 ( θ ) = h 1 ( θ ) − w i , t ( f i ( θ ) − f i + 1 ( θ ) ) h_2(\theta) = h_1(\theta) - w_{i,t}(f_i(\theta) - f_{i+1}(\theta)) h 2 ( θ ) = h 1 ( θ ) − w i , t ( f i ( θ ) − f i + 1 ( θ ))
T V ( h 1 ( θ ) , h 2 ( θ ) ) = 1 2 ∫ ∣ h 1 ( θ ) − ( h 1 ( θ ) − w i , t ( f i ( θ ) − f i + 1 ( θ ) ) ) ∣ d θ = 1 2 ∫ ∣ w i , t ( f i ( θ ) − f i + 1 ( θ ) ) ∣ d θ = w i , t T V ( f i ( θ ) , f i + 1 ( θ ) ) . \begin{aligned}
TV(h_1(\theta), h_2(\theta)) &= \frac{1}{2} \int\left|h_1(\theta) - \left(h_1(\theta) - w_{i,t}(f_i(\theta) - f_{i+1}(\theta))\right)\right|\mathrm{d}\theta\\
&= \frac{1}{2} \int\left|w_{i,t}(f_i(\theta) - f_{i+1}(\theta))\right|\mathrm{d}\theta\\
&= w_{i,t}TV(f_i(\theta), f_{i+1}(\theta)).
\end{aligned} T V ( h 1 ( θ ) , h 2 ( θ )) = 2 1 ∫ ∣ h 1 ( θ ) − ( h 1 ( θ ) − w i , t ( f i ( θ ) − f i + 1 ( θ )) ) ∣ d θ = 2 1 ∫ ∣ w i , t ( f i ( θ ) − f i + 1 ( θ )) ∣ d θ = w i , t T V ( f i ( θ ) , f i + 1 ( θ )) . Thus we will choose to merge component i i i
w i , t T V ( f i ( θ ) , f i + 1 ( θ ) ) . (1) w_{i,t}TV(f_i(\theta), f_{i+1}(\theta)). \tag*{(1)} w i , t T V ( f i ( θ ) , f i + 1 ( θ )) . (1) This will maintain the same support for the changepoint locations, but will reduce the computational cost. At any time, if we have M M M θ \theta θ
So at any iteration t t t θ \theta θ T = { j 1 , j 2 , … , j M } \mathcal{T} = \{j_1,j_2,\ldots,j_M\} T = { j 1 , j 2 , … , j M } w i , t w_{i,t} w i , t i = 0 , … , M i = 0,\ldots,M i = 0 , … , M i ∈ { 1 , … , M } i \in\{1,\ldots,M\} i ∈ { 1 , … , M } w i , t w_{i,t} w i , t { j i − 1 + 1 , … , j i } \{j_{i-1}+1,\ldots,j_i\} { j i − 1 + 1 , … , j i } j 0 = 0 j_0 = 0 j 0 = 0 M M M R = { r ( 1 ) , … , r ( M ) } \mathcal{R} = \{r^{(1)},\ldots,r^{(M)}\} R = { r ( 1 ) , … , r ( M ) } w i w_i w i j i − 1 + 1 j_{i-1}+1 j i − 1 + 1 j i j_i j i w i r ( i ) w_i r^{(i)} w i r ( i )
( w 0 , w 1 r ( 1 ) , … , w M r ( M ) ) , (w_0,w_1r^{(1)},\ldots,w_Mr^{(M)}), ( w 0 , w 1 r ( 1 ) , … , w M r ( M ) ) , with, after normalisation, the first component being the probability of no change, and the remaining components being the probability of a change at 1 , … , t − 1 1,\ldots,t-1 1 , … , t − 1
The idea is that we can update this representation by just updating the weights, and the parameters associated with each distinct posterior – which has a computational cost proportional to M M M M M M θ \theta θ t > M t > M t > M f i ( θ ) f_i(\theta) f i ( θ ) i i i f i ( θ ) f_i(\theta) f i ( θ ) f i + 1 ( θ ) f_{i+1}(\theta) f i + 1 ( θ ) i ≥ 1 i \ge1 i ≥ 1 w i TV ( f i ( θ ) , f i + 1 ( θ ) ) w_i\operatorname{TV}(f_i(\theta), f_{i+1}(\theta)) w i TV ( f i ( θ ) , f i + 1 ( θ ))
Calculating the total variation distance
Now we will show how to calculate the total variation, which is needed for our merging step. This can be calculated exactly if θ \theta θ
One dimensional parameter
Algorithm 1: one iteration of Bagel
Input: existing candidates set T = { 0 , i 1 , … , i N } \mathcal{T} = \{0, i_1, \ldots, i_N\} T = { 0 , i 1 , … , i N } { r ( 1 ) , … , r ( N ) } ∗ \{r^{(1)}, \ldots, r^{(N)}\}^* { r ( 1 ) , … , r ( N ) } ∗ { w i , t − 1 } \{w_{i,t-1}\} { w i , t − 1 } Θ i \Theta_i Θ i i ∈ { 0 , … , N } i \in\{0, \ldots, N\} i ∈ { 0 , … , N }
Data: Observe y t y_t y t t t t
Introducing new changepoint candidate
Calculate Θ N + 1 \Theta_{N+1} Θ N + 1 τ = t − 1 \tau= t - 1 τ = t − 1 Theorem 2.2 .
T ← T ∪ { t − 1 } \mathcal{T} \leftarrow\mathcal{T} \cup\{t - 1\} T ← T ∪ { t − 1 }
Updating weights and posterior distributions
for i ∈ { 0 , … , N + 1 } i \in\{0, \ldots, N + 1\} i ∈ { 0 , … , N + 1 } do
Calculate w i , t w_{i,t} w i , t Theorems 2.1 and 2.4 .
Update parameters of posterior distribution, Θ 0 , … , Θ N + 1 \Theta_0, \ldots, \Theta_{N+1} Θ 0 , … , Θ N + 1 Theorem 2.3 .
Normalize w i , t w_{i,t} w i , t i ∈ { 0 , … , N + 1 } i \in\{0, \ldots, N + 1\} i ∈ { 0 , … , N + 1 }
Decide if there is a changepoint
if 1 − w 0 , t > c 1 - w_{0,t} > c 1 − w 0 , t > c then
output: Change detected; T \mathcal{T} T { r ( 1 ) , … , r ( N + 1 ) } ∗ \{r^{(1)}, \ldots, r^{(N+1)}\}^* { r ( 1 ) , … , r ( N + 1 ) } ∗ { Θ 0 , … , Θ N + 1 } \{\Theta_0, \ldots, \Theta_{N+1}\} { Θ 0 , … , Θ N + 1 } w 0 : N + 1 , t w_{0:N+1,t} w 0 : N + 1 , t
Merging step
if N + 1 > M N + 1 > M N + 1 > M then
Search the index i i i i ∈ { 1 , … , N } i \in\{1, \ldots, N\} i ∈ { 1 , … , N }
*Update the ratio vector { r ( i + 1 ) } ← ( w i , t + w i + 1 , t ) − 1 ( w i , t r ( i ) , w i + 1 , t r ( i + 1 ) ) \{r^{(i+1)}\} \leftarrow(w_{i,t} + w_{i+1,t})^{-1}(w_{i,t}r^{(i)}, w_{i+1,t}r^{(i+1)}) { r ( i + 1 ) } ← ( w i , t + w i + 1 , t ) − 1 ( w i , t r ( i ) , w i + 1 , t r ( i + 1 ) )
Combine w i + 1 , t ← w i , t + w i + 1 , t w_{i+1,t} \leftarrow w_{i,t} + w_{i+1,t} w i + 1 , t ← w i , t + w i + 1 , t
Remove candidate i i i Θ i \Theta_i Θ i r ( i ) r^{(i)} r ( i ) w i , t w_{i,t} w i , t
N ← N − 1 N \leftarrow N - 1 N ← N − 1
output: T \mathcal{T} T { r ( 1 ) , … , r ( N + 1 ) } ∗ \{r^{(1)}, \ldots, r^{(N+1)}\}^* { r ( 1 ) , … , r ( N + 1 ) } ∗ { Θ 0 , … , Θ N + 1 } \{\Theta_0, \ldots, \Theta_{N+1}\} { Θ 0 , … , Θ N + 1 } w 0 : N + 1 , t w_{0:N+1,t} w 0 : N + 1 , t
Note: The set marked with * can be omitted, resulting in constant memory complexity.
Calculating the total variation between two univariate distributions, such as in the change-in-mean case with known variance, is straightforward. In this case, θ \theta θ
Proposition 1. The total variation between two univariate Gaussian distributions with means μ i \mu_i μ i μ i + 1 \mu_{i+1} μ i + 1 σ 2 Σ i \sigma^2\Sigma_i σ 2 Σ i σ 2 Σ i + 1 \sigma^2\Sigma_{i+1} σ 2 Σ i + 1
[ Φ ( c 1 − μ i + 1 σ 2 Σ i + 1 ) − Φ ( c 2 − μ i + 1 σ 2 Σ i + 1 ) + Φ ( c 2 − μ i σ 2 Σ i ) − Φ ( c 1 − μ i σ 2 Σ i ) ] . \left[\Phi\left(\frac{c_1-\mu_{i+1}}{\sqrt{\sigma^2\Sigma_{i+1}}}\right)-\Phi\left(\frac{c_2-\mu_{i+1}}{\sqrt{\sigma^2\Sigma_{i+1}}}\right)+\Phi\left(\frac{c_2-\mu_i}{\sqrt{\sigma^2\Sigma_i}}\right)-\Phi\left(\frac{c_1-\mu_i}{\sqrt{\sigma^2\Sigma_i}}\right)\right]. [ Φ ( σ 2 Σ i + 1 c 1 − μ i + 1 ) − Φ ( σ 2 Σ i + 1 c 2 − μ i + 1 ) + Φ ( σ 2 Σ i c 2 − μ i ) − Φ ( σ 2 Σ i c 1 − μ i ) ] . where ( c 1 , c 2 ) = μ i σ 2 Σ i + 1 − μ i + 1 σ 2 Σ i ∓ σ 2 Σ i σ 2 Σ i + 1 ( μ i − μ i + 1 ) 2 + ( σ 2 Σ i + 1 − σ 2 Σ i ) log σ 2 Σ i + 1 σ 2 Σ i σ 2 Σ i + 1 − σ 2 Σ i (c_1,c_2)=\frac{\mu_i\sigma^2\Sigma_{i+1}-\mu_{i+1}\sigma^2\Sigma_i\mp\sqrt{\sigma^2\Sigma_i\sigma^2\Sigma_{i+1}}\sqrt{(\mu_i-\mu_{i+1})^2+(\sigma^2\Sigma_{i+1}-\sigma^2\Sigma_i)\log\frac{\sigma^2\Sigma_{i+1}}{\sigma^2\Sigma_i}}}{\sigma^2\Sigma_{i+1}-\sigma^2\Sigma_i} ( c 1 , c 2 ) = σ 2 Σ i + 1 − σ 2 Σ i μ i σ 2 Σ i + 1 − μ i + 1 σ 2 Σ i ∓ σ 2 Σ i σ 2 Σ i + 1 ( μ i − μ i + 1 ) 2 + ( σ 2 Σ i + 1 − σ 2 Σ i ) l o g σ 2 Σ i σ 2 Σ i + 1 Φ ( c ) \Phi(c) Φ ( c )
Proof. See Appendix § 9 .
Multi-dimensional parameters
In higher dimensions, it is often intractable to find the exact closed form for the total variation between two distributions. Instead, we use the following tractable bound as an approximation to the total variation distance.
For any pair of distributions for distributions θ \theta θ R d \mathbb{R}^d R d f i ( ⋅ ) f_i(\cdot) f i ( ⋅ ) f i + 1 ( ⋅ ) f_{i+1}(\cdot) f i + 1 ( ⋅ ) [19 ] states that
T V ( f i ( θ ) , f i + 1 ( θ ) ) ≤ 1 2 KL ( f i ( θ ) ∥ f i + 1 ( θ ) ) = 1 2 ∫ f i ( θ ) log f i ( θ ) f i + 1 ( θ ) d θ . TV(f_i(\theta),f_{i+1}(\theta))\leq\frac{1}{\sqrt{2}}\sqrt{\operatorname{KL}(f_i(\theta)\mathbin{\|}f_{i+1}(\theta))}=\frac{1}{\sqrt{2}}\sqrt{\int f_i(\theta)\log\frac{f_i(\theta)}{f_{i+1}(\theta)}\,\mathrm{d}\theta}. T V ( f i ( θ ) , f i + 1 ( θ )) ≤ 2 1 KL ( f i ( θ ) ∥ f i + 1 ( θ )) = 2 1 ∫ f i ( θ ) log f i + 1 ( θ ) f i ( θ ) d θ . where KL denotes the Kullback–Leibler divergence. The KL divergence for our posteriors, which will be either multivariate Gaussian or Normal Inverse Gamma, is given by the following results.
Lemma 2.5. (Williams and Rasmussen, [31 ] ) The KL between two multivariate Gaussian distributions with means μ i \mu_i μ i μ i + 1 \mu_{i+1} μ i + 1 σ 2 Σ i \sigma^2\Sigma_i σ 2 Σ i σ 2 Σ i + 1 \sigma^2\Sigma_{i+1} σ 2 Σ i + 1
D KL ( N ( μ i , σ 2 Σ i ) ∥ N ( μ i + 1 , σ 2 Σ i + 1 ) ) = 1 2 ( tr ( Σ i + 1 − 1 Σ i − I ) + 1 σ 2 ( μ i + 1 − μ i ) T Σ i + 1 − 1 ( μ i + 1 − μ i ) + log det ( Σ i − 1 Σ i + 1 ) ) D_{\operatorname{KL}}\left(\mathcal{N}(\mu_i,\sigma^2\Sigma_i)\mathbin{\|}\mathcal{N}(\mu_{i+1},\sigma^2\Sigma_{i+1})\right)=\frac{1}{2}\left(\operatorname{tr}(\Sigma_{i+1}^{-1}\Sigma_i-I)+\frac{1}{\sigma^2}(\mu_{i+1}-\mu_i)^{\mathsf{T}}\Sigma_{i+1}^{-1}(\mu_{i+1}-\mu_i)+\log\det(\Sigma_i^{-1}\Sigma_{i+1})\right) D KL ( N ( μ i , σ 2 Σ i ) ∥ N ( μ i + 1 , σ 2 Σ i + 1 ) ) = 2 1 ( tr ( Σ i + 1 − 1 Σ i − I ) + σ 2 1 ( μ i + 1 − μ i ) T Σ i + 1 − 1 ( μ i + 1 − μ i ) + log det ( Σ i − 1 Σ i + 1 ) ) Lemma 2.6. [25 ] The KL distance between two normal inverse gamma distributions with means μ i \mu_i μ i μ i + 1 \mu_{i+1} μ i + 1 Σ i \Sigma_i Σ i Σ i + 1 \Sigma_{i+1} Σ i + 1 ν i \nu_i ν i ν i + 1 \nu_{i+1} ν i + 1 ι i \iota_i ι i ι i + 1 \iota_{i+1} ι i + 1
D K L [ N I G ( μ i , Σ i , ν i , ι i ) ∥ N I G ( μ i + 1 , Σ i + 1 , ν i + 1 , ι i + 1 ) ] = 1 2 [ ( μ i + 1 − μ i ) T Σ i + 1 − 1 ( μ i + 1 − μ i ) ν i ι i + tr ( Σ i + 1 − 1 Σ i − I ) − log det ( Σ i ) det ( Σ i + 1 ) D_{\mathrm{KL}}[\mathrm{NIG}(\mu_i,\Sigma_i,\nu_i,\iota_i)\Vert\mathrm{NIG}(\mu_{i+1},\Sigma_{i+1},\nu_{i+1},\iota_{i+1})] =
\frac{1}{2}\left[
(\mu_{i+1}-\mu_i)^{\mathsf{T}}\Sigma_{i+1}^{-1}(\mu_{i+1}-\mu_i)\frac{\nu_i}{\iota_i}
+\operatorname{tr}(\Sigma_{i+1}^{-1}\Sigma_i-I)
-\log\frac{\det(\Sigma_i)}{\det(\Sigma_{i+1})}
\right. D KL [ NIG ( μ i , Σ i , ν i , ι i ) ∥ NIG ( μ i + 1 , Σ i + 1 , ν i + 1 , ι i + 1 )] = 2 1 [ ( μ i + 1 − μ i ) T Σ i + 1 − 1 ( μ i + 1 − μ i ) ι i ν i + tr ( Σ i + 1 − 1 Σ i − I ) − log det ( Σ i + 1 ) det ( Σ i ) + ν i + 1 log ι i ι i + 1 − log Γ ( ν i ) Γ ( ν i + 1 ) + ( ν i − ν i + 1 ) ψ ( ν i ) − ( ι i − ι i + 1 ) ν i ι i ] . \left.
+\nu_{i+1}\log\frac{\iota_i}{\iota_{i+1}}
-\log\frac{\Gamma(\nu_i)}{\Gamma(\nu_{i+1})}
+(\nu_i-\nu_{i+1})\psi(\nu_i)
-(\iota_i-\iota_{i+1})\frac{\nu_i}{\iota_i}
\right]. + ν i + 1 log ι i + 1 ι i − log Γ ( ν i + 1 ) Γ ( ν i ) + ( ν i − ν i + 1 ) ψ ( ν i ) − ( ι i − ι i + 1 ) ι i ν i ] . where Γ ( ⋅ ) \Gamma(\cdot) Γ ( ⋅ ) ψ ( ⋅ ) \psi(\cdot) ψ ( ⋅ )
Our algorithm, Bagel, implements merging with this criteria, based on the TV-distance or its approximation, to choose which pair of changepoint sets to merge. To understand the separate benefits of merging rather than pruning from the choice of criteria to merge, in the following we will also compare Bagel to an algorithm that merges the set of changepoints with lowest weight. We call this latter method LP-M. Thus, comparing LP-M with LP-P indicates any advantage of merging over pruning and comparing Bagel with LP-M shows any advantage of our new merging criteria.
We now revisit the example in Figure 2 . Figure 3 presents the cumulative density function (CDF) of the exact approach, Bagel, LP-P and LP-M. It shows that merging preserves a better posterior approximation than pruning, and the total-variation-based merge criterion improves approximation accuracy.
Figure 3. The CDFs of the exact approach (grey), Bagel with recovered posterior distribution (blue), LP-M (orange) and LP-P (black) when M = 50 M = 50 M = 50
Simulation Results
We conduct simulations to compare the posterior approximation accuracy of Bagel, LP-P and LP-M to the exact posterior, and examine whether Bagel accurately approximates the exact Bayesian approach at reduced computational cost.
Accuracy of the posterior approximation
To evaluate posterior approximation accuracy, we first compare the posteriors of changepoint location from Bagel, LP-M, and LP-P to the exact approach on identical streams under identical priors. For each example, we generate T = 1000 T = 1000 T = 1000 τ ∼ U { 200 , … , 800 } \tau\sim U\{200, \ldots, 800\} τ ∼ U { 200 , … , 800 } M = 50 , 100 , 200 M = 50, 100, 200 M = 50 , 100 , 200 § 10 .
To assess the similarity between the exact approach and the approximate approaches, we measure the average total variation distance over the posterior distribution
1 2 R ∑ r = 1 R [ ∣ Pr m r ( τ ≥ T ∣ y 1 : T ) − Pr ∞ r ( τ ≥ T ∣ y 1 : T ) ∣ + ∑ i = 1 T − 1 ∣ Pr m r ( τ = i ∣ y 1 : T ) − Pr ∞ r ( τ = i ∣ y 1 : T ) ∣ ] \frac{1}{2R}\sum_{r=1}^{R}\left[\left|\Pr_{m}^{r}(\tau\ge T\mid y_{1:T})-\Pr_{\infty}^{r}(\tau\ge T\mid y_{1:T})\right|+\sum_{i=1}^{T-1}\left|\Pr_{m}^{r}(\tau=i\mid y_{1:T})-\Pr_{\infty}^{r}(\tau=i\mid y_{1:T})\right|\right] 2 R 1 r = 1 ∑ R [ m Pr r ( τ ≥ T ∣ y 1 : T ) − ∞ Pr r ( τ ≥ T ∣ y 1 : T ) + i = 1 ∑ T − 1 m Pr r ( τ = i ∣ y 1 : T ) − ∞ Pr r ( τ = i ∣ y 1 : T ) ] and the average Wasserstein distance on the conditional posterior of changepoint location
1 R ∑ r = 1 R ∑ i = 1 T − 1 ∣ Pr m r ( τ ≤ i ∣ τ < T , y 1 : T ) − Pr ∞ r ( τ ≤ i ∣ τ < T , y 1 : T ) ∣ \frac{1}{R}\sum_{r=1}^{R}\sum_{i=1}^{T-1}\left|\Pr_{m}^{r}(\tau\le i\mid\tau<T,y_{1:T})-\Pr_{\infty}^{r}(\tau\le i\mid\tau<T,y_{1:T})\right| R 1 r = 1 ∑ R i = 1 ∑ T − 1 m Pr r ( τ ≤ i ∣ τ < T , y 1 : T ) − ∞ Pr r ( τ ≤ i ∣ τ < T , y 1 : T ) over R R R Pr m r = Pr ( τ ∣ y 1 : T , M = m ) \Pr_{m}^{r}=\Pr(\tau\mid y_{1:T},M=m) Pr m r = Pr ( τ ∣ y 1 : T , M = m ) m m m Pr ∞ r = Pr ( τ ∣ y 1 : T , M = ∞ ) \Pr_{\infty}^{r}=\Pr(\tau\mid y_{1:T},M=\infty) Pr ∞ r = Pr ( τ ∣ y 1 : T , M = ∞ ) r r r
Simulation results are shown in Table 1 . The approximation accuracy improves as the budget M M M
Table 1. Simulation on a fixed horizon T = 1000 T = 1000 T = 1000 W 1 , c h a n g e W_{1,\mathrm{change}} W 1 , change T V f u l l \mathrm{TV}_{\mathrm{full}} TV full
Accuracy under online setting
We now compare the performance of methods in terms of whether and how quickly they detect a change, and, once they have detected one, how accurate inference about where that change is. We simulate data with a fixed size of change and changepoint location of 1000. We set the pre-change prior mean to match the data-generating process. As before, we consider a main prior and a diffuse one. In each case, the prior for the change size γ \gamma γ § 11 . Different from the previous study, each method stops at the first time its detection statistic exceeds a predefined threshold. The threshold is calibrated on data of length 2000 without a change and set to achieve 5% false-alarm rate over 500 replications.
We look at a range of properties of each method, including the accuracy of the posterior distribution of the location of a change when a change is detected. To do this, for each approximate method, we rerun the exact method on the same data at the stopping time of the approximate method, to obtain the exact posterior. For replicate r r r T a r T_a^r T a r τ ^ a r \hat{\tau}_a^r τ ^ a r H a r H_a^r H a r τ ^ e r ( T a r ) \hat{\tau}_e^r(T_a^r) τ ^ e r ( T a r ) H e r ( T a r ) H_e^r(T_a^r) H e r ( T a r ) T a r T_a^r T a r Table 2 shows. The invalid percentage shows the proportion of replicates that trigger false alarms or fail to alarm. The average detection delay provides a reference detection speed. The MAP distance measures the distance between the approximate and exact MAP estimates. The Jaccard index and HPD mass error measure the accuracy of HPD sets at the nominated level.
Table 2. Evaluation metrics used in Table 3 . V = { r : τ < T a r ≤ 2000 } V=\{r:\tau<T_a^r\le2000\} V = { r : τ < T a r ≤ 2000 } P e P_e P e R R R
Table 3 presents the simulation results for two representative examples under the main prior. Additional results under both main and diffuse priors are provided in Appendix § 11 . Bagel and LP-M have fewer invalid replicates and shorter delays compared to the pruned approach. When evaluating the posterior approximation accuracy, Bagel consistently provides a good approximation, with the lowest HPD error and the highest Jaccard. This indicates that the HPD set inferred by Bagel accurately preserves the changepoint candidates and their corresponding probabilities. On the other hand, LP-P can preserve the accuracy of the MAP estimator in most cases but has poor estimated HPD sets. The lower Jaccard overlap and larger HPD error indicate that LP-P can lead to misleading uncertainty quantification. The performance of LP-M sits between LP-P and Bagel in most cases. It can provide a better approximation than LP-P but has lower Jaccard and larger HPD error than Bagel, and can introduce small shifts in the MAP estimator. In the continuous slope change, Bagel and LP-M have similar performance. This is because neighbouring candidates may already have highly similar post-change distributions due to continuity constraints.
Table 3. Simulation result on a fixed horizon T = 2000 T = 2000 T = 2000 τ = 1000 \tau= 1000 τ = 1000
Empirical runtime
We evaluate the empirical runtime of the exact approach and Bagel. The main code of Bagel is implemented in C++, and invoked per step from the R interface. Figure 4 shows the runtime when varying M M M M M M § 12 . Bagel is substantially faster and more stable than the exact approach over long monitoring horizons. These results show that merging preserves the accuracy of the exact method while remaining fast to implement, making Bagel computationally practical for high-frequency, long-run and real-world monitoring tasks.
Figure 4. Median per-step runtime of Bagel on T = 1000 T = 1000 T = 1000 M = 50 , 100 , 200 M = 50, 100, 200 M = 50 , 100 , 200 t = M t = M t = M
Real data example - Azure LLM inference traces
Large language model (LLM) inference workloads are dynamic over time. Changes in generated-token count reflect the change decoding demand, which further affects latency. Detecting these changes is useful for workload monitoring, diagnosis and capacity planning. We analyse the generated token count for code traffic from the Azure LLM inference trace ([28 ] ) available from https://github.com/Azure/AzurePublicDataset/ . The data records the request timestamps and generated-token counts over one week in 2024 from Friday to Thursday. We aggregate the requests into five-minute intervals (2016 bins), and monitor the mean generated-token count for code traffic within each interval as a summary of decoding demand. We assume the first day (288 bins) has no change to calibrate the alarming threshold, and the remaining six days are for monitoring. The model is
y t = β 0 + β 1 t − 1 288 + β sin sin 2 π h t 288 + β cos cos 2 π h t 288 + 1 ( t > τ ) [ γ 0 + γ 1 t − τ 288 + γ sin sin 2 π h t 288 + γ cos cos 2 π h t 288 ] + σ ϵ t , \begin{aligned}
y_t &= \beta_0 + \beta_1 \frac{t-1}{288} + \beta_{\sin} \sin\frac{2\pi h t}{288} + \beta_{\cos} \cos\frac{2\pi h t}{288} \\
&\quad+ \mathbf{1}(t > \tau)\left[\gamma_0 + \gamma_1 \frac{t-\tau}{288} + \gamma_{\sin} \sin\frac{2\pi h t}{288} + \gamma_{\cos} \cos\frac{2\pi h t}{288}\right] + \sigma\epsilon_t,
\end{aligned} y t = β 0 + β 1 288 t − 1 + β s i n sin 288 2 π h t + β c o s cos 288 2 π h t + 1 ( t > τ ) [ γ 0 + γ 1 288 t − τ + γ s i n sin 288 2 π h t + γ c o s cos 288 2 π h t ] + σ ϵ t , where σ 2 ∼ IG ( ν , ι ) \sigma^2 \sim\operatorname{IG}(\nu,\iota) σ 2 ∼ IG ( ν , ι ) β ∣ σ 2 ∼ N ( 0 , 100 σ 2 I 4 ) \beta\mid\sigma^2 \sim N(0,100\sigma^2 I_4) β ∣ σ 2 ∼ N ( 0 , 100 σ 2 I 4 ) σ 2 ∼ IG ( 2 , 1 ) \sigma^2 \sim\operatorname{IG}(2,1) σ 2 ∼ IG ( 2 , 1 )
Figure 5 presents the detection results. The MAP estimator of changepoint locations is at bins 400, 675 and 981, with alarms at indices 407, 953 and 1001. The estimated duration between MAP and alarm is 7, 278 and 20 observations. The estimated changepoint locations divide the monitored observations into four segments: early Saturday morning, Saturday morning to Sunday morning, Sunday morning to Monday morning, and Monday through Thursday. The second change is associated with changes in the daily seasonality, so it requires longer periods to accumulate sufficient evidence. The first and third alarms are easier to detect, as both of them have level change as well as seasonality change. Although there is no ground truth label for anomalies, the detected changepoint matches visible shifts in the level and daily pattern.
Figure 5. The detection results of Bagel on Azure LLM inference trace data with M = 100 M = 100 M = 100 − log w 0 -\log w_0 − log w 0 w 0 w_0 w 0 − log w 0 > c -\log w_0 > c − log w 0 > c 1 − w 0 > 1 − e − c 1-w_0 > 1-e^{-c} 1 − w 0 > 1 − e − c
Discussion
We propose Bagel, a computationally efficient Bayesian online changepoint detection algorithm for linear models that has the ability to detect changes in mean and slope and can incorporate seasonality. By merging similar neighbouring changepoint candidates, Bagel limits the number of posterior particles that must be updated so the computational complexity is reduced to a constant per observation after reaching the limit. The merged posterior can be recovered later to provide uncertainty quantification. Through the simulation, we show that Bagel maintains a better approximation of the exact posterior compared to LP-P and LP-M. Therefore, under strict latency requirements and limited computational budgets, Bagel can provide stable and fast sequential updates while retaining rich posterior information. This is particularly useful in practice, as information about when the change happens and its posterior probability can help downstream root-cause analysis and decision-making, and its efficient computation also makes it practical to deploy without demanding substantial resources. The Bagel package is available at https://github.com/grosed/bagel.git .
Acknowledgements This work is supported by EPSRC grants EP/N031938/1, EP/Z531327/1, EP/Y028783/1, EP/S022252/1 and BT Group.
AI Declaration: During the preparation of this work, Claude Code Opus 4.8 and Codex 5.6-sol were used to provide code assistance for the simulation and real data experiment. Specifically, they were used to assist with code refactoring, drafting, debugging, and formatting the R output to LaTeX for analyses specified by the authors. After using this tool, the first author reviewed and tested the code as needed and takes full responsibility for the final outputs.
Properties of the Prior in Section 2.1.
Example 1
To see that the second prior for Example 1 gives the same distribution for the pre-change and post-change mean, and that these are independent, we apply a change of variable. We have ( β , γ ) T (\beta,\gamma)^{\mathsf{T}} ( β , γ ) T
[ μ 1 0 ] \begin{bmatrix} \mu_1 \\ 0 \end{bmatrix} [ μ 1 0 ] and variance σ 2 \sigma^2 σ 2
[ δ 1 − δ 1 − δ 1 2 δ 1 ] . \begin{bmatrix} \delta_1 & -\delta_1 \\ -\delta_1 & 2\delta_1 \end{bmatrix}. [ δ 1 − δ 1 − δ 1 2 δ 1 ] . Let the pre-change and post-change means be ( β , θ ) T (\beta,\theta)^{\mathsf{T}} ( β , θ ) T
[ β θ ] = [ 1 0 1 1 ] [ β γ ] . \begin{bmatrix} \beta\\ \theta\end{bmatrix}
=
\begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}
\begin{bmatrix} \beta\\ \gamma\end{bmatrix}. [ β θ ] = [ 1 1 0 1 ] [ β γ ] . This is a linear transformation, so ( β , θ ) T (\beta,\theta)^{\mathsf{T}} ( β , θ ) T
[ 1 0 1 1 ] [ μ 1 0 ] = [ μ 1 μ 1 ] , \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}
\begin{bmatrix} \mu_1 \\ 0 \end{bmatrix}
=
\begin{bmatrix} \mu_1 \\ \mu_1 \end{bmatrix}, [ 1 1 0 1 ] [ μ 1 0 ] = [ μ 1 μ 1 ] , and variance
σ 2 [ 1 0 1 1 ] [ δ 1 − δ 1 − δ 1 2 δ 1 ] [ 1 0 1 1 ] T = σ 2 [ δ 1 0 0 δ 1 ] , \sigma^2
\begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}
\begin{bmatrix} \delta_1 & -\delta_1 \\ -\delta_1 & 2\delta_1 \end{bmatrix}
\begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}^{\mathsf{T}}
=
\sigma^2
\begin{bmatrix} \delta_1 & 0 \\ 0 & \delta_1 \end{bmatrix}, σ 2 [ 1 1 0 1 ] [ δ 1 − δ 1 − δ 1 2 δ 1 ] [ 1 1 0 1 ] T = σ 2 [ δ 1 0 0 δ 1 ] , as required.
Example 2
In the change-in-slope case, first consider the continuous change model. We can re-parameterise the model to the pre-change intercept and slope, the change in the mean at the changepoint and the change in slope. Denote these parameters as ( β 1 , β 2 , ψ , γ 2 ) T (\beta_1,\beta_2,\psi,\gamma_2)^{\mathsf{T}} ( β 1 , β 2 , ψ , γ 2 ) T
[ β 1 β 2 ψ γ 2 ] = [ 1 0 0 0 0 1 0 0 0 0 1 τ 0 0 0 1 ] [ β 1 β 2 γ 1 γ 2 ] . \begin{bmatrix} \beta_1 \\ \beta_2 \\ \psi\\ \gamma_2 \end{bmatrix}
=
\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & \tau\\ 0 & 0 & 0 & 1 \end{bmatrix}
\begin{bmatrix} \beta_1 \\ \beta_2 \\ \gamma_1 \\ \gamma_2 \end{bmatrix}. β 1 β 2 ψ γ 2 = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 τ 1 β 1 β 2 γ 1 γ 2 . We can calculate the mean and variance of ( β 1 , β 2 , ψ , γ 2 ) T (\beta_1,\beta_2,\psi,\gamma_2)^{\mathsf{T}} ( β 1 , β 2 , ψ , γ 2 ) T
[ 1 0 0 0 0 1 0 0 0 0 1 τ 0 0 0 1 ] [ μ 1 μ 2 0 0 ] = [ μ 1 μ 2 0 0 ] , \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & \tau\\ 0 & 0 & 0 & 1 \end{bmatrix}
\begin{bmatrix} \mu_1 \\ \mu_2 \\ 0 \\ 0 \end{bmatrix}
=
\begin{bmatrix} \mu_1 \\ \mu_2 \\ 0 \\ 0 \end{bmatrix}, 1 0 0 0 0 1 0 0 0 0 1 0 0 0 τ 1 μ 1 μ 2 0 0 = μ 1 μ 2 0 0 , and
σ 2 ( 1 0 0 0 0 1 0 0 0 0 1 τ 0 0 0 1 ) ( δ 1 0 0 0 0 δ 2 0 0 0 0 τ 2 δ 3 − τ δ 3 0 0 − τ δ 3 δ 3 ) ( 1 0 0 0 0 1 0 0 0 0 1 τ 0 0 0 1 ) T = σ 2 ( δ 1 0 0 0 0 δ 2 0 0 0 0 0 0 0 0 0 δ 3 ) . \sigma^2
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & \tau\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\delta_1 & 0 & 0 & 0 \\
0 & \delta_2 & 0 & 0 \\
0 & 0 & \tau^2\delta_3 & -\tau\delta_3 \\
0 & 0 & -\tau\delta_3 & \delta_3
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & \tau\\
0 & 0 & 0 & 1
\end{pmatrix}^{\mathsf{T}}
=
\sigma^2
\begin{pmatrix}
\delta_1 & 0 & 0 & 0 \\
0 & \delta_2 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & \delta_3
\end{pmatrix}. σ 2 1 0 0 0 0 1 0 0 0 0 1 0 0 0 τ 1 δ 1 0 0 0 0 δ 2 0 0 0 0 τ 2 δ 3 − τ δ 3 0 0 − τ δ 3 δ 3 1 0 0 0 0 1 0 0 0 0 1 0 0 0 τ 1 T = σ 2 δ 1 0 0 0 0 δ 2 0 0 0 0 0 0 0 0 0 δ 3 . This shows that, ψ \psi ψ
Now consider a discontinuous change. This time we re-parameterise in terms of the pre-change and post-change intercept and slope, which we denote ( β 1 , β 2 , θ 1 , θ 2 ) T (\beta_1,\beta_2,\theta_1,\theta_2)^{\mathsf{T}} ( β 1 , β 2 , θ 1 , θ 2 ) T
( β 1 β 2 θ 1 θ 2 ) = ( 1 0 0 0 0 1 0 0 1 τ 1 τ 0 1 0 1 ) ( β 1 β 2 γ 1 γ 2 ) . \begin{pmatrix}
\beta_1 \\
\beta_2 \\
\theta_1 \\
\theta_2
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & \tau& 1 & \tau\\
0 & 1 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\beta_1 \\
\beta_2 \\
\gamma_1 \\
\gamma_2
\end{pmatrix}. β 1 β 2 θ 1 θ 2 = 1 0 1 0 0 1 τ 1 0 0 1 0 0 0 τ 1 β 1 β 2 γ 1 γ 2 . With this linear transformation, and using the prior specified in Section 2.1 for the discontinuous model, we have that the mean of ( β 1 , β 2 , θ 1 , θ 2 ) T (\beta_1,\beta_2,\theta_1,\theta_2)^{\mathsf{T}} ( β 1 , β 2 , θ 1 , θ 2 ) T
( 1 0 0 0 0 1 0 0 1 τ 1 τ 0 1 0 1 ) ( μ 1 μ 2 − μ 2 τ 0 ) = ( μ 1 μ 2 μ 1 μ 2 ) , \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & \tau& 1 & \tau\\
0 & 1 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\mu_1 \\
\mu_2 \\
-\mu_2\tau\\
0
\end{pmatrix}
=
\begin{pmatrix}
\mu_1 \\
\mu_2 \\
\mu_1 \\
\mu_2
\end{pmatrix}, 1 0 1 0 0 1 τ 1 0 0 1 0 0 0 τ 1 μ 1 μ 2 − μ 2 τ 0 = μ 1 μ 2 μ 1 μ 2 , and its variance is
σ 2 ( 1 0 0 0 0 1 0 0 1 τ 1 τ 0 1 0 1 ) ( δ 1 0 − δ 1 0 0 δ 2 0 − δ 2 − δ 1 0 2 δ 1 + τ 2 δ 2 − τ δ 2 0 − δ 2 − τ δ 2 2 δ 2 ) ( 1 0 0 0 0 1 0 0 1 τ 1 τ 0 1 0 1 ) T = σ 2 ( δ 1 0 0 0 0 δ 2 0 0 0 0 δ 1 0 0 0 0 δ 2 ) . \sigma^2
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & \tau& 1 & \tau\\
0 & 1 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\delta_1 & 0 & -\delta_1 & 0 \\
0 & \delta_2 & 0 & -\delta_2 \\
-\delta_1 & 0 & 2\delta_1+\tau^2\delta_2 & -\tau\delta_2 \\
0 & -\delta_2 & -\tau\delta_2 & 2\delta_2
\end{pmatrix}
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & \tau& 1 & \tau\\
0 & 1 & 0 & 1
\end{pmatrix}^{\mathsf{T}}
=
\sigma^2
\begin{pmatrix}
\delta_1 & 0 & 0 & 0 \\
0 & \delta_2 & 0 & 0 \\
0 & 0 & \delta_1 & 0 \\
0 & 0 & 0 & \delta_2
\end{pmatrix}. σ 2 1 0 1 0 0 1 τ 1 0 0 1 0 0 0 τ 1 δ 1 0 − δ 1 0 0 δ 2 0 − δ 2 − δ 1 0 2 δ 1 + τ 2 δ 2 − τ δ 2 0 − δ 2 − τ δ 2 2 δ 2 1 0 1 0 0 1 τ 1 0 0 1 0 0 0 τ 1 T = σ 2 δ 1 0 0 0 0 δ 2 0 0 0 0 δ 1 0 0 0 0 δ 2 . So the pre-change and post-change parameters are independent and follow the same distribution, as required.
Example 3. We consider the prior for Example 3 that the pre-change and post-change distributions are independent and identically distributed. Our prior is that ( β sin , β cos , γ sin , γ cos ) T (\beta_{\sin},\beta_{\cos},\gamma_{\sin},\gamma_{\cos})^{\mathsf{T}} ( β s i n , β c o s , γ s i n , γ c o s ) T [ μ sin μ cos 0 0 ] T [\mu_{\sin}\ \ \mu_{\cos}\ \ 0\ \ 0]^{\mathsf{T}} [ μ s i n μ c o s 0 0 ] T
σ 2 ( δ sin 0 − δ sin 0 0 δ cos 0 − δ cos − δ sin 0 2 δ sin 0 0 − δ cos 0 2 δ cos ) . \sigma^2
\begin{pmatrix}
\delta_{\sin} & 0 & -\delta_{\sin} & 0 \\
0 & \delta_{\cos} & 0 & -\delta_{\cos} \\
-\delta_{\sin} & 0 & 2\delta_{\sin} & 0 \\
0 & -\delta_{\cos} & 0 & 2\delta_{\cos}
\end{pmatrix}. σ 2 δ s i n 0 − δ s i n 0 0 δ c o s 0 − δ c o s − δ s i n 0 2 δ s i n 0 0 − δ c o s 0 2 δ c o s . Let the pre-change and post-change coefficients of the sin and cos terms be ( β sin , β cos , θ sin , θ cos ) T (\beta_{\sin},\beta_{\cos},\theta_{\sin},\theta_{\cos})^{\mathsf{T}} ( β s i n , β c o s , θ s i n , θ c o s ) T
[ β sin β cos θ sin θ cos ] = [ 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 1 ] [ β sin β cos γ sin γ cos ] . \begin{bmatrix}
\beta_{\sin}\\
\beta_{\cos}\\
\theta_{\sin}\\
\theta_{\cos}
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
1 & 0 & 1 & 0\\
0 & 1 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\beta_{\sin}\\
\beta_{\cos}\\
\gamma_{\sin}\\
\gamma_{\cos}
\end{bmatrix}. β s i n β c o s θ s i n θ c o s = 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 β s i n β c o s γ s i n γ c o s . So ( β sin , β cos , θ sin , θ cos ) T (\beta_{\sin},\beta_{\cos},\theta_{\sin},\theta_{\cos})^{\mathsf{T}} ( β s i n , β c o s , θ s i n , θ c o s ) T
[ 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 1 ] [ μ sin μ cos 0 0 ] = [ μ sin μ cos μ sin μ cos ] , \begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
1 & 0 & 1 & 0\\
0 & 1 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\mu_{\sin}\\
\mu_{\cos}\\
0\\
0
\end{bmatrix}
=
\begin{bmatrix}
\mu_{\sin}\\
\mu_{\cos}\\
\mu_{\sin}\\
\mu_{\cos}
\end{bmatrix}, 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 μ s i n μ c o s 0 0 = μ s i n μ c o s μ s i n μ c o s , and variance
σ 2 [ 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 1 ] [ δ sin 0 − δ sin 0 0 δ cos 0 − δ cos − δ sin 0 2 δ sin 0 0 − δ cos 0 2 δ cos ] [ 1 0 0 0 0 1 0 0 1 0 1 0 0 1 0 1 ] T = σ 2 [ δ sin 0 0 0 0 δ cos 0 0 0 0 δ sin 0 0 0 0 δ cos ] , \sigma^2
\begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
1 & 0 & 1 & 0\\
0 & 1 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\delta_{\sin} & 0 & -\delta_{\sin} & 0\\
0 & \delta_{\cos} & 0 & -\delta_{\cos}\\
-\delta_{\sin} & 0 & 2\delta_{\sin} & 0\\
0 & -\delta_{\cos} & 0 & 2\delta_{\cos}
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
1 & 0 & 1 & 0\\
0 & 1 & 0 & 1
\end{bmatrix}^{\mathsf{T}}
=
\sigma^2
\begin{bmatrix}
\delta_{\sin} & 0 & 0 & 0\\
0 & \delta_{\cos} & 0 & 0\\
0 & 0 & \delta_{\sin} & 0\\
0 & 0 & 0 & \delta_{\cos}
\end{bmatrix}, σ 2 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 δ s i n 0 − δ s i n 0 0 δ c o s 0 − δ c o s − δ s i n 0 2 δ s i n 0 0 − δ c o s 0 2 δ c o s 1 0 1 0 0 1 0 1 0 0 1 0 0 0 0 1 T = σ 2 δ s i n 0 0 0 0 δ c o s 0 0 0 0 δ s i n 0 0 0 0 δ c o s , as required.
Proofs for sequential updating
Proof of Theorem 2.1 For i = 1 , … , t − 2 i = 1,\ldots,t - 2 i = 1 , … , t − 2
w i , t = Pr t ( τ = i ) p ( y 1 : t ∣ τ = i ) = Pr t − 1 ( τ = i ) p ( y 1 : t − 1 ∣ τ = i ) Pr t ( τ = i ) Pr t − 1 ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) = w i , t − 1 Pr t ( τ = i ) Pr t − 1 ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) . \begin{aligned}
w_{i,t} &= \Pr_t(\tau= i)p(y_{1:t}\mid\tau= i)\\
&= \Pr_{t-1}(\tau= i)p(y_{1:t-1}\mid\tau= i)\frac{\Pr_t(\tau= i)}{\Pr_{t-1}(\tau= i)}p(y_t\mid y_{1:t-1},\tau= i)\\
&= w_{i,t-1}\frac{\Pr_t(\tau= i)}{\Pr_{t-1}(\tau= i)}p(y_t\mid y_{1:t-1},\tau= i).
\end{aligned} w i , t = t Pr ( τ = i ) p ( y 1 : t ∣ τ = i ) = t − 1 Pr ( τ = i ) p ( y 1 : t − 1 ∣ τ = i ) Pr t − 1 ( τ = i ) Pr t ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) = w i , t − 1 Pr t − 1 ( τ = i ) Pr t ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = i ) . A similar simple argument applies for i = 0 i = 0 i = 0 i = t − 1 i = t - 1 i = t − 1
w t − 1 , t = Pr t ( τ = t − 1 ) p ( y 1 : t ∣ τ = t − 1 ) = Pr t − 1 ( τ ≥ t − 1 ) p ( y 1 : t − 1 ∣ τ ≥ t − 1 ) Pr t ( τ = i ) Pr t − 1 ( τ ≥ t − 1 ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) = w 0 , t − 1 Pr t ( τ = i ) Pr t − 1 ( τ ≥ t − 1 ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) . \begin{aligned}
w_{t-1,t} &= \Pr_t(\tau= t - 1)p(y_{1:t}\mid\tau= t - 1)\\
&= \Pr_{t-1}(\tau\ge t - 1)p(y_{1:t-1}\mid\tau\ge t - 1)\frac{\Pr_t(\tau= i)}{\Pr_{t-1}(\tau\ge t - 1)}p(y_t\mid y_{1:t-1},\tau= t - 1)\\
&= w_{0,t-1}\frac{\Pr_t(\tau= i)}{\Pr_{t-1}(\tau\ge t - 1)}p(y_t\mid y_{1:t-1},\tau= t - 1).
\end{aligned} w t − 1 , t = t Pr ( τ = t − 1 ) p ( y 1 : t ∣ τ = t − 1 ) = t − 1 Pr ( τ ≥ t − 1 ) p ( y 1 : t − 1 ∣ τ ≥ t − 1 ) Pr t − 1 ( τ ≥ t − 1 ) Pr t ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) = w 0 , t − 1 Pr t − 1 ( τ ≥ t − 1 ) Pr t ( τ = i ) p ( y t ∣ y 1 : t − 1 , τ = t − 1 ) . Proof of Theorem 2.2
If we have a change at t − 1 t-1 t − 1 ( β , γ ) (\beta,\gamma) ( β , γ ) μ t − 1 \mu^{t-1} μ t − 1 σ 2 Σ t − 1 \sigma^2\Sigma^{t-1} σ 2 Σ t − 1 β \beta β γ \gamma γ β \beta β t − 1 t-1 t − 1
The posterior distribution given y 1 : t − 1 y_{1:t-1} y 1 : t − 1 t − 1 t-1 t − 1
p ( β , γ , σ ∣ y 1 : t − 1 , τ = t − 1 ) ∝ p ( σ ) p ( β ∣ σ ) p ( γ ∣ β , σ , τ = t − 1 ) p ( y 1 : t − 1 ∣ β , γ , σ , τ = t − 1 ) = p ( σ ) p ( β ∣ σ ) p ( γ ∣ β , σ , τ = t − 1 ) p ( y 1 : t − 1 ∣ β , σ , τ = t − 1 ) ∝ p ( β , σ ∣ y 1 : t − 1 , τ = t − 1 ) p ( γ ∣ β , σ , τ = t − 1 ) = p ( β , σ ∣ y 1 : t − 1 , τ ≥ t − 1 ) p ( γ ∣ β , σ , τ = t − 1 ) . \begin{aligned}
p(\beta,\gamma,\sigma\mid y_{1:t-1},\tau=t-1)
&\propto p(\sigma)p(\beta\mid\sigma)p(\gamma\mid\beta,\sigma,\tau=t-1)p(y_{1:t-1}\mid\beta,\gamma,\sigma,\tau=t-1)\\
&=p(\sigma)p(\beta\mid\sigma)p(\gamma\mid\beta,\sigma,\tau=t-1)p(y_{1:t-1}\mid\beta,\sigma,\tau=t-1)\\
&\propto p(\beta,\sigma\mid y_{1:t-1},\tau=t-1)p(\gamma\mid\beta,\sigma,\tau=t-1)\\
&=p(\beta,\sigma\mid y_{1:t-1},\tau\ge t-1)p(\gamma\mid\beta,\sigma,\tau=t-1).
\end{aligned} p ( β , γ , σ ∣ y 1 : t − 1 , τ = t − 1 ) ∝ p ( σ ) p ( β ∣ σ ) p ( γ ∣ β , σ , τ = t − 1 ) p ( y 1 : t − 1 ∣ β , γ , σ , τ = t − 1 ) = p ( σ ) p ( β ∣ σ ) p ( γ ∣ β , σ , τ = t − 1 ) p ( y 1 : t − 1 ∣ β , σ , τ = t − 1 ) ∝ p ( β , σ ∣ y 1 : t − 1 , τ = t − 1 ) p ( γ ∣ β , σ , τ = t − 1 ) = p ( β , σ ∣ y 1 : t − 1 , τ ≥ t − 1 ) p ( γ ∣ β , σ , τ = t − 1 ) . Here we have used the fact that the observations prior to the change do not depend on γ \gamma γ β \beta β σ \sigma σ y 1 : t − 1 y_{1:t-1} y 1 : t − 1 τ = t − 1 \tau=t-1 τ = t − 1 τ ≥ t − 1 \tau\ge t-1 τ ≥ t − 1
The posterior for ( β , σ ) (\beta,\sigma) ( β , σ ) y 1 : t − 1 y_{1:t-1} y 1 : t − 1 τ ≥ t − 1 \tau\ge t-1 τ ≥ t − 1 t − 1 t-1 t − 1 t − 1 t-1 t − 1 σ \sigma σ β \beta β μ t − 1 , 0 \mu^{t-1,0} μ t − 1 , 0 Σ t − 1 , 0 \Sigma^{t-1,0} Σ t − 1 , 0
As p ( β , γ ∣ σ , τ = t − 1 ) p(\beta,\gamma\mid\sigma,\tau=t-1) p ( β , γ ∣ σ , τ = t − 1 ) p ( γ ∣ β , σ , τ = t − 1 ) p(\gamma\mid\beta,\sigma,\tau=t-1) p ( γ ∣ β , σ , τ = t − 1 ) β \beta β
γ ∣ β , σ , τ = t − 1 ∼ N ( μ γ + Σ γ , β Σ β , β − 1 ( β − μ β ) , σ 2 ( Σ γ , γ − Σ γ , β Σ β , β − 1 Σ β , γ ) ) \gamma\mid\beta,\sigma,\tau=t-1\sim\mathcal{N}\left(\mu_\gamma+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}(\beta-\mu_\beta),\sigma^2\left(\Sigma_{\gamma,\gamma}-\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}\right)\right) γ ∣ β , σ , τ = t − 1 ∼ N ( μ γ + Σ γ , β Σ β , β − 1 ( β − μ β ) , σ 2 ( Σ γ , γ − Σ γ , β Σ β , β − 1 Σ β , γ ) ) Therefore, the marginal mean of γ \gamma γ β \beta β
E ( γ ) = E ( E ( γ ∣ β ) ) = μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) . E(\gamma)=E\left(E(\gamma\mid\beta)\right)=\mu_\gamma+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\left(\mu^{t-1,0}-\mu_\beta\right). E ( γ ) = E ( E ( γ ∣ β ) ) = μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) . The marginal scaled variance is
Var ( γ ) = E ( Var ( γ ∣ β ) ) + Var ( E ( γ ∣ β ) ) = σ 2 ( Σ γ , γ − Σ γ , β Σ β , β − 1 Σ β , γ + Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ ) = σ 2 ( Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ) . \begin{aligned}
\operatorname{Var}(\gamma)&=E\left(\operatorname{Var}(\gamma\mid\beta)\right)+\operatorname{Var}\left(E(\gamma\mid\beta)\right)\\
&=\sigma^2\left(\Sigma_{\gamma,\gamma}-\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma^{t-1,0}\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}\right)\\
&=\sigma^2\left(\Sigma_{\gamma,\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\left(\Sigma^{t-1,0}-\Sigma_{\beta,\beta}\right)\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}\right).
\end{aligned} Var ( γ ) = E ( Var ( γ ∣ β ) ) + Var ( E ( γ ∣ β ) ) = σ 2 ( Σ γ , γ − Σ γ , β Σ β , β − 1 Σ β , γ + Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ ) = σ 2 ( Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ) . The covariance can be calculated as
Cov ( γ , β ) = σ 2 Cov ( Σ γ , β Σ β , β − 1 β , β ) = σ 2 Σ γ , β Σ β , β − 1 Σ t − 1 , 0 . \operatorname{Cov}(\gamma,\beta)=\sigma^2\operatorname{Cov}\left(\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\beta,\beta\right)=\sigma^2\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma^{t-1,0}. Cov ( γ , β ) = σ 2 Cov ( Σ γ , β Σ β , β − 1 β , β ) = σ 2 Σ γ , β Σ β , β − 1 Σ t − 1 , 0 . Thus, the full updated joint distribution of ( β , γ ) (\beta,\gamma) ( β , γ )
N ( [ μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) ] , σ 2 [ Σ t − 1 , 0 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ] ) . N\left(
\begin{bmatrix}
\mu_{\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}(\mu^{t-1,0}-\mu_{\beta})\\
\end{bmatrix},
\sigma^2
\begin{bmatrix}
\Sigma^{t-1,0} & \Sigma^{t-1,0}\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}\\
\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}\Sigma^{t-1,0} & \Sigma_{\gamma,\gamma}+\Sigma_{\gamma,\beta}\Sigma_{\beta,\beta}^{-1}(\Sigma^{t-1,0}-\Sigma_{\beta,\beta})\Sigma_{\beta,\beta}^{-1}\Sigma_{\beta,\gamma}
\end{bmatrix}
\right). N ( [ μ γ + Σ γ , β Σ β , β − 1 ( μ t − 1 , 0 − μ β ) ] , σ 2 [ Σ t − 1 , 0 Σ γ , β Σ β , β − 1 Σ t − 1 , 0 Σ t − 1 , 0 Σ β , β − 1 Σ β , γ Σ γ , γ + Σ γ , β Σ β , β − 1 ( Σ t − 1 , 0 − Σ β , β ) Σ β , β − 1 Σ β , γ ] ) . as required.
Proof of Theorems 2.3 and 2.4
The results in Theorems 2.3 and 2.4 follow from standard results. One way to see this is to view our model as a dynamic linear model (West and Harrison, [30 ] ), where the hidden state is the parameter, but we define the dynamics such that the hidden state does not change over time.
For the case of no change, the dynamic linear model is
β t = β t − 1 y t = h 0 T β t + σ ϵ t , \begin{aligned}
\beta_t &= \beta_{t-1}\\
y_t &= h_0^{T}\beta_t+\sigma\epsilon_t,
\end{aligned} β t y t = β t − 1 = h 0 T β t + σ ϵ t , where ϵ t ∼ N ( 0 , 1 ) \epsilon_t\sim\mathcal{N}(0,1) ϵ t ∼ N ( 0 , 1 ) h 0 h_0 h 0 t − 1 t-1 t − 1 y 1 : t − 1 y_{1:t-1} y 1 : t − 1
β t − 1 ∣ σ 2 , y 1 : t − 1 , τ ≥ t − 1 ∼ N ( μ t − 1 , 0 , σ 2 Σ t − 1 , 0 ) , \beta_{t-1}\mid\sigma^2,y_{1:t-1},\tau\geq t-1\sim\mathcal{N}(\mu^{t-1,0},\sigma^2\Sigma^{t-1,0}), β t − 1 ∣ σ 2 , y 1 : t − 1 , τ ≥ t − 1 ∼ N ( μ t − 1 , 0 , σ 2 Σ t − 1 , 0 ) , σ 2 ∣ y 1 : t − 1 , τ ≥ t − 1 ∼ I G ( ν t − 1 , 0 , ι t − 1 , 0 ) . \sigma^2\mid y_{1:t-1},\tau\geq t-1\sim IG(\nu^{t-1,0},\iota^{t-1,0}). σ 2 ∣ y 1 : t − 1 , τ ≥ t − 1 ∼ I G ( ν t − 1 , 0 , ι t − 1 , 0 ) . After observing y t y_t y t ( β , σ 2 ) (\beta,\sigma^2) ( β , σ 2 ) [30 ] ; Fearnhead and Liu, [8 ] ):
Σ t , 0 = Σ t − 1 , 0 − A A T Q , \Sigma^{t,0}=\Sigma^{t-1,0}-AA^{T}Q, Σ t , 0 = Σ t − 1 , 0 − A A T Q , μ t , 0 = μ t − 1 , 0 + A e 0 , \mu^{t,0}=\mu^{t-1,0}+Ae_0, μ t , 0 = μ t − 1 , 0 + A e 0 , ν t , 0 = ν t − 1 , 0 + 1 2 , \nu^{t,0}=\nu^{t-1,0}+\frac{1}{2}, ν t , 0 = ν t − 1 , 0 + 2 1 , ι t , 0 = ι t − 1 , 0 + e 0 2 2 Q , \iota^{t,0}=\iota^{t-1,0}+\frac{e_0^2}{2Q}, ι t , 0 = ι t − 1 , 0 + 2 Q e 0 2 , where e 0 = y t − h 0 T μ t − 1 , 0 e_0=y_t-h_0^{T}\mu^{t-1,0} e 0 = y t − h 0 T μ t − 1 , 0 Q = h 0 T Σ t − 1 , 0 h 0 + 1 Q=h_0^{T}\Sigma^{t-1,0}h_0+1 Q = h 0 T Σ t − 1 , 0 h 0 + 1 A = Σ t − 1 , 0 h 0 / Q A=\Sigma^{t-1,0}h_0/Q A = Σ t − 1 , 0 h 0 / Q y t y_t y t y 1 : t − 1 y_{1:t-1} y 1 : t − 1
y t ∣ y 1 : t − 1 , τ ≥ t − 1 ∼ t 2 ν t − 1 , 0 ( h 0 T μ t − 1 , 0 , ι t − 1 , 0 ν t − 1 , 0 ( h 0 T Σ t − 1 , 0 h 0 + 1 ) ) , y_t\mid y_{1:t-1},\tau\geq t-1\sim t_{2\nu^{t-1,0}}\left(h_0^{T}\mu^{t-1,0},\frac{\iota^{t-1,0}}{\nu^{t-1,0}}\left(h_0^{T}\Sigma^{t-1,0}h_0+1\right)\right), y t ∣ y 1 : t − 1 , τ ≥ t − 1 ∼ t 2 ν t − 1 , 0 ( h 0 T μ t − 1 , 0 , ν t − 1 , 0 ι t − 1 , 0 ( h 0 T Σ t − 1 , 0 h 0 + 1 ) ) , where t 2 ν t − 1 , 0 t_{2\nu^{t-1,0}} t 2 ν t − 1 , 0 t t t 2 ν t − 1 , 0 2\nu^{t-1,0} 2 ν t − 1 , 0
A similar derivation works for the case where there is a change at time i i i
( β t ⊤ , γ t ⊤ ) ⊤ = ( β t − 1 ⊤ , γ t − 1 ⊤ ) ⊤ y t = h i ⊤ ( β t , γ t ) ⊤ + σ ϵ t . \begin{aligned}
(\beta_t^\top,\gamma_t^\top)^\top&= (\beta_{t-1}^\top,\gamma_{t-1}^\top)^\top\\
y_t &= h_i^\top(\beta_t,\gamma_t)^\top+\sigma\epsilon_t.
\end{aligned} ( β t ⊤ , γ t ⊤ ) ⊤ y t = ( β t − 1 ⊤ , γ t − 1 ⊤ ) ⊤ = h i ⊤ ( β t , γ t ) ⊤ + σ ϵ t . The model follows the same dynamic linear structures, but with an augmented parameter vector. As such, the posterior and predictive quantities follow immediately by applying the standard DLM update formulas introduced earlier.
Proof of Proposition 1
For candidate i i i j = i + 1 j=i+1 j = i + 1 σ 2 Σ i < σ 2 Σ j \sigma^2\Sigma_i < \sigma^2\Sigma_j σ 2 Σ i < σ 2 Σ j Figure 6 .
Figure 6. Illustration of two univariate Gaussian densities with σ 2 Σ i < σ 2 Σ j \sigma^2\Sigma_i < \sigma^2\Sigma_j σ 2 Σ i < σ 2 Σ j
Lemma C.1. Denote f i ( θ ) f_i(\theta) f i ( θ ) and f j ( θ ) f_j(\theta) f j ( θ ) as the density of N ( μ i , σ 2 Σ i ) \mathcal{N}(\mu_i,\sigma^2\Sigma_i) N ( μ i , σ 2 Σ i ) and N ( μ j , σ 2 Σ j ) \mathcal{N}(\mu_j,\sigma^2\Sigma_j) N ( μ j , σ 2 Σ j ) f i ( θ ) f_i(\theta) f i ( θ ) and f j ( θ ) f_j(\theta) f j ( θ ) intersect at
( c 1 , c 2 ) = μ i σ 2 Σ j − μ j σ 2 Σ i ∓ σ 2 Σ i σ 2 Σ j ( μ i − μ j ) 2 + ( σ 2 Σ j − σ 2 Σ i ) log σ 2 Σ j σ 2 Σ i σ 2 Σ j − σ 2 Σ i . (c_1,c_2)=\frac{\mu_i\sigma^2\Sigma_j-\mu_j\sigma^2\Sigma_i\mp\sqrt{\sigma^2\Sigma_i\sigma^2\Sigma_j}\sqrt{(\mu_i-\mu_j)^2+(\sigma^2\Sigma_j-\sigma^2\Sigma_i)\log\frac{\sigma^2\Sigma_j}{\sigma^2\Sigma_i}}}{\sigma^2\Sigma_j-\sigma^2\Sigma_i}. ( c 1 , c 2 ) = σ 2 Σ j − σ 2 Σ i μ i σ 2 Σ j − μ j σ 2 Σ i ∓ σ 2 Σ i σ 2 Σ j ( μ i − μ j ) 2 + ( σ 2 Σ j − σ 2 Σ i ) log σ 2 Σ i σ 2 Σ j . Proof. The intersection points satisfy f i ( θ ) = f j ( θ ) f_i(\theta)=f_j(\theta) f i ( θ ) = f j ( θ )
1 σ 2 Σ i exp ( − ( θ − μ i ) 2 2 σ 2 Σ i ) = 1 σ 2 Σ j exp ( − ( θ − μ j ) 2 2 σ 2 Σ j ) . \frac{1}{\sqrt{\sigma^2\Sigma_i}}\exp\left(-\frac{(\theta-\mu_i)^2}{2\sigma^2\Sigma_i}\right)=\frac{1}{\sqrt{\sigma^2\Sigma_j}}\exp\left(-\frac{(\theta-\mu_j)^2}{2\sigma^2\Sigma_j}\right). σ 2 Σ i 1 exp ( − 2 σ 2 Σ i ( θ − μ i ) 2 ) = σ 2 Σ j 1 exp ( − 2 σ 2 Σ j ( θ − μ j ) 2 ) . Taking logarithms we have
( θ − μ j ) 2 σ 2 Σ j − ( θ − μ i ) 2 σ 2 Σ i + log σ 2 Σ j σ 2 Σ i = 0. \frac{(\theta-\mu_j)^2}{\sigma^2\Sigma_j}-\frac{(\theta-\mu_i)^2}{\sigma^2\Sigma_i}+\log\frac{\sigma^2\Sigma_j}{\sigma^2\Sigma_i}=0. σ 2 Σ j ( θ − μ j ) 2 − σ 2 Σ i ( θ − μ i ) 2 + log σ 2 Σ i σ 2 Σ j = 0. Solving this quadratic equation in θ \theta θ c 1 c_1 c 1 c 2 c_2 c 2
Therefore, the total variation distance can be written as
T V ( f i ( θ ) , f j ( θ ) ) = 1 2 ∫ − ∞ ∞ ∣ f i ( θ ) − f j ( θ ) ∣ d θ = 1 2 [ F j ( c 1 ) − F i ( c 1 ) + F i ( c 2 ) − F i ( c 1 ) − F j ( c 2 ) + F j ( c 1 ) + F i ( c 2 ) − F j ( c 2 ) ] = F j ( c 1 ) − F j ( c 2 ) + F i ( c 2 ) − F i ( c 1 ) , \begin{aligned}
TV(f_i(\theta),f_j(\theta)) &= \frac{1}{2}\int_{-\infty}^{\infty}|f_i(\theta)-f_j(\theta)|d\theta\\
&=\frac{1}{2}\left[F_j(c_1)-F_i(c_1)+F_i(c_2)-F_i(c_1)-F_j(c_2)+F_j(c_1)+F_i(c_2)-F_j(c_2)\right]\\
&=F_j(c_1)-F_j(c_2)+F_i(c_2)-F_i(c_1),
\end{aligned} T V ( f i ( θ ) , f j ( θ )) = 2 1 ∫ − ∞ ∞ ∣ f i ( θ ) − f j ( θ ) ∣ d θ = 2 1 [ F j ( c 1 ) − F i ( c 1 ) + F i ( c 2 ) − F i ( c 1 ) − F j ( c 2 ) + F j ( c 1 ) + F i ( c 2 ) − F j ( c 2 ) ] = F j ( c 1 ) − F j ( c 2 ) + F i ( c 2 ) − F i ( c 1 ) , where F i ( θ ) = Φ ( θ − μ i σ 2 Σ i ) F_i(\theta)=\Phi\left(\frac{\theta-\mu_i}{\sqrt{\sigma^2\Sigma_i}}\right) F i ( θ ) = Φ ( σ 2 Σ i θ − μ i ) F j ( θ ) = Φ ( θ − μ j σ 2 Σ j ) F_j(\theta)=\Phi\left(\frac{\theta-\mu_j}{\sqrt{\sigma^2\Sigma_j}}\right) F j ( θ ) = Φ ( σ 2 Σ j θ − μ j )
[ Φ ( c 1 − μ j σ 2 Σ j ) − Φ ( c 2 − μ j σ 2 Σ j ) + Φ ( c 2 − μ i σ 2 Σ i ) − Φ ( c 1 − μ i σ 2 Σ i ) ] . (2) \left[\Phi\left(\frac{c_1-\mu_j}{\sqrt{\sigma^2\Sigma_j}}\right)-\Phi\left(\frac{c_2-\mu_j}{\sqrt{\sigma^2\Sigma_j}}\right)+\Phi\left(\frac{c_2-\mu_i}{\sqrt{\sigma^2\Sigma_i}}\right)-\Phi\left(\frac{c_1-\mu_i}{\sqrt{\sigma^2\Sigma_i}}\right)\right]. \tag*{(2)} [ Φ ( σ 2 Σ j c 1 − μ j ) − Φ ( σ 2 Σ j c 2 − μ j ) + Φ ( σ 2 Σ i c 2 − μ i ) − Φ ( σ 2 Σ i c 1 − μ i ) ] . (2) Priors for the Section 3.1
For the calibration experiment, the changepoint location follows a uniform distribution τ ∼ U ( 200 , 800 ) \tau\sim U(200,800) τ ∼ U ( 200 , 800 ) Pr ( τ ≥ 1000 ) = 0.9 \Pr(\tau\geq1000)=0.9 Pr ( τ ≥ 1000 ) = 0.9
For mean change with known variance case, data follow y t ∼ N ( β , 1 ) y_t\sim N(\beta,1) y t ∼ N ( β , 1 ) t ≤ τ t\leq\tau t ≤ τ y t ∼ N ( β + γ , 1 ) y_t\sim N(\beta+\gamma,1) y t ∼ N ( β + γ , 1 ) t > τ t>\tau t > τ ( β , γ ) ∼ N ( 0 , 0.25 2 [ 1 − 1 − 1 2 ] ) (\beta,\gamma)\sim N\left(0,0.25^2\begin{bmatrix}1&-1\\-1&2\end{bmatrix}\right) ( β , γ ) ∼ N ( 0 , 0.2 5 2 [ 1 − 1 − 1 2 ] ) ( β , γ ) ∼ N ( 0 , 0.25 2 × 10 [ 1 − 1 − 1 2 ] ) (\beta,\gamma)\sim N\left(0,0.25^2\times10\begin{bmatrix}1&-1\\-1&2\end{bmatrix}\right) ( β , γ ) ∼ N ( 0 , 0.2 5 2 × 10 [ 1 − 1 − 1 2 ] )
For mean change with unknown variance case, data follow y t ∼ N ( β , σ 2 ) y_t\sim N(\beta,\sigma^2) y t ∼ N ( β , σ 2 ) t ≤ τ t\leq\tau t ≤ τ y t ∼ N ( β + γ , σ 2 ) y_t\sim N(\beta+\gamma,\sigma^2) y t ∼ N ( β + γ , σ 2 ) t > τ t>\tau t > τ σ 2 ∼ I G ( 6 , 5 ) \sigma^2\sim IG(6,5) σ 2 ∼ I G ( 6 , 5 ) ( β , γ ) ∣ σ 2 ∼ N ( 0 , 0.25 2 σ 2 [ 1 − 1 − 1 2 ] ) (\beta,\gamma)\mid\sigma^2 \sim N\left(\mathbf{0},0.25^2\sigma^2\begin{bmatrix}1 & -1 \\ -1 & 2\end{bmatrix}\right) ( β , γ ) ∣ σ 2 ∼ N ( 0 , 0.2 5 2 σ 2 [ 1 − 1 − 1 2 ] ) ( β , γ ) ∣ σ 2 ∼ N ( 0 , 10 × 0.25 2 σ 2 [ 1 − 1 − 1 2 ] ) (\beta,\gamma)\mid\sigma^2 \sim N\left(\mathbf{0},10\times0.25^2\sigma^2\begin{bmatrix}1 & -1 \\ -1 & 2\end{bmatrix}\right) ( β , γ ) ∣ σ 2 ∼ N ( 0 , 10 × 0.2 5 2 σ 2 [ 1 − 1 − 1 2 ] )
For continuous change in slope with known variance case, data follow y t = β 1 + β 2 t ~ + ϵ t y_t=\beta_1+\beta_2\tilde{t}+\epsilon_t y t = β 1 + β 2 t ~ + ϵ t t ≤ τ t\leq\tau t ≤ τ ϵ t ∼ N ( 0 , 1 ) \epsilon_t\sim N(0,1) ϵ t ∼ N ( 0 , 1 ) t ~ \tilde{t} t ~ t ~ = t / T \tilde{t}=t/T t ~ = t / T y t = β 1 + β 2 t ~ + γ 2 ( t ~ − τ ~ ) + ϵ t y_t=\beta_1+\beta_2\tilde{t}+\gamma_2(\tilde{t}-\tilde{\tau})+\epsilon_t y t = β 1 + β 2 t ~ + γ 2 ( t ~ − τ ~ ) + ϵ t t > τ t>\tau t > τ
( β 1 , β 2 , γ 1 , γ 2 ) T ∼ N ( 0 , v [ 1 0 0 0 0 1 0 0 0 0 τ ~ 2 − τ ~ 0 0 − τ ~ 1 ] ) (\beta_1,\beta_2,\gamma_1,\gamma_2)^{\mathsf{T}}\sim N\left(\mathbf{0},v\begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \tilde{\tau}^2 & -\tilde{\tau} \\ 0 & 0 & -\tilde{\tau} & 1\end{bmatrix}\right) ( β 1 , β 2 , γ 1 , γ 2 ) T ∼ N 0 , v 1 0 0 0 0 1 0 0 0 0 τ ~ 2 − τ ~ 0 0 − τ ~ 1 where v = 0.25 2 v=0.25^2 v = 0.2 5 2 v = 10 × 0.25 2 v=10\times0.25^2 v = 10 × 0.2 5 2
For the discontinuous change in slope with known variance case, data follow y t = β 1 + β 2 t ~ + ϵ t y_t=\beta_1+\beta_2\tilde{t}+\epsilon_t y t = β 1 + β 2 t ~ + ϵ t t ≤ τ t\leq\tau t ≤ τ ϵ t ∼ N ( 0 , 1 ) \epsilon_t\sim N(0,1) ϵ t ∼ N ( 0 , 1 ) t ~ \tilde{t} t ~ t ~ = t / T \tilde{t}=t/T t ~ = t / T y t = ℓ + s ( t ~ − τ ~ ) + ϵ t y_t=\ell+s(\tilde{t}-\tilde{\tau})+\epsilon_t y t = ℓ + s ( t ~ − τ ~ ) + ϵ t t > τ t>\tau t > τ
( β 1 , β 2 , γ 1 , γ 2 ) T ∼ N ( 0 , v [ 1 0 − 1 0 0 1 0 − 1 − 1 0 2 + τ ~ 2 − τ ~ 0 − 1 − τ ~ 2 ] ) (\beta_1,\beta_2,\gamma_1,\gamma_2)^{\mathsf{T}}\sim N\left(\mathbf{0},v\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ -1 & 0 & 2+\tilde{\tau}^2 & -\tilde{\tau} \\ 0 & -1 & -\tilde{\tau} & 2\end{bmatrix}\right) ( β 1 , β 2 , γ 1 , γ 2 ) T ∼ N 0 , v 1 0 − 1 0 0 1 0 − 1 − 1 0 2 + τ ~ 2 − τ ~ 0 − 1 − τ ~ 2 where v = 0.25 2 v=0.25^2 v = 0.2 5 2 v = 10 × 0.25 2 v=10\times0.25^2 v = 10 × 0.2 5 2
For change in seasonality within the known-variance case, data follow y t = β s sin ( 2 π t / 24 ) + β c cos ( 2 π t / 24 ) + [ γ s sin ( 2 π t / 24 ) + γ c cos ( 2 π t / 24 ) ] 1 { t > τ } + ϵ t y_t=\beta_s\sin(2\pi t/24)+\beta_c\cos(2\pi t/24)+[\gamma_s\sin(2\pi t/24)+\gamma_c\cos(2\pi t/24)]1\{t>\tau\}+\epsilon_t y t = β s sin ( 2 π t /24 ) + β c cos ( 2 π t /24 ) + [ γ s sin ( 2 π t /24 ) + γ c cos ( 2 π t /24 )] 1 { t > τ } + ϵ t ϵ t ∼ N ( 0 , 1 ) \epsilon_t\sim N(0,1) ϵ t ∼ N ( 0 , 1 )
( β s , β c , γ s , γ c ) T ∼ N ( 0 , v [ 1 0 − 1 0 0 1 0 − 1 − 1 0 2 0 0 − 1 0 2 ] ) (\beta_s,\beta_c,\gamma_s,\gamma_c)^{\mathsf{T}}\sim N\left(\mathbf{0},v\begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ -1 & 0 & 2 & 0 \\ 0 & -1 & 0 & 2\end{bmatrix}\right) ( β s , β c , γ s , γ c ) T ∼ N 0 , v 1 0 − 1 0 0 1 0 − 1 − 1 0 2 0 0 − 1 0 2 where v = 0.25 2 v=0.25^2 v = 0.2 5 2 v = 10 ( 0.25 2 ) v=10(0.25^2) v = 10 ( 0.2 5 2 )
Priors and full results for the Section 3.2
Prior settings
In this simulation, we simulate data with a length of 2000 and have a changepoint at 1000. The prior probability of no change is set at P r ( τ ≥ 2000 ) = 0.9 Pr(\tau\geq2000)=0.9 P r ( τ ≥ 2000 ) = 0.9
For the change-in-mean with known variance case, the data follow N ( 0 , 1 ) N(0,1) N ( 0 , 1 ) N ( 0.25 , 1 ) N(0.25,1) N ( 0.25 , 1 ) ( β , γ ) (\beta,\gamma) ( β , γ ) ( 0 , 0 ) (0,0) ( 0 , 0 ) 0.25 2 [ 1 − 1 − 1 2 ] 0.25^2 \begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix} 0.2 5 2 [ 1 − 1 − 1 2 ] 10 × 0.25 2 [ 1 − 1 − 1 2 ] 10 \times0.25^2 \begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix} 10 × 0.2 5 2 [ 1 − 1 − 1 2 ]
For the change-in-mean with unknown variance case, the data follow N ( 0 , 1 ) \mathcal{N}(0,1) N ( 0 , 1 ) N ( 0.25 , 1 ) \mathcal{N}(0.25,1) N ( 0.25 , 1 ) ( β , γ ) (\beta,\gamma) ( β , γ ) ( 0 , 0 ) (0,0) ( 0 , 0 ) v σ 2 [ 1 − 1 − 1 2 ] v\sigma^2 \begin{bmatrix} 1 & -1 \\ -1 & 2 \end{bmatrix} v σ 2 [ 1 − 1 − 1 2 ] σ 2 ∼ IG ( 6 , 5 ) \sigma^2 \sim\operatorname{IG}(6,5) σ 2 ∼ IG ( 6 , 5 ) v = 0.25 2 v=0.25^2 v = 0.2 5 2 v = 10 × 0.25 2 v=10 \times0.25^2 v = 10 × 0.2 5 2
For the continuous slope change with known variance case, data follow N ( 0 , 1 ) \mathcal{N}(0,1) N ( 0 , 1 ) 1.73 ( t ~ − 0.5 ) + N ( 0 , 1 ) 1.73(\tilde{t}-0.5)+\mathcal{N}(0,1) 1.73 ( t ~ − 0.5 ) + N ( 0 , 1 ) t ~ = t / T \tilde{t}=t/T t ~ = t / T ( β 1 , β 2 , γ 1 , γ 2 ) (\beta_1,\beta_2,\gamma_1,\gamma_2) ( β 1 , β 2 , γ 1 , γ 2 ) ( 0 , 0 , 0 , 0 ) (0,0,0,0) ( 0 , 0 , 0 , 0 ) [ 0.25 2 0 0 0 0 0.25 2 0 0 0 0 τ ~ 2 1.73 2 − τ ~ 1.73 2 0 0 − τ ~ 1.73 2 1.73 2 ] \begin{bmatrix} 0.25^2 & 0 & 0 & 0 \\ 0 & 0.25^2 & 0 & 0 \\ 0 & 0 & \tilde{\tau}^2 1.73^2 & -\tilde{\tau}1.73^2 \\ 0 & 0 & -\tilde{\tau}1.73^2 & 1.73^2 \end{bmatrix} 0.2 5 2 0 0 0 0 0.2 5 2 0 0 0 0 τ ~ 2 1.7 3 2 − τ ~ 1.7 3 2 0 0 − τ ~ 1.7 3 2 1.7 3 2
For the discontinuous slope change with known variance case, data follow N ( 0 , 1 ) \mathcal{N}(0,1) N ( 0 , 1 ) 0.15 + 1.65 ( t ~ − 0.5 ) + N ( 0 , 1 ) 0.15+1.65(\tilde{t}-0.5)+\mathcal{N}(0,1) 0.15 + 1.65 ( t ~ − 0.5 ) + N ( 0 , 1 ) t ~ = t / T \tilde{t}=t/T t ~ = t / T ( β 1 , β 2 , γ 1 , γ 2 ) (\beta_1,\beta_2,\gamma_1,\gamma_2) ( β 1 , β 2 , γ 1 , γ 2 ) ( 0 , 0 , 0 , 0 ) (0,0,0,0) ( 0 , 0 , 0 , 0 ) [ 0.15 2 0 − 0.15 2 0 0 1.65 2 0 − 1.65 2 − 0.15 2 0 2 × 0.15 2 + 1.65 2 τ ~ 2 − τ ~ 1.65 2 0 − 1.65 2 − τ ~ 1.65 2 2 × 1.65 2 ] \begin{bmatrix} 0.15^2 & 0 & -0.15^2 & 0 \\ 0 & 1.65^2 & 0 & -1.65^2 \\ -0.15^2 & 0 & 2 \times0.15^2+1.65^2\tilde{\tau}^2 & -\tilde{\tau}1.65^2 \\ 0 & -1.65^2 & -\tilde{\tau}1.65^2 & 2 \times1.65^2 \end{bmatrix} 0.1 5 2 0 − 0.1 5 2 0 0 1.6 5 2 0 − 1.6 5 2 − 0.1 5 2 0 2 × 0.1 5 2 + 1.6 5 2 τ ~ 2 − τ ~ 1.6 5 2 0 − 1.6 5 2 − τ ~ 1.6 5 2 2 × 1.6 5 2
For the seasonality change with known variance case, data follow y t = sin ( 2 π t / 24 ) + cos ( 2 π t / 24 ) + [ 0.25 sin ( 2 π t / 24 ) + 0.25 cos ( 2 π t / 24 ) ] 1 { t > τ } + ϵ t y_t=\sin(2\pi t/24)+\cos(2\pi t/24)+[0.25\sin(2\pi t/24)+0.25\cos(2\pi t/24)]1\{t>\tau\}+\epsilon_t y t = sin ( 2 π t /24 ) + cos ( 2 π t /24 ) + [ 0.25 sin ( 2 π t /24 ) + 0.25 cos ( 2 π t /24 )] 1 { t > τ } + ϵ t ϵ t ∼ N ( 0 , 1 ) \epsilon_t\sim\mathcal{N}(0,1) ϵ t ∼ N ( 0 , 1 ) ( β s , β c , γ s , γ c ) (\beta_s,\beta_c,\gamma_s,\gamma_c) ( β s , β c , γ s , γ c ) ( 1 , 1 , 0 , 0 ) (1,1,0,0) ( 1 , 1 , 0 , 0 ) v [ 1 0 − 1 0 0 1 0 − 1 − 1 0 2 0 0 − 1 0 2 ] v\begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \\ -1 & 0 & 2 & 0 \\ 0 & -1 & 0 & 2 \end{bmatrix} v 1 0 − 1 0 0 1 0 − 1 − 1 0 2 0 0 − 1 0 2 v = 0.25 2 v=0.25^2 v = 0.2 5 2 v = 10 ( 0.25 2 ) v=10(0.25^2) v = 10 ( 0.2 5 2 )
Results
Table 4 provides the full results for five examples under both main and diffuse priors here. Denote V = { r : τ < T a r ≤ 2000 } V=\{r:\tau<T_a^r\le2000\} V = { r : τ < T a r ≤ 2000 } P e P_e P e R R R ( 1 − ∣ V ∣ R ) (1-\frac{|V|}{R}) ( 1 − R ∣ V ∣ ) ( 1 ∣ V ∣ ∑ r ∈ V ( T a r − τ ) ) (\frac{1}{|V|}\sum_{r\in V}(T_a^r-\tau)) ( ∣ V ∣ 1 ∑ r ∈ V ( T a r − τ )) ( median r ∈ V ∣ τ ^ a r ( T a r ) − τ ^ e r ( T a r ) ∣ ) (\operatorname{median}_{r\in V}|\hat{\tau}_a^r(T_a^r)-\hat{\tau}_e^r(T_a^r)|) ( median r ∈ V ∣ τ ^ a r ( T a r ) − τ ^ e r ( T a r ) ∣ ) ( 1 ∣ V ∣ ∑ r ∈ V H a r ( T a r ) ∩ H e r ( T a r ) H a r ( T a r ) ∪ H e r ( T a r ) ) (\frac{1}{|V|}\sum_{r\in V}\frac{H_a^r(T_a^r)\cap H_e^r(T_a^r)}{H_a^r(T_a^r)\cup H_e^r(T_a^r)}) ( ∣ V ∣ 1 ∑ r ∈ V H a r ( T a r ) ∪ H e r ( T a r ) H a r ( T a r ) ∩ H e r ( T a r ) ) ( 100 ∣ V ∣ ∑ r ∈ V ∣ P e ( H a r ( T a r ) ∣ y 1 : T a r ) − 0.95 ∣ ) (\frac{100}{|V|}\sum_{r\in V}|P_e(H_a^r(T_a^r)\mid y_{1:T_a^r})-0.95|) ( ∣ V ∣ 100 ∑ r ∈ V ∣ P e ( H a r ( T a r ) ∣ y 1 : T a r ) − 0.95∣ )
Table 4. Simulation result on a fixed horizon T = 2000 T = 2000 T = 2000 τ = 1000 \tau= 1000 τ = 1000
Computational results for Section 3.3
Here we calculate the median of computational speed at each step of Bagel. The main code is implemented in C++, and it is invoked per step from R. The code was executed on a desktop machine with a 13th-generation Intel(R) Core(TM) i7, featuring 16 physical cores with a 3.4 GHz base clock, and 64 GB of RAM, running on Windows 11 and R 4.4.3 with a single thread. The running time is recorded through the system wall clock read immediately before and after each call to the recursion, so that the cost of every step is measured separately.
Figure 7 shows the median running time of Bagel when varying the budget M M M d 1 = d 2 = 1 d_1 = d_2 = 1 d 1 = d 2 = 1 M = 100 M = 100 M = 100
Figure 7. Median of per-step runtime of the Bagel over 200 replications. The left figure shows the running time for different M M M d 1 = d 2 = 1 d_1 = d_2 = 1 d 1 = d 2 = 1 ( d 1 , d 2 ) (d_1,d_2) ( d 1 , d 2 ) M = 100 M = 100 M = 100 T = 400 T = 400 T = 400 T = 1000 T = 1000 T = 1000
There is a short period during which Bagel can be slower than the exact approach when the number of particles reaches its maximum length M M M t t t
