12
$\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$
  • $\begingroup$ Please post your data $\endgroup$ – Conor Cosnett Aug 5 '16 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$ – J. M. will be back soon Aug 5 '16 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 '16 at 11:08
  • $\begingroup$ @J.M. Would be very grateful if you could provide some details in an answer. $\endgroup$ – მამუკა ჯიბლაძე Aug 5 '16 at 11:09
11
$\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$

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$
3
$\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$
  • 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 '16 at 11:18
  • $\begingroup$ They are related. They all are functions of the same independent variable t. $\endgroup$ – Conor Cosnett Aug 5 '16 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 '16 at 11:27
  • 1
    $\begingroup$ ok! I understand now! I will try to find a solution without knowledge of t! $\endgroup$ – Conor Cosnett Aug 5 '16 at 11:30

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.