13
$\begingroup$

The Theil-Sen estimator finds the slope and intercept of a line passing through a set of points by calculating the median slope and median intercept of the set of lines passing through all possible distinct point pairs. It is spectacular at fitting a line through data containing outliers.

Here is a reasonably efficient way of calculating this for medium-sized sets of points:

slope[data:{{_,_}..}]:=Median[
 Join@@Table[
  (data[[;;-(n+1),2]] - data[[n+1;;,2]]) / (data[[;;-(n+1),1]] - data[[n+1;;,1]]),
  {n,Length[data]-1}
 ]
]
intercept[data:{{_,_}..}]:=Median[
 Join@@Table[
  (data[[;;-(n+1),1]] data[[n+1;;,2]] - data[[n+1;;,1]] data[[;;-(n+1),2]]) /
  (data[[;;-(n+1),1]] - data[[n+1;;,1]]),
  {n,Length[data]-1}
 ]
]

However, on a machine with 16 GB of RAM, this runs out of memory when processing somewhere between 20,000 points and 50,000 points. How can the code be made more memory-efficient to operate on bigger datasets?

Here is a way of generating an example value for data:

sampleData[nPoints_Integer?Positive]:=(
 SeedRandom[0];
 datax = 10 N[Normalize[Range[nPoints],Max]];
 Transpose[{
  datax,
  1 + 0.001 datax + RandomReal[NormalDistribution[0,0.01],nPoints] +
  datax RandomChoice[{0.01,1-0.01}->{0.02,0},nPoints]
 }]
)
$\endgroup$
6
  • $\begingroup$ Maybe some useful information here. $\endgroup$ – b.gates.you.know.what Jan 31 '15 at 10:44
  • $\begingroup$ TAOCP offers a collection of extremely efficient methods for such tasks. Truth is: Basic principles of computer sciences cannot be outsmarted by any kind of "standard software". $\endgroup$ – Jinxed Jan 31 '15 at 23:04
  • 2
    $\begingroup$ Is there some reason not to use a sampling of point pairs, say n log(n) or so for n points? I think I know of an approach that uses maybe O(n) storage but the time complexity would be n^2 log(n) and that's quite steep for the size range in question. And fishing through the details for an actual implementation would be not so easy. $\endgroup$ – Daniel Lichtblau Feb 4 '15 at 16:55
  • $\begingroup$ @DanielLichtblau, is a random sample guaranteed to give the same answer? One of the main benefits of the technique is its insensitivity to outliers, even if a substantial fraction of the data is composed of outliers. I'd hate to get unlucky and randomly sample only pairs that include the outlier data. $\endgroup$ – ArgentoSapiens Feb 4 '15 at 18:26
  • $\begingroup$ No, unfortunately there is no such guarantee beyond what statistics might indicate (e.g. in terms of variance from "correct" median). $\endgroup$ – Daniel Lichtblau Feb 4 '15 at 19:36
5
$\begingroup$

As a starting point, here is a compact implementation of Theil-Sen:

theilSen[data_?MatrixQ] := Median[With[{df = First[Differences[#]]},
                             {df[[2]], -Det[#]}/df[[1]]] & /@ Subsets[data, {2}]]

It is a bit more efficient than the routine in the OP when I tested it on small sets of points.


In fact, this formulation shows why it may have trouble "when processing somewhere between 20,000 points and 50,000 points." For data of length $n$, you will be generating $\binom{n}{2}$ subsets. $\binom{50\,000}{2}\approx1.25\times 10^9$, and having that many point pairs will indeed give your computer a hard time. Thus, to process that many points, one might consider resorting to lazy subset generation. But, this presents another problem: updating the median when a slope and intercept are generated. There has been some work (e.g. this) on median updating, but I haven't gotten around to fully digesting the literature. I might edit this post later if I figure out a nice implementation of median updating.

$\endgroup$
2
  • $\begingroup$ (+1) It is amazing what you packed into a single line of code! $\endgroup$ – JimB Mar 25 '18 at 1:50
  • $\begingroup$ I was amused at how compact it was, but seeing it break on moderately-sized datasets was a cold shower. ;) $\endgroup$ – J. M.'s torpor Mar 25 '18 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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