8
$\begingroup$

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

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

Improve this question
$\endgroup$
4
  • $\begingroup$ Maybe it would be useful to consider change-point detection techniques. Circular binary segmentation might be just what you need. $\endgroup$
    – Mur1lo
    Aug 15, 2016 at 18:22
  • $\begingroup$ Chen and Liu 1993 describe a method of detecting level shift in a time series. Their algorithm is implemented in the tsoutliers R package. $\endgroup$
    – pionpi_
    Aug 16, 2016 at 13:34
  • 1
    $\begingroup$ 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.) $\endgroup$
    – Richard Hardy
    Aug 16, 2016 at 18:33
  • $\begingroup$ 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/… $\endgroup$
    – francois-xavier sireta
    Oct 4, 2018 at 16:19

5 Answers 5

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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 $\endgroup$
    – jm22b
    Aug 15, 2016 at 14:13
  • 1
    $\begingroup$ I have implemented your idea and it works like a charm, thanks a lot! $\endgroup$
    – jm22b
    Aug 16, 2016 at 9:40
  • $\begingroup$ That is great to hear. Thank you for letting me know. $\endgroup$
    – jason
    Aug 16, 2016 at 21:10
  • $\begingroup$ Can you explain what you mean by front-to-back and back-to-front filters? $\endgroup$
    – thentangler
    Oct 28, 2020 at 17:48
  • 1
    $\begingroup$ You filter your data. Then you take the results and flip them. The beginning becomes the end and the end the beginning. Then you apply your filter again. Then you flip your data back to its original orientation. The cutoff will be slightly sharper then the original design and there will be no delay $\endgroup$
    – jason
    Oct 29, 2020 at 18:30
4
$\begingroup$

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 enter image description here 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.enter image description here enter image description here

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-enter image description herefit-forecast graph for that approach and the identified robust transfer function model . enter image description here with residual plot enter image description here and residual acf here enter image description 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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – IrishStat
    Aug 16, 2016 at 10:23
  • $\begingroup$ 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! $\endgroup$
    – jm22b
    Aug 16, 2016 at 10:30
3
$\begingroup$

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

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ A Pandas implementation of this. First load data with test_df = pd.read_csv('test.csv') then calculate the "before" standard deviation on y series with test_df['std_before'] = test_df['y'].expanding().std() and calculate the "after" standard deviation on the reversed y series with test_df['std_after'] = test_df['y'][::-1].expanding().std()[::-1]; finally sum and plot with test_df['std_before+after'] = test_df['std_after'] + test_df['std_before'] and test_df.plot(x='x', y = 'std_before+after', figsize=(14,4)); $\endgroup$
    – MyCarta
    Jan 12, 2021 at 21:59
2
$\begingroup$

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.

enter image description here

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – jm22b
    Mar 24, 2017 at 12:56
0
$\begingroup$

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

Look here for some example: http://it.mathworks.com/help/wavelet/examples/detecting-discontinuities-and-breakdown-points.html?requestedDomain=www.mathworks.com

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

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.