# Step Change Detection

I am using a non linear least squares method to fit an analytical function to some experimental data. I have to provide some initial guess values to the algorithm, so I am trying to figure out how to do this automatically (rather than by eye, which is what I have been doing).

This is some simulated data, created by adding normally distributed random noise to the analytical function I am trying reliably detect the position of this step change in the data. I have had some limited success by calculating the mean variance in the data points and looking for points in the data that differ significantly from this value, but this approach seems very limited by the signal to noise ratio.

I am hoping for some direction on what I need to look into to solve my problem, as I don't know much statistics at all.

Thank you!

-Edit paste bin link to xy data

http://pastebin.com/QTawFex3

• Maybe it would be useful to consider change-point detection techniques. Circular binary segmentation might be just what you need. – Mur1lo Aug 15 '16 at 18:22
• Chen and Liu 1993 describe a method of detecting level shift in a time series. Their algorithm is implemented in the tsoutliers R package. – pionpi_ Aug 16 '16 at 13:34
• A new powerful change detection technique is impulse-indicator saturation and step-indicator saturation in the works of Castle, Doornik, Hendry & Co. (Google these keywords to find more.) – Richard Hardy Aug 16 '16 at 18:33
• The best reply I could find on that is there, works very well in my case and looking at your data it should work as well. stackoverflow.com/questions/48000663/… – francois-xavier sireta Oct 4 '18 at 16:19

There may be some more sophisticated methods for this but here is my first thought.

You basically want to take the derivative of the function and find where it is the largest. Numerically, you can just take the difference between data points and find which two points have the biggest difference. Then the midpoint of the x-values for these two points is your location of biggest change.

This simple method is susceptible to the noise. So you can first filter the data using a filter that does not shift the data to the right or left. If you use a simple FIR filter, then filter front to back and then filter the result from back to front. The result is a doubly filtered and NON shifted data set. Then follow the procedure above to find the point with the largest difference between values.

You can also use more sophisticated numerical differential calculations that use more then the difference of two points.

• Thank you for the answer. I was hung up on the fact that this was discrete data, and didn't think to just differentiate it! I'm going to look into using the FIR filter on some of my smaller peaks, to see if this works – jm22b Aug 15 '16 at 14:13
• I have implemented your idea and it works like a charm, thanks a lot! – jm22b Aug 16 '16 at 9:40
• That is great to hear. Thank you for letting me know. – jason Aug 16 '16 at 21:10

Your data visually suggests an asymptotic (gradual) change to the new level. Time series methods can often be used to detect these kinds of structures even if the data is not time series. Please post your data and I may be able to demonstrate this with "toys" at my disposal. If your data is time series then as @jason reflected one needs to deal effectively with the noise model to correctly "see" the structure.

EDITED UPON RECEIPT OF DATA:

Modelling is often an iterative approach with interim steps providing valuable clues to a useful model. I took your data and introduced it to AUTOBOX (one of my toys which I have helped develop). An initial graph strongly suggested a longitudinal (chronological) data set where the X series is reported at fixed intervals. AUTOBOX automatically suggested a standard ARIMA model (with Intervention Detection) replacing the non-stationary X with a differencing operator. Here is the actual/fit/forecast graph and the suggested model.  Upon examination another possible model incorporating a lag structure for an indicator variable suggested itself. I introduced a Pulse at time period 76 (a Dynamic Predictor expressly allowing up to a possible lag effect of 50 periods) (the beginning of the transition) to deal with the relationship between the original Y and the user-suggested X in order to more fully investigate the effect of X than accept the total setting-aside of X.

Following is the actual- fit-forecast graph for that approach and the identified robust transfer function model . with residual plot and residual acf here The final model captures the dynamics in certain lags of the Dynamic Predictor and a few pulses and a reasonable memory structure.

Even the most powerful analysis packages often need some guidance when dealing with complex real world data sets like this one as nothing compares to the creative human mind.

• What's the best way to share the data? – jm22b Aug 15 '16 at 15:03
• I have edited my question to include a pastebin link containing some data – jm22b Aug 15 '16 at 15:15
• I am just curious ...1) What was you analytic function and 2) what was the variance of the error series (your normally distributed random noise) : AUTOBOX delivered : Variance =SOS/(n) .402675E-04 Adjusted Variance =SOS/(n-m) .443625E-04 Standard Deviation RMSE =SQRT(Adj Var) .666052E-02 – IrishStat Aug 16 '16 at 10:23
• My analytical function is $T(\lambda) = C_1 + C_2 [\mathrm{erfc(\frac{\lambda_0 - \lambda}{\sqrt{2} \sigma})} - \mathrm{Exp}(\frac{\lambda_0 - \lambda}{\tau} + \frac{\sigma^2}{2 \tau^2}) * erfc(\frac{\lambda_0 - \lambda}{\sqrt{2} \sigma} + \frac{\sigma}{\sqrt{2} \tau})]$ Where, $C_1, C_2, \lambda_0, \sigma, \tau$ are parameters to be fitted The image of the plot in my original post was simply this function + normally distributed random noise with a stdev of 0.1 (I think). The data I shared with you was real data, and I haven't the foggiest what it's properties are! – jm22b Aug 16 '16 at 10:30

One technique is to test all values of the x variable for the standard deviation of the data before and after it. For a true step function, the sum of those two will be minimal at the step location, and the minimum should be a good starting parameter for your nonlinear function.

Here is a plot of your original data (black), the standard deviation before x (blue), after x (red) and the sum of the last two (green). I recognize that this question is old. But I wanted to throw another method out there. Canny wrote a paper (A Computational Approach to Edge Detection) where he solved this problem in the two dimensional case for edge detection in images. You can read the paper if you like, but to cut to the chase, you can get a very good approximation to the change point by doing the following:

1. Perform the convolution of your signal with the derivative of the Gaussian.

$$f(x) = \frac{-x}{\chi^2} e^{-\frac{x^2}{\sigma^2}}$$

where $\chi$ is a scaling factor.

1. The peak of the response is where the change point occurs. In my experience with this method, I have found that choosing the correct scaling parameters is difficult. But there may be more work on this that I am not aware of.

I have summarized Canny's paper and provided an example here.

• Hi, I appreciate the response and interesting method presented. Interesting blog post too! The project that generated this question involved some image processing, so I had read a bit about canny edge detection then. – jm22b Mar 24 '17 at 12:56

You may investigate Wavelet tranformed time-series using short period types haar/db4. I have no pointers but just some search terms, try 'wavelet change point detection'.

There are several R packages on Wavelets, see time-series task view: https://cran.r-project.org/web/views/TimeSeries.html

For theory look for Mallat et Hwang paper: "Singularity detection and processing with wavelets"

See related answer: Application of wavelets to time-series-based anomaly detection algorithms