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?
3 Answers
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"}]]
-
$\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
-
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.
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]
]
-
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$– conorAug 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
t
s 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$– conorAug 5, 2016 at 11:30
FindMinimum[]
yourself to do it, since this method is not built-in. $\endgroup$Table[{.5 t + .1, 1.2 t - .3, -2.2 t + 1.2}+RandomReal[{-1,1},3], {t, RandomReal[{-10, 10}, 100]}]
$\endgroup$