Fitting 3D points to a straight line

Question

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?

Answer 1

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"}]]

Answer 2

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:

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.

Answer 3

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]
]