# How to split data into clusters based on fitting to function

I have data which looks somewhat like this:

On the picture, you can see that the data can be described by 2 linear functions - if you manage to split it into 2 groups, each can easily be fitted with a x + b

However, when I run FindClusters on that data, I get this:

It seems, I need to adjust it somehow

Can anyone show me how to properly use FindClusters to split this data?

UPDATE

The functions are all known and the same.

To play with, the data can be emulated with:

dd = Join[Table[ {x, .3 x + 2 + RandomReal[{-1, 1}]}, {x, 0, 100}],
Table[ {x, .6 x - 20 + RandomReal[{-1, 1}]}, {x, 0,
100}]]; ListPlot[dd]; cc = FindClusters[dd]; ListPlot[cc]

• Different clustering metrics can "discover" differently shaped clusters. For non-spherical clusters you might consider looking at science.sciencemag.org/content/344/6191/1492. And remember that clustering algorithms are (usually) based solely on the data and no outside information and have a more limited capability to recognize clusters than one's brain can do. So if your eye/brain system sees two lines, then you should fit two lines and not rely on a canned clustering method to do so.
– JimB
May 15 '17 at 16:58
• @JimBaldwin so how would one select the data to fit any of the lines (I mean, without manual intervention)? May 15 '17 at 17:28
• You should ask such a question on CrossValidated. It's important to lay out what exactly you know about the clusters. For example, do you know if the clusters are straight lines? Might they be quadratics? Might they be more like meandering rivers? Are there linear features and spherical clusters? The answers below all assume a specific functional form. You need to give what you know.
– JimB
May 15 '17 at 21:56

If we can assume we know more or less the functions for the two groups, in this case two linear functions (two parameters each) we can just NMinimize for the Min distance of the two functions.

{w, sol} = NMinimize[
Total[
Map[
Min[Abs[#[[2]] - (n1 + m1 #[[1]])],
Abs[#[[2]] - (n2 + m2 #[[1]])]] &
, dd
]
], {n1, n2, m1, m2}]

{95.3241, {n1 -> 1.84575, n2 -> -19.8623, m1 -> 0.300732,  m2 -> 0.598053}}


Once we have extracted the parameters for the two functions, then we just GatherBy each point, based on which function is closer.

ListPlot@Module[{n1, n2, m1, m2},
{n1, n2, m1, m2} = Values[sol];
GatherBy[dd,
LessEqual[Abs[#[[2]] - (n1 + m1 #[[1]])],
Abs[#[[2]] - (n2 + m2 #[[1]])]] &]
]


## More general

Define your function with arbitrary number of parameters

func[parms_List][x_] := parms[[1]] + parms[[2]]  x


A distance function

dist[f_][{x_, y_}] := Abs[y - f[x]]


Or (after comment by @SampoSmolander)

dist[f_][{x_, y_}] := (y - f[x])^2


Fake data

dd = Table[
With[
{
a = RandomChoice[{1, 4}],
b = RandomChoice[{1/2, 2}]
}, {x, func[{a, b}][x] + RandomReal[{-1, 1}]}], {x, 0, 10,
0.01}];


Minimisation, The i index goes through the number of parameters of each function (here 2), the index j over the number of functions to fit (in this case 4).

sol2 = Last@With[{n = 2, m = 4},
NMinimize[
Total@Map[Function[{L}, Min @@ Table[
dist[func[Array[c[#, k] &, n]]][L]
, {k, m}]], dd]
, Flatten[Table[c[i, j], {i, n}, {j, m}]]]
]


Plot

Show[
ListPlot@GatherBy[
dd
, Position[#, Min[#]] &[
Table[dist[func[{c[1, k], c[2, k]} /. sol2]][#], {k, 4}]] &
],
Plot[
Evaluate[
Table[
func[{c[1, k], c[2, k]}][x]
, {k, 4}] /. sol2
]
, {x, 0, 10}
]]


Obviously there is no way to be sure this will converge to reasonable clusters if the lines are too mixed or the functions too complicated.

## Fitting different number of functions

Let's see how does this perform when there are four groups of points and we ask for m different number of clusters, from 1 to 6.

plots = Table[
sol2 = Last@With[{n = 2},
NMinimize[
Total@Map[Function[{L}, Min @@ Table[
dist[func[Array[c[#, k] &, n]]][L]
, {k, m}]], dd]
, Flatten[Table[c[i, j], {i, n}, {j, m}]]]
];
ListPlot[
GatherBy[
dd
, Position[#, Min[#]] &[
Table[dist[func[{c[1, k], c[2, k]} /. sol2]][#], {k, m}]] &
], PlotLabel -> m], {m, 6}];

Grid@Partition[plots, 2]


• I've played around with the code, and there is one thing missing. Let's say function 1's trace contains 100 points, while function 2's trace contains 10 points. In that case, the code provided miss-classifies the points : screencast.com/t/fhvZ3CRT2ow May 15 '17 at 21:21
• It's somewhat more robust if you use quadratic cost function instead of Abs. May 16 '17 at 3:09
• @SampoSmolander By principle, or do you have a benchmark? What about higher powers? I would be interested in understanding the reasons behind. Thanks for your comment. May 16 '17 at 10:27
• No, I just played around a bit, and the quadratic seems to be more robust. But if you make tricky enough data, it can be fooled, too. May 16 '17 at 17:28
• @ArsenZahray I think it's not fair to ask now to fine tune the answer to specific cases ("1's trace contains 100 points, while function 2's trace contains 10 points"). You already moved the goalpost once by asking for more functions. Now you have to do your part. We are trying to provide general help, not a free coding service. Please read this meta post May 17 '17 at 7:12

Another approach to this problem is to use the image processing capabilities. This code creates a binary image from the data, then uses ImageLines, which implements the Radon transform to locate the most prominent lines. These are then highlighted in orange.

dd = Join[Table[{x, .3 x + 2 + RandomReal[{-1, 1}]}, {x, 0, 100}],
Table[{x, .6 x - 20 + RandomReal[{-1, 1}]}, {x, 0, 100}]];
img = ListPlot[dd, Axes -> False] // Image;
binImg = Binarize[ColorNegate@img];
lines = ImageLines[binImg, 0.2];
HighlightImage[binImg, {Orange, Line /@ lines}]


In the event of more lines, you can change the second argument to ImageLines.

This does the same as the other answer (fit two lines to the data, cluster data points based on whether they are closer to one line or the other), but with a bit different coding style:

dd = Join[Table[ {x, .3 x + 2 + RandomReal[{-1, 1}]}, {x, 0, 100}],
Table[ {x, .6 x - 20 + RandomReal[{-1, 1}]}, {x, 0, 100}]];

cost1[a1_, b1_, a2_, b2_, {x_, y_}] :=
Min[(a1 + b1 x - y)^2, (a2 + b2 x - y)^2]
cost[a1_, b1_, a2_, b2_] :=
Total[Map[cost1[a1, b1, a2, b2, #]&, dd]]

sol = Minimize[cost[a1, b1, a2, b2], {a1, b1, a2, b2}][[2]]

set1q[{x_, y_}] := (a1 + b1 x - y)^2 < (a2 + b2 x - y)^2 /. sol
set1 = Select[dd, set1q];
set2 = Complement[dd, set1];

ListPlot[{set1, set2}]

• Nice, but I still don't know in advance how many traces there can be. It may be possible that the data contains a single trace, or multiple, and I need to separate them May 15 '17 at 18:11
• If you can estimate the typical magnitude of the residuals (difference between data point and the line fitted to the corresponding "trace"), you can just try to fit and increasing number of lines to the data, and when you find the correct number, the residuals should drop to that level. Before you find the correct number, the residuals will be, on average, much larger. And if you fit more than the correct number, it will not decrease the residuals much. That's how you could automate finding the correct number of lines. May 15 '17 at 18:18

If there are no crossing points, Method->"Agglomerate"option and specifying the number of clusters at 2nd argument works great.

dd = Join[Table[{x, .3 x + 2 + RandomReal[{-1, 1}]}, {x, 0, 100}],
Table[{x, .6 x - 20 + RandomReal[{-1, 1}]}, {x, 0, 100}]];
ListPlot[dd]
cc = FindClusters[Select[dd, #[[1]] < 60 &], 2,
Method -> "Agglomerate"];
cc2 = FindClusters[Select[dd, #[[1]] > 85 &], 2,
Method -> "Agglomerate"];
ListPlot[Join[cc, cc2, 2]]