15
$\begingroup$

Are there any documented algorithms to separate sections of a given dataset into different curves of best fit?

For example, most humans looking at this chart of data would readily divide it into 3 parts: a sinusoidal segment, a linear segment, and the inverse exponential segment. In fact, I made this particular one with a sine wave, a line and a simple exponential formula.

Chart of data with three distinct parts visible

Are there existing algorithms for finding parts like that, which can then be separately fitted to various curves/lines to make a kind of compound series of best-fits of subsets of the data?

Note that although the example has the ends of the segments pretty much line up, this won't necessarily be the case; there may also be a sudden jolt in the values at a segment cutoff. Perhaps those cases will be easier to detect.

Update: Here is an image of a small bit of real-world data: Real world chart

Update 2: here is an unusually small real-world set of data (only 509 data points):

4,53,53,53,53,58,56,52,49,52,56,51,44,39,39,39,37,33,27,21,18,12,19,30,45,66,92,118,135,148,153,160,168,174,181,187,191,190,191,192,194,194,194,193,193,201,200,199,199,199,197,193,190,187,176,162,157,154,144,126,110,87,74,57,46,44,51,60,65,66,90,106,99,87,84,85,83,91,95,99,101,102,102,103,105,110,107,108,135,171,171,141,120,78,42,44,52,54,103,128,82,103,46,27,73,123,125,77,24,30,27,36,42,49,32,55,20,16,21,31,78,140,116,99,58,139,70,22,44,7,48,32,18,16,25,16,17,35,29,11,13,8,8,18,14,0,10,18,2,1,4,0,61,87,91,2,0,2,9,40,21,2,14,5,9,49,116,100,114,115,62,41,119,191,190,164,156,109,37,15,0,5,1,0,0,2,4,2,0,48,129,168,112,98,95,119,125,191,241,209,229,230,231,246,249,240,99,32,0,0,2,13,28,39,15,15,19,31,47,61,92,91,99,108,114,118,121,125,129,129,125,125,131,135,138,142,147,141,149,153,152,153,159,161,158,158,162,167,171,173,174,176,178,184,190,190,185,190,200,199,189,196,197,197,196,199,200,195,187,191,192,190,186,184,184,179,173,171,170,164,156,155,156,151,141,141,139,143,143,140,146,145,130,126,127,127,125,122,122,127,131,134,140,150,160,166,175,192,208,243,251,255,255,255,249,221,190,181,181,181,181,179,173,165,159,153,162,169,165,154,144,142,145,136,134,131,130,128,124,119,115,103,78,54,40,25,8,2,7,12,25,13,22,15,33,34,57,71,48,16,1,2,0,2,21,112,174,191,190,152,153,161,159,153,71,16,28,3,4,0,14,26,30,26,15,12,19,21,18,53,89,125,139,140,142,141,135,136,140,159,170,173,176,184,180,170,167,168,170,167,161,163,170,164,161,160,163,163,160,160,163,169,166,161,156,155,156,158,160,150,149,149,151,154,156,156,156,151,149,150,153,154,151,146,144,149,150,151,152,151,150,148,147,144,141,137,133,130,128,128,128,136,143,159,180,196,205,212,218,222,225,227,227,225,223,222,222,221,220,220,220,220,221,222,223,221,223,225,226,227,228,232,235,234,236,238,240,241,240,239,237,238,240,240,237,236,239,238,235

Here it is, charted, with the appoximate position of some known real-world element edges marked with dotted lines, a luxury we won't normally have:

enter image description here

One luxury we do have, however, is hindsight: the data in my case is not a time series, but is rather spatially related; it only makes sense to analyse a whole dataset (usually 5000 - 15000 data points) at once, not in an ongoing manner.

Improve this question
$\endgroup$
12
  • $\begingroup$ One possibility would be to fit the whole family of curves at once, using a meta-model. To make things more precise, suppose your ultimate objective is to smooth that histogram, say using a KDE. Then, your smooth estimation from the KDE will be more precise if you use a model in which the width of the kernel is allowed to vary over the range of values of $x$ as in the model used here, equations (2)-(3) $\endgroup$
    – user603
    Commented Jun 12, 2015 at 10:51
  • 1
    $\begingroup$ You constructed the example so that the idea makes sense: so far, so good. With real histograms, it is much more common that a complicated shape reflects a mixture of overlapping distributions: the interest is not then in changepoints on the observed histogram which generally don't exist convincingly or are not the right way to think about mixtures. It's possible, however, that you are using "histogram" in a much broader way than is standard in statistical science where it means bar chart of frequency or probability distribution (only). $\endgroup$
    – Nick Cox
    Commented Jun 12, 2015 at 11:36
  • $\begingroup$ @IrishStat - the usual datasets have 5000 to 15000 entries. I was trying to prepare a summarised real one for here, but it turned out to be a bad example, and I had to start over. On the other hand, doing that did suggest a partial answer to me in terms of simply smoothing and averaging clumps of data for looking initially for patterns, to be finessed later so thanks for that :) I have a real one thats only 509 wide that looks like it could be good; I'll add that to the question when I can. $\endgroup$
    – whybird
    Commented Jun 12, 2015 at 12:31
  • $\begingroup$ @user603 - thank you, that looks very useful, especially section 2.2, and I will be reading it in detail in the morning (I am no doubt in a very different timezone than a statistically significant number of you on here :) ) $\endgroup$
    – whybird
    Commented Jun 12, 2015 at 12:34
  • $\begingroup$ @NickCox - Thanks, yes, I used the word histogram incorrectly. I'll edit the question appropriately. I'm not a statistician, and it is entirely likely that I am inadvertently using other terms incorrectly as well. What I have is a stream of measurements with values which will be 0-255 in practice in my specific example. Patterns observable in the data, in order, will represent real-world phenomena. $\endgroup$
    – whybird
    Commented Jun 12, 2015 at 12:40

3 Answers 3

Reset to default
2
+50
$\begingroup$

My interpretation of the question is that the OP is looking for methodologies that would fit the shape(s) of the examples provided, not the HAC residuals. In addition, automated routines that don't require significant human or analyst intervention are desired. Box-Jenkins may not be appropriate, despite their emphasis in this thread, since they do require substantial analyst involvement.

R modules exist for this type of non-moment based, pattern matching. Permutation distribution clustering is such a pattern matching technique developed by a Max Planck Institute scientist that meets the criteria you've outlined. Its application is to time series data, but it's not limited to that. Here's a citation for the R module that's been developed:

pdc: An R Package for Complexity-Based Clustering of Time Series by Andreas Brandmaier

In addition to PDC, there's the machine learning, iSax routine developed by Eamon Keogh at UC Irvine that's also worth comparison.

Finally, there's this paper on Data Smashing: Uncovering Lurking Order in Data by Chattopadhyay and Lipson. Beyond the clever title, there is a serious purpose at work. Here's the abstract: "From automatic speech recognition to discovering unusual stars, underlying almost all automated discovery tasks is the ability to compare and contrast data streams with each other, to identify connections and spot outliers. Despite the prevalence of data, however, automated methods are not keeping pace. A key bottleneck is that most data comparison algorithms today rely on a human expert to specify what ‘features’ of the data are relevant for comparison. Here, we propose a new principle for estimating the similarity between the sources of arbitrary data streams, using neither domain knowledge nor learning. We demonstrate the application of this principle to the analysis of data from a number of real-world challenging problems, including the disambiguation of electro-encephalograph patterns pertaining to epileptic seizures, detection of anomalous cardiac activity fromheart sound recordings and classification of astronomical objects from raw photometry. In all these cases and without access to any domain knowledge, we demonstrate performance on a par with the accuracy achieved by specialized algorithms and heuristics devised by domain experts. We suggest that data smashing principles may open the door to understanding increasingly complex observations, especially when experts do not know what to look for."

This approach goes way beyond curvilinear fit. It's worth checking out.

Improve this answer
$\endgroup$
13
  • $\begingroup$ Thank you - you are correct that what I want is to find clusters automatically, without analyst intervention. For what I am wanting to do to work, I will need to break datasets of 5000-15000 data points into clusters that each conform well to simple formulae (including repetitive ones) without human intervention over groups of around 50000 such datasets in a timeframe tolerable by humans on domestic computer hardware. $\endgroup$
    – whybird
    Commented Jun 15, 2015 at 23:22
  • $\begingroup$ As for which curve to fit to each cluster, once I have detected the boundaries by whatever means, it is simple enough I think to just try different models (sine wave, polynomial, exponential) and see which gives a better ordinary r^2. $\endgroup$
    – whybird
    Commented Jun 15, 2015 at 23:26
  • 2
    $\begingroup$ OK, I think the miscommunication arises from this: Sax and iSax are representations formats for storing and processing time series, they are not clustering or segment/pattern detection algorithms (per OP's post). My understanding from your answer was that Keogh had come up with an algorithm that is based on the SAX representation format and happens to address the OP's problem. But I think this is not what you meant? $\endgroup$
    – Zhubarb
    Commented Jun 17, 2015 at 14:12
  • 2
    $\begingroup$ OK, no need to reach out to Keogh, I know about iSax and Sax, they are representation formats for efficient mining of time series. The links explain them. iSax is the newer version. I was excited by my misunderstanding of your answer, hence the questions (not trying to be pedantic) :) . $\endgroup$
    – Zhubarb
    Commented Jun 17, 2015 at 14:24
  • 2
    $\begingroup$ i wasn't trying to conceal anything, i interpreted 'isax routine' as an algorithm operating on isax. I suggest your answer needs re-wording / modification after the clarification. $\endgroup$
    – Zhubarb
    Commented Jun 17, 2015 at 14:57
3
$\begingroup$

Detecting change points in a time series requires the construction of a robust global ARIMA model (certainly flawed by model changes and parameter changes over time in your case ) and then identifying the most significant change point in the parameters of that model. Using your 509 values the most significant change point was around period 353. I used some proprietary algorithms available in AUTOBOX (which I have helped develop) that could possibly be licensed for your customized application. The basic idea is to separate the data into two parts and upon finding the most important change point re-analyze each of the two time ranges separately (1-352 ; 353-509 ) to determine further change points within each of the two sets. This is repeated until you have k subsets. I have attached the first step using this approach. Visually it appears to me that the most important point change point has been identified. enter image description here

enter image description here

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ Why does 353 get flagged when 153 and 173 have lower P-values? $\endgroup$
    – Nick Cox
    Commented Jun 13, 2015 at 11:55
  • 1
    $\begingroup$ @NickCox Good question ! Great comment For forecasting purposes the whole idea is to separate the most recent (significant) subset from the older subset which is why 353 won .... For purposes here one would indeed select 173 . $\endgroup$
    – IrishStat
    Commented Jun 13, 2015 at 13:33
  • $\begingroup$ The title "MOST RECENT SIGNIFICANT BREAK POINT" attempts to tell the story $\endgroup$
    – IrishStat
    Commented Jun 13, 2015 at 13:41
  • 1
    $\begingroup$ Thanks for the explanation: the idea is indeed explicit in the last note. (incidentally, I haven't seen so much UPPER CASE in program output since about the early 1990s. I would recommend changing "95% confidence level" to "5% significance level" assuming that is what is meant.) $\endgroup$
    – Nick Cox
    Commented Jun 14, 2015 at 13:37
  • 1
    $\begingroup$ Your wit is matched by your scholarship. Thanks for both . $\endgroup$
    – IrishStat
    Commented Jun 14, 2015 at 13:41
2
$\begingroup$

I think that the title of the thread is misleading: You are not looking to compare density functions but you are actually looking for structural breaks in a time series. However, you do not specify whether these structural breaks are supposed to be found in a rolling time window or in hindsight by looking at the total history of the time series. In this sense your question is actually a duplicate to this: What method to detect structural breaks on time series?

As mentioned by Rob Hyndman in this link, R offers the strucchange package for this purpose. I played around with your data but I must say that the results are dissappointing [is the first data point really 4 or supposed to be 54?]: 

raw = c(54,53,53,53,53,58,56,52,49,52,56,51,44,39,39,39,37,33,27,21,18,12,19,30,45,66,92,118,135,148,153,160,168,174,181,187,191,190,191,192,194,194,194,193,193,201,200,199,199,199,197,193,190,187,176,162,157,154,144,126,110,87,74,57,46,44,51,60,65,66,90,106,99,87,84,85,83,91,95,99,101,102,102,103,105,110,107,108,135,171,171,141,120,78,42,44,52,54,103,128,82,103,46,27,73,123,125,77,24,30,27,36,42,49,32,55,20,16,21,31,78,140,116,99,58,139,70,22,44,7,48,32,18,16,25,16,17,35,29,11,13,8,8,18,14,0,10,18,2,1,4,0,61,87,91,2,0,2,9,40,21,2,14,5,9,49,116,100,114,115,62,41,119,191,190,164,156,109,37,15,0,5,1,0,0,2,4,2,0,48,129,168,112,98,95,119,125,191,241,209,229,230,231,246,249,240,99,32,0,0,2,13,28,39,15,15,19,31,47,61,92,91,99,108,114,118,121,125,129,129,125,125,131,135,138,142,147,141,149,153,152,153,159,161,158,158,162,167,171,173,174,176,178,184,190,190,185,190,200,199,189,196,197,197,196,199,200,195,187,191,192,190,186,184,184,179,173,171,170,164,156,155,156,151,141,141,139,143,143,140,146,145,130,126,127,127,125,122,122,127,131,134,140,150,160,166,175,192,208,243,251,255,255,255,249,221,190,181,181,181,181,179,173,165,159,153,162,169,165,154,144,142,145,136,134,131,130,128,124,119,115,103,78,54,40,25,8,2,7,12,25,13,22,15,33,34,57,71,48,16,1,2,0,2,21,112,174,191,190,152,153,161,159,153,71,16,28,3,4,0,14,26,30,26,15,12,19,21,18,53,89,125,139,140,142,141,135,136,140,159,170,173,176,184,180,170,167,168,170,167,161,163,170,164,161,160,163,163,160,160,163,169,166,161,156,155,156,158,160,150,149,149,151,154,156,156,156,151,149,150,153,154,151,146,144,149,150,151,152,151,150,148,147,144,141,137,133,130,128,128,128,136,143,159,180,196,205,212,218,222,225,227,227,225,223,222,222,221,220,220,220,220,221,222,223,221,223,225,226,227,228,232,235,234,236,238,240,241,240,239,237,238,240,240,237,236,239,238,235)
raw = log(raw+1)
d = as.ts(raw,frequency = 12)
dd = ts.intersect(d = d, d1 = lag(d, -1),d2 = lag(d, -2),d3 = lag(d, -3),d4 = lag(d, -4),d5 = lag(d, -5),d6 = lag(d, -6),d7 = lag(d, -7),d8 = lag(d, -8),d9 = lag(d, -9),d10 = lag(d, -10),d11 = lag(d, -11),d12 = lag(d, -12))

(breakpoints(d ~d1 + d2+ d3+ d4+ d5+ d6+ d7+ d8+ d9+ d10+ d11+ d12, data = dd))
>Breakpoints at observation number:
>151 
>Corresponding to breakdates:
>163 

(breakpoints(d ~d1 + d2, data = dd))
>Breakpoints at observation number:
>95 178 
>Corresponding to breakdates:
>107 190

I am not a regular user of the package. As you can see it depends on the model you fit on the data. You can experiment with 

library(forecast)
auto.arima(raw)

which gives you the best fitting ARIMA model.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you! I have edited out the word 'histogram' from the title; I had usud it incorrectly initially, and forgot to edit the title when I removed it from the body in an earlier edit in response to a comment. $\endgroup$
    – whybird
    Commented Jun 14, 2015 at 0:06
  • $\begingroup$ My data is actually a series of spatially related data, it is not time-based and will usually not exist on a straight line or even in a plane often enough - but you are right that at some fundamental level it can be considered in the same way; I guess that may be part of why my earlier searches didn't find the answers I was expecting. $\endgroup$
    – whybird
    Commented Jun 14, 2015 at 0:08
  • $\begingroup$ The first data point in that example is really a 4, but it could well be that we happened to hit the end of a previous structure or perhaps it was noise; I'd be happy to leave it out as an outlier, but whatever system I come up with will have to cope with things like that too. $\endgroup$
    – whybird
    Commented Jun 14, 2015 at 0:10
  • $\begingroup$ Oh, and the analysis is in hindsight. I'll edit the question to clarify. $\endgroup$
    – whybird
    Commented Jun 14, 2015 at 0:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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