13
$\begingroup$

I have 3D data - a bunch of triples like {{x1, y1, z1}, {x2, y2, z2}, ...}, and I know they lie on a curve rather than a surface; in fact, I need a least squares fit of these points to a 3D straight line. That is, I am looking for six numbers ax, bx, ay, by, az, bz such that my points are as close as possible to the line {ax*t + bx, ay*t + by, az*t + bz}, with t running through reals. I could not find a way to do it in Mathematica. Does anybody know a way?

$\endgroup$
4
  • $\begingroup$ Please post your data $\endgroup$
    – conor
    Aug 5, 2016 at 10:06
  • $\begingroup$ Orthogonal fitting is pretty much the only viable option in this case, but you'll have to set up FindMinimum[] yourself to do it, since this method is not built-in. $\endgroup$ Aug 5, 2016 at 10:09
  • 1
    $\begingroup$ @JohnConorCosnett I have 1285 points, and unfortunately I don't know how to take a sensible small sample from them, sorry. You could take something like Table[{.5 t + .1, 1.2 t - .3, -2.2 t + 1.2}+RandomReal[{-1,1},3], {t, RandomReal[{-10, 10}, 100]}] $\endgroup$ Aug 5, 2016 at 11:08
  • $\begingroup$ @J.M. Would be very grateful if you could provide some details in an answer. $\endgroup$ Aug 5, 2016 at 11:09

3 Answers 3

13
$\begingroup$

As it turns out, you don't need FindMinimum[] in the linear case of total least squares/orthogonal distance regression; all that is needed is a clever application of SVD:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            p = RandomReal[{-2, 2}, 3]; (* point on true line *)
            (* direction cosines *)
            q = Normalize[RandomVariate[NormalDistribution[], 3]];
            (* random points clustered near the line *)
            pts = Table[p + t q + RandomVariate[NormalDistribution[0, 1/10], 3],
                        {t, 0, 1, 1/90}];]

(* orthogonal fit *)
lin = InfiniteLine[Mean[pts], Flatten[Last[
                   SingularValueDecomposition[Standardize[pts, Mean, 1 &], 1]]]];

Legended[Graphics3D[{{Directive[AbsolutePointSize[6], Brown], Point[pts]},
                     {Directive[AbsoluteThickness[4], ColorData[97, 1]], 
                      lin},
                     {Directive[AbsoluteThickness[4], ColorData[97, 3]], 
                      InfiniteLine[p, q]}}, Axes -> True], 
         LineLegend[{ColorData[97, 1], ColorData[97, 3]},
                    {"orthogonal fit", "true line"}]]

orthogonally-fitted line

$\endgroup$
4
  • $\begingroup$ Is it easy to see that the best linear fit line must go through the mass center? $\endgroup$ Aug 5, 2016 at 12:24
  • $\begingroup$ Also, it might be quite nontrivial to estimate errors for the comparison of various possible approaches... $\endgroup$ Aug 5, 2016 at 12:37
  • 3
    $\begingroup$ Here is a nice survey paper you might want to look at. If you want the residuals, you can use the usual method for finding the distance of a point to a line. $\endgroup$ Aug 5, 2016 at 12:54
  • $\begingroup$ Thanks a lot, very useful link. $\endgroup$ Aug 5, 2016 at 13:55
4
$\begingroup$

Inspired by the answer by John Conor Cosnett I came up with something:

data = Table[{.5, 1.2, -2.2} t + {.1, -.3, 1.2} + RandomReal[{-1, 1}, 3],  
 {t, RandomReal[{-10, 10}, 100]}];
xyfit = FindFit[data[[All, {1, 2}]], axy x + bxy, {axy, bxy}, x]
xzfit = FindFit[data[[All, {1, 3}]], axz x + bxz, {axz, bxz}, x]

Seems to give good results,

Show[
 ParametricPlot3D[{x, axy x + bxy /. xyfit, axz x + bxz /. xzfit},
  {x, -10000, 10000}, PlotStyle -> Red],
 ListPointPlot3D[data],
 BoxRatios -> {1, 1, 1}, PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10}}]

produces this:

enter image description here

Still I am in doubt since this approach somehow breaks symmetry, as it treats $x$ as an independent variable, with $y$ and $z$ as its functions. I somehow suspect this might introduce some bias, with errors not uniform wrt the variables. So I am leaving this unaccepted, maybe somebody can come up with something better.

$\endgroup$
2
$\begingroup$

Complete code for finding {ax, bx, ay, by, az, by} using your example data from your comment:

t = RandomReal[{-10, 10}, 100];

points = pts = 
 Table[{.5 t + .1 + RandomReal[], 
 1.2 t - .3 + RandomReal[], -2.2 t + 1.2 + RandomReal[]}, {t, 
 RandomReal[{-10, 10}, 100]}];

 x = #[[1]] & /@ points; 
 y = #[[2]] & /@ points;
 z = #[[3]] & /@ points;


 Xdata = Thread[{t, x}];
 Ydata = Thread[{t, y}];
 Zdata = Thread[{t, z}];

 Clear[t]
  Join[
   FindFit[Xdata, ax*t + bx, {ax, bx}, t],
   FindFit[Ydata, ay*t + by, {ay, by}, t],
   FindFit[Zdata, az*t + bz, {az, bz}, t]
  ]

Contrived ExampleData:

f[t_] := {1 t + 2, 3*t + 4, 5*t + 6}

points = f /@ Range[0, 10, 0.1]

Split up data into 3 linear regressions:

$x(t)=ax* t+ bx$

$y(t)=ay* t+ by$

$z(t)=az* t+ bz$

Make up $t$ data.

t = Range[0, 10, 0.1]

Extract x, y, z:

x = #[[1]] & /@ points;
y = #[[2]] & /@ points;
z = #[[3]] & /@ points;

Combine into {input, output} lists:

Xdata = Thread[{t, x}];
Ydata = Thread[{t, y}];
Zdata = Thread[{t, z}];

Use 3 separate FindFit operations:

Clear[t]
Join[
  FindFit[Xdata, ax*t + bx, {ax, bx}, t],
  FindFit[Ydata, ay*t + by, {ay, by}, t],
  FindFit[Zdata, az*t + bz, {az, bz}, t]
]
$\endgroup$
4
  • 1
    $\begingroup$ Hmmmm are you sure this is the right way? I mean, it somehow treats $x$, $y$ and $z$ coordinates as entirely unrelated separate random variables, is it OK? $\endgroup$ Aug 5, 2016 at 11:18
  • $\begingroup$ They are related. They all are functions of the same independent variable t. $\endgroup$
    – conor
    Aug 5, 2016 at 11:23
  • $\begingroup$ That's my trouble precisely, I just failed to formulate it well enough. What I have is a bunch of triples, I don't really know the corresponding ts for them. In your solution you use this information, but it must be presupposed to be hidden. $\endgroup$ Aug 5, 2016 at 11:27
  • 1
    $\begingroup$ ok! I understand now! I will try to find a solution without knowledge of t! $\endgroup$
    – conor
    Aug 5, 2016 at 11:30

Your Answer

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

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