# Maximizing Log-Likelihood Estimation for Changepoint Detection

I'm trying to code the changepoint detection algo described here:

Link to original ppt: http://www.slideshare.net/kuma0177/velocity-ny-2014v5-39160794
Slides: 13-16


Slide 16: http://imgur.com/fF69Q3Q

It is my understanding that the slides use a rolling median as input to calculate lambda on slide 16. However, I'm not entirely sure how you plug the median into that lambda equation. My stats background is pretty weak.

My guess is that this algorithm is not an online one, but a batch one. I think what happens is that you have a time series of medians (x0,y0), (x1,y1), ... (xn,yn) then you split this into 2 time series in all possible ways:

change_point_Tau1@(x0,y0), (x1,y1), (x2,y2), ... (xn,yn)
(x0,y0), change_point_Tau1@(x1,y1), (x2,y2), ... (xn,yn)
(x0,y0), (x1,y1), change_point_Tau1@(x2,y2), ... (xn,yn)
...
...
(x0,y0), (x1,y1), (x2,y2), ... change_point_Tau1@(x(n-1),y(n-1)), (xn, yn)
(x0,y0), (x1,y1), (x2,y2), ... (x(n-1),y(n-1)), change_point_Tau1@(xn, yn)


Then for each of these 2 halves you want to calculate some Theta. Since we can't get theta we use an approximation for it, shown as Theta hat. Theta is sth called a probability density function. I'd need this function to be simpler enough so that it can be approximated using a programming language like Scala or Java.

Q: How do you calculate this theta hat? There are 2 of them in equation ML.

Q: Is my understanding of this algorithm correct?

Q: How to calculate all those probabilities?


For example I see P(y(1:Tau1) | Theta_hat) - probability no change points exists at Tau1 given Theta hat. Using Bayes, this becomes:

P(y(1:Tau1) | Theta_hat) = P(Theta_hat | y(1:Tau1))*P(y(1:Tau1)) / P(Theta_hat)


Now where would I get those new probabilities like P(Theta_hat | y(1:Tau1)), P(y(1:Tau1)) and P(Theta_hat)?

I know, lots of questions all crammed in here..

• Hello, you could also try this package in R changepoint – forecaster Nov 11 '14 at 16:22

It is not clear from the presentation what distributional assumptions are being made in order to calculate the likelihoods.

It might be simpler for you to look at the recently published BreakoutDetection package published by the same authors: https://blog.twitter.com/2014/breakout-detection-in-the-wild.

But if you are more interested in learning about changepoint detection and how to use likelihoods then read on.

Normally you have a time series y which we assumes has n observations. In order to use likelihoods we need to make some assumptions about the distribution that y comes from. It is often assumed that the data come from Normal distribution (although this isn't always appropriate). Following a distributional assumption you need to decide which parameters of the distribution are allowed to change, e.g. mean, variance, both.

If all the parameters change then you proceed by splitting your data into 2 halves, before change and after change, and use maximum likelihoods to fit the parameters to each half. In this way it is like a normal analysis where you just have data points and you fit a model to that data. If not all parameters can change then you need to estimate those that don't change from the whole data (not always an easy task, especially when parameters are linked).

The trick with changepoint analysis is that you don't know where the change is, so you have to calculate the likelihood for each possible changepoint location and take the most likely as the hypothesized changepoint location. It is this location and likelihood that you then test to see if the change is significant by comparing the likelihood ratio to a threshold to see if the change is significant (comparing to c in slide 16).

• The code has now been published to R and on github. See github.com/twitter/BreakoutDetection Though I find the code does not work that well. I'm experimenting w trying a moving window median and then finding the StDev of last 24 medians and creating 2 bands around the median: med+1.5StDev, med-1.5StDev. When bands diverge I will mark that point as change point. – Adrian Nov 11 '14 at 21:52

Not sure if there is still interest in answering this question, anyways here is my take on this.

Adrian, you are right that what is in the slide above is for a given sample. For cases in which the sample size varies in time we can talk about an "online" changepoint detection algorithm. The big difference is that an online changepoint detection algorithm has to take into account the multiple testing problem.

But we are not in that situation here, which makes things easier.

Consistent with the slide above and with the notation of the R changepoint package (for the documentation, see here https://cran.r-project.org/web/packages/changepoint/changepoint.pdf), let's define $\tau$ as the changepoint time that we want to test. Each data point in the time series is assumed to be drawn from some probability distribution function which is fully described by the quantity $\theta$ (in general, a set of parameters). For example, $\theta$ could be the probability $p$ of success in a binomial distribution, or the mean and variance in a normal distribution.

At this point the test goes like this: the null hypothesis is that there is no changepoint, while the alternative hypothesis assumes that there is a changepoint at the time $t=\tau$. More formally, here is our hypothesis test: \begin{eqnarray} H_0 &:& \theta_1=\theta_2=\dots=\theta_{N-1}=\theta_N \nonumber \\ H_1 &:& \theta_1=\theta_2=\dots=\theta_{\tau-1}=\theta_\tau\neq\theta_{\tau+1}=\theta_{\tau+2}=\dots=\theta_{N-1}=\theta_N \nonumber \end{eqnarray}

The key in the expression $H_1$ above is in the inequality $\theta_\tau\neq\theta_{\tau+1}$: at some point in the time series, and precisely between $t=\tau$ and $t=\tau+1$, the underlying distribution changes.

The algorithm for detection is based on the log-likelihood ratio. Let's first define what is the likelihood. The likelihood is nothing else than the probability of observing the data that we have (in the time series), assuming that the null (or the alternative) hypothesis are true. It is a measure of how good the hypothesis is: the highest the likelihood, the higher the data are well fit by the $H_0$ (or $H_1$) assumption.

Assuming independent random variables, under the null hypothesis $H_0$ the likelihood $\mathcal{L}(H_0)$ is given by the probability of observing the data $\mathbf{x}=x_1,\dots,x_N$ conditional on $H_0$. In other words, \begin{equation} \mathcal{L}(H_0)=p(\mathbf{x}|H_0)=\prod_{i=1}^{N}p(x_i|\theta_0) \label{eq:L1} \end{equation}

Now let us define the likelihood of the alternative hypothesis, \begin{equation} \mathcal{L}(H_1)=p(\mathbf{x}|H_1)=\prod_{i=1}^{\tau}p(x_i|\theta_1)\prod_{j=\tau+1}^{N}p(x_j|\theta_2) \label{eq:L2} \end{equation}

Note that in the above expressions I have written $\theta_0$ as the set of parameters that define the distribution under $H_0$ and with $\theta_{1,2}$ the parameters that define the distribution under $H_1$ before and after the changepoint, respectively.

The log-likelihood ratio $\mathcal{R}_\tau$ is then \begin{equation} \mathcal{R}_\tau=\log\left(\frac{\mathcal{L}_{H_1}}{\mathcal{L}_{H_0}}\right)=\sum_{i=1}^{\tau}\log p(x_i|\theta_1) + \sum_{j=\tau+1}^{N}\log p(x_j|\theta_2) - \sum_{k=1}^{N}\log p(x_k|\theta_0) \label{eq:LRatio} \end{equation}

Since $\tau$ is not known, the previous equation becomes a function of $\tau$. We can then define a generalized log-likelihood ratio $G$, which is the maximum of $\mathcal{R}_\tau$ for all the possible values of $\tau$, \begin{equation} G = \max_{1\leq\tau\leq N}\mathcal{R}_\tau \label{eq_likelihood} \end{equation}

If the null hypothesis is rejected, then the maximum likelihood estimate of the changepoint is the value $\hat{\tau}$ that maximizes the generalized likelihood ratio, \begin{equation} \hat{\tau} = \underset{1\leq\tau\leq N}{\mathrm{argmax}} \ \mathcal{R}_\tau \end{equation}

In general, the null hypothesis is rejected for a sufficiently large value of $G$. In other words, there is a critical value $\lambda^*$ such that $H_0$ is rejected if \begin{equation} \bbox[lightblue,5px,border:2px solid red] { 2G=2R(\hat{\tau})>\lambda^* } \label{eq:criterion} \end{equation}

The factor 2 is retained to be consistent with the $\texttt{changepoint}$ package.

The problem is how to define this critical value $\lambda^*$. The package $\texttt{changepoint}$ has several ways of defining the "penalty" factor $\lambda^*$. Going through the R code, I managed to find the definition of a few of them,

• BIC (Bayesian Information Criterion): $\lambda^* = k \log n$
• MBIC (Modified Bayesian Information Criterion): $\lambda^* = (k+1)\log n + \log(\tau) + \log(n-\tau+1)$. NOTE: this is the only one that does not seem to be consistent with other definitions of MBIC that I found in the literature.
• AIC (Akaike Information Criterion): $\lambda^* = 2k$
• Hannan-Quinn: $\lambda^* = 2k\log(\log n)$

where $k$ is the number of extra parameters that are added as a result of defining a changepoint. For example, it is $k=1$ if there is just a shift in the mean or a shift in the variance.

### EXAMPLE. Normally distributed random variables - change in mean

Assume the variables that compose the time series are drawn from independent normal random distributions. We want to test the hypothesis that there is a change in the mean of the distribution at some discrete point in time $\tau$, while we assume that the variance $\sigma^2$ does not change. The probability density function is

\begin{equation} f(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}\textrm{e}^{-(x-\mu)^2/2\sigma^2} \nonumber \end{equation}

Let us call $\mu_1$ the mean before the changepoint, $\mu_2$ the mean after the changepoint and $\mu_0$ the global mean.

Under the null hypothesis there is no change in mean so the likelihood looks like

\begin{equation} \mathcal L_{H_0} = \frac{1}{\sqrt{2\pi\sigma^{2N}}}\prod_{i=1}^{N}\exp\left[-\frac{(x_i-\mu_0)^2}{2\sigma^2}\right] \label{eq:L3} \end{equation}

Under the alternative hypothesis, a changepoint occurs at the time $\tau$ and the corresponding likelihood will then be

\begin{equation} \mathcal L_{H_1} = \frac{1}{\sqrt{2\pi \sigma^{2N}}}\prod_{i=1}^{\tau}\exp\left[-\frac{(x_i-\mu_1)^2}{2\sigma^2}\right] \prod_{j=\tau+1}^{N}\exp\left[-\frac{(x_j-\mu_2)^2}{2\sigma^2}\right]\ . \label{eq:L4} \end{equation}

Now it is time to write the log-likelihood ratio, obtaining \begin{equation} \bbox[white,5px,border:2px solid red]{ \mathcal{R}_\tau=\log\left(\frac{\mathcal{L}_{H_1}}{\mathcal{L}_{H_0}}\right)=-\frac{1}{2\sigma^2}\left[\sum_{i=1}^{\tau}(x_i-\mu_1)^2+\sum_{j=\tau+1}^{N}(x_j-\mu_2)^2-\sum_{k=1}^{N}(x_k-\mu_0)^2\right] } \label{eq:LR1} \end{equation}

NOTE: this answer is based on a longer post that can be found at http://www.claudiobellei.com/2016/11/15/changepoint-frequentist/

• Thank you for your answer! While I solved my change point in some other way a while ago I like your post a lot - in fact I did use log likelihood. It makes me want to experiment a bit now :) – Adrian Dec 2 '16 at 14:19
• Thank you Adrian, glad to be of some help to the community :) – cbellei Dec 2 '16 at 16:38