8
$\begingroup$

Suppose I have a dataset:

testdata = 
SortBy[#, #[[1]] &] &@
Flatten[{{#, #/(1 + Exp[#/20] // Log) + RandomReal[0.3]} & /@ 
  Range[100], {#, 10 + 20/#^0.3 + RandomReal[0.3]} & /@ 
  Range[100]}, 1];

ListPlot[testdata] looks like:

enter image description here

We can see it contains "two curves". If I have the testdata, then how can I separate these two curves? How can I extract them from the mixed dataset?

The dataset above is a simple example. In reality, the mixed dataset comes from the solution of a complicated dispersion equation which is hard to factorization.

$\endgroup$
6
  • $\begingroup$ I am wondering, if the actual data points are equally dense and have similar variance per curve. Is it a problem of isolating a linear and a quadratic component or are there perhaps points that can profitably use other fitting functions? Here is the result of my first attempt but it depends on the information I am seeking through this comment. $\endgroup$
    – Syed
    Oct 15, 2021 at 13:38
  • $\begingroup$ If the points of your 2 curves are adjacent to each other like in the example, you may write a simple Table command to extract the curves. E.g. first curve: Table[testdata[[i]], {i, 1, Length[testdata], 2}] and second curve: Table[testdata[[i]], {i, 2, Length[testdata], 2}] $\endgroup$ Oct 15, 2021 at 14:06
  • $\begingroup$ @DanielHuber No..that's cheating...the real testdata I have to deal with is unordered and contains more than two "curves". $\endgroup$
    – Harry
    Oct 15, 2021 at 16:57
  • $\begingroup$ @Syed That's really nice... The actual data points are equally dense and data error is very small. however, I don't have prior knowledge about what fitting functions should be used (Though we can guess it when I plot the dataset out). "linear and a quadratic" is just an oversimplified example. $\endgroup$
    – Harry
    Oct 15, 2021 at 17:07
  • $\begingroup$ To extract near-collinear points, it is important to have a linear component in the data. If there are oscillations or interference from other polynomials or saturation artifacts, things will only go downhill from there. Is there any chance you can load a part/picture of the data so that readers can gauge their options? $\endgroup$
    – Syed
    Oct 15, 2021 at 17:19

1 Answer 1

9
$\begingroup$

This is a very crude, inefficient and not very Mathematica-like answer, meant only to provide an algorithmic example of approaching this problem.

This algorithm can be thought of as some kind of an "Euler-like integrator".

  1. Pick a random point $r_0$ from the convex hull
  2. Find the nearest point $r_1$
  3. Approximate the derivative $\Delta r \approx r_1 - r_0$
  4. Estimate the next point by $r_3^* = r_0 + 2\Delta r$
  5. Find $\tilde{r}$ as the nearest point to the $r_3^*$:
    • If this point is close enough ($|\tilde{r}-r_3^*|<\varepsilon$), let $\tilde{r}$ become $r_3$. Go to 3.
    • Otherwise, complete the curve. Go to 1.

The interpolation could be also of higher order to provide a better estimate $r_3^*$.

Note 1: The algorithm doesn't complete the curve immediately, but first tries to continue it at it other end.

Note 2: There is some randomness to it, and some runs may give worse results.

Note 3: Up to 2 points can miss in the curves.

findNewStartingPoint[data_ /; Length@data <= 2] = {0, 0};
findNewStartingPoint[data_] := Module[{hull, first, second, newdata},
   hull = MeshCoordinates[ConvexHullMesh[data]];
   first = First@RandomSample[hull, 1];
   newdata = Delete[data, Position[data, first]];
   second = First@Nearest[newdata -> {"Index", "Element"}, first];
   newdata = Delete[newdata, second[[1]]];
   {{first, second[[2]]}, newdata}
   ];

findCurves[pts0_, eps_] := 
  Module[{curves, currCurve, pts, pnt, direction, dr},
   direction = 1;
   curves = {};
   {currCurve, pts} = findNewStartingPoint[pts0];
   
   While[Length@pts > 0,
    dr = currCurve[[-1]] - currCurve[[-2]];
    pnt = First@Nearest[pts -> {"Index", "Element", "Distance"},
       currCurve[[-2]] + 2 dr
       ];
    
    If[(pnt[[3]] > eps && direction == -1) || Length@pts == 1,
     (* Curve completed, start a new one *)
     AppendTo[curves, currCurve];
     {currCurve, pts} = findNewStartingPoint[pts];
     direction = 1
     ,
     If[pnt[[3]] > eps && direction == 1,
      (* Continue at the other end *)
      direction = -1;
      currCurve = Reverse@currCurve;
      ,
      (* Add point to current curve *)
      AppendTo[currCurve, pnt[[2]]];
      pts = Delete[pts, pnt[[1]]];
      ]]];
   curves
   ];

curves = findCurves[N@testdata, 3];
ListPlot[curves]

Mathematica graphics

testdata2 = 
  N@RandomSample[#, Length@#] &@
   Catenate[
    Table[{i, # + RandomReal[.5]}, {i, 0, 100}] & /@ {i - 20, i^2/200,
       200/(i + 5), -i^2/180 + 50, 70 - i, 10 Sin[i/5] + 45}];
ListPlot[findCurves[N@testdata2, 3], Axes -> False, 
    FrameTicks -> False, Frame -> True, AspectRatio -> 1] & /@ Range[3] // Row

Mathematica graphics

$\endgroup$
3
  • $\begingroup$ This is really nice already....thank you very much! $\endgroup$
    – Harry
    Oct 16, 2021 at 6:15
  • $\begingroup$ @Domen, in step 5 of the algorithm, "If this point is close enough ($|\tilde{r}-r_3^*|<\varepsilon$), let $\tilde{r}$ become $r_3$", should it be "let $\tilde{r}$ become $r_1$ and $r_1$ become $r_0$"? $\endgroup$
    – lxy
    Mar 1, 2023 at 14:40
  • $\begingroup$ @Domen, btw, in your algorithm there is no $r_3$ only $r_3^\ast$ $\endgroup$
    – lxy
    Mar 2, 2023 at 1:36

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.