A message from our CEO about the future of Stack Overflow and Stack Exchange. Read now.
21

There are 3 main functions in the changepoint package, cpt.mean, cpt.var and cpt.meanvar. As a practitioner these are the only functions in the package that you should need. If you think that your data may contain a change in mean then you use the cpt.mean function, etc. The next question you should ask yourself if whether you are looking for a single or ...


20

Another approach would be to wrap the call to lmer in a function that is passed the breakpoint as a parameter, then minimize the deviance of the fitted model conditional upon the breakpoint using optimize. This maximizes the profile log likelihood for the breakpoint, and, in general (i.e., not just for this problem) if the function interior to the wrapper (...


17

You could use time series outlier detection to detect changes in time series. Tsay's or Chen and Liu's procedures are popular time series outlier detection methods . See my earlier question on this site. R's tsoutlier package uses Chen and Liu's method for detection outliers. SAS/SPSS/Autobox can also do this. See below for the R code to detect changes in ...


14

It appears you are looking for spikes within intervals of relative quiet. "Relative" means compared to typical nearby values, which suggests smoothing the series. A robust smooth is desirable precisely because it should not be influenced by a few local spikes. "Quiet" means variation around that smooth is small. Again, a robust estimate of local ...


11

If the observations of your time series data are correlated with the immediately previous observations, the paper by Chen and Liu (1993)$^{[1]}$ may interest you. It describes a method to detect level shifts and temporary changes in the framework of autoregressive moving-average time series models. [1]: Chen, C. and Liu, L-M. (1993), "Joint Estimation of ...


10

For simplicity, I would suggest analyzing the sizes (absolute values) of the residuals relative to a robust smooth of the data. For automated detection, consider replacing those sizes by an indicator: 1 when they exceed some high quantile, say at level $1-\alpha$, and 0 otherwise. Smooth this indicator and highlight any smoothed values that exceed $\alpha$....


9

If you are wanting to test "significance" then I suggest you use the Asymptotic penalty option, i.e. penalty='Asymptotic' and pen.value=0.05 for 95% confidence. This automatically sets the penalty based on the cost function you are using. I find that this works well for smaller data sets <1000 but not too small <100. If you want to use the manual ...


8

In general, it's a bit odd to want to fit something as piece-wise linear. However, if you really wish to do so, then the MARS algorithm is the most direct. It will build up a function one knot at a time; and then usually prunes back the number of knots to combat over-fitting ala decision trees. You can access the MARS algotithm in R via earth or mda. In ...


8

Would MARS be applicable? R has the package earth that implements it.


7

You can try out the changefinder library on PyPI. The description says that it's an online Change Detection Library based on the ChangeFinder algorithm There are also some Python implementations of Michele Basseville's Statistical Change Point Detection techniques available in tutorial format on this Github repo.


7

I would approach this problem from the following perspectives. These are just some ideas off the top of my head - please take them with a grain of salt. Nevertheless, I hope that this will be useful. Time series clustering. For example, by using popular dynamic time warping (DTW) or alternative approaches. Please see my related answers: on DTW for ...


7

Some remarks: You need to estimate two parameters in each segment (intercept and slope). Hence breakpoints() requires that there are at least three observations in each segment...otherwise you cannot estimate the parameters (without getting a perfect fit). But three is the minimal value that is technically possible. Whether or not it leads to meaningful ...


7

The two good papers on this subject are below: 1) Bayesian Online Change Point Detection 2) Modeling changing dependency structure in multivariate time series These do not apply a clustering algorithm but take the interval (since the last change point) into account as you have asked for. And they work with parametric distributions. The paper by Adams and ...


6

This will be an R-centric answer. One approach is to wrap the call to lm in a function which is passed the breakpoint and constructs a regression conditional upon that breakpoint, then minimize the deviance of the fitted model conditional upon the breakpoint by just iterating over the possible values for the breakpoint. This maximizes the profile log ...


6

A "Jump" in a time series is a permanent change in the equation's implied intercept. Consider a series 1,1,1,1,4,4,4,4,4 where the basic model is y(t)=1 + 3*x(t) and x is the level shift/step shift series 0,0,0,0,1,1,1,1,1 . Thus suggests an intercept change at period 5 from a "1" to a "4" . If a second series is 1,2,3,4,8,9,10,11,12 then the model is [1-B]y(...


6

There are still some gaps in the Python library for using advanced statistics packages. Have you tried using the RPy module? When using RPy you can load R modules. brief tutorial on RPy: http://www.sciprogblog.com/2012/08/using-r-from-within-python.html strucchange


6

SeanEaster has some good advice. Bayes factor can be difficult to compute, but there are some good blog posts specifically for Bayes factor in PyMC2. A closly related question is goodness-of-fit of a model. A fair method for this is just inspection - posteriors can give us evidence of goodness-of-fit. Like quoted: "Had no change occurred, or had the ...


6

This problem in Stats is referred to as the (univariate) Temporal Event Detection. The simplest idea is to use a moving average and standard deviation. Any reading that is "out of" 3-standard deviations (rule-of-thumb) is considered an "event". There of course, are more advanced models that use HMMs, or Regression. Here is an introductory overview of the ...


6

My response using AUTOBOX is quite similar to @forecaster but with a much simpler model. Box and Einstein and others have reflected on keeping solutions simple but not too simple. The model that was automatically developed was . The actual and cleansed plot is very similar . A plot of the residuals (which should always be shown ) is here along with the ...


6

I think you are on the right track and I guess that in principle the changepoint package should be usable. You can simply order your response variable (e.g., a species abundance) by the gradient of interest (e.g., water temperature) and then apply the functions to the as if they were a time series. You may need to map the changepoint indexes back to the ...


6

Your model implicitly is this. You have data $(x_1,y_1), \ldots, (x_m,y_m)$ and other data $(x_{m+1},y_{m+1}), \ldots, (x_n,y_n)$ for which $x_j \ge x_i$ whenever $m \lt j \le n$ and $1 \le i \le m$: the first data set is the "left portion" and the second is the "right portion". You wish to estimate coefficients $\alpha, \alpha^\prime, \beta, \beta^\prime$ ...


5

The solution proposed by jbowman is very good, just adding a few theoretical remarks: Given the discontinuity of the indicator function used, the profile-likelihood might be highly erratic, with multiple local minima, so usual optimizers might not work. The usual solution for such "threshold models" is to use instead the more cumbersome grid search, ...


5

This sounds to me more like a change point problem. You should investigate the bcp package, http://cran.r-project.org/web/packages/bcp/index.html, as well as strucchange, http://cran.r-project.org/web/packages/strucchange/index.html John Emerson maintains the bcp package, and when I met him, he seemed very willing to engage the community and potential users:...


5

There is a lot of literature for testing the change in mean. If it is known that mean does not change, and you need to test the variance, you can convert the problem of testing for change in variance to the one of testing for change in mean with simple transformation. Suppose your initial data is $X_i$, then define $Y_i=(X_i-\mu)^2$, where $\mu$ is the ...


5

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 ...


5

Both (1) and (1b) are correct. The OP has it right that (in this model) there might be a changepoint at $t+1$, and $x_{t+1}$ depends on whether there is a changepoint. This does not imply any problems with (1) as the possible values of $r_{t+1}$ are fully "covered" by $P(x_{t+1} \mid r_t, x_{1:t})$. $P(x_{t+1} | r_t, x_{1:t})$ means the conditional ...


4

To quote from the intro in the Efron/Tibshirani text on the subject: The message of this book can therefore be summarized by paraphrasing Tukey: "The bootstrap, like a shotgun, can blow the head off any problem if the statistician can stand the resulting mess''. It sounds like you're trying to make inference based on the assumption that the ...


4

If you are using the changepoint package (which from the output I presume you are) then the changepoint locations always end with n. Thus the length of your data is 118 and there is a single change in mean at 15. For info AMOC means At Most One Change and thus the maximum number of changepoints the method will identify is 1. As 1 has been identified then ...


4

This implementation of the Python package rpy2 worked for me: import numpy as np from rpy2.robjects.packages import importr import rpy2.robjects as robjects r = robjects.r #allows access to r object with r. bcp = importr('bcp') #import bayesian change point package in python values = bcp.bcp( r.c( r.rnorm(50) , r.rnorm(50,5,1), r.rnorm(50) ) ) #use bcp ...


4

Your use case is similar to a hidden-markov-model (HMM) but not quite since your state names are known and your training data is already labelled. Based on your labelled dataset, you can learn not only the transition matrix between states (stable/unstable) but also the observation probabilities per state. In other words, given [speed1=3, speed2=2, ...


Only top voted, non community-wiki answers of a minimum length are eligible