4
$\begingroup$

I'm trying to display the plane of best fit in the same 3D-box. Currently, I have the following code which spits out a plot of the points in 3D. I want to have a plane of best fit, can Mathematica do this automatically?

t1 = {{10.26`, 43.`, 330.`}, {985.69`, 28.`, 120.`}, {1445.56`, 35.`, 
360.`}, {1188.19`, 33.`, 270.`}, {574.51`, 44.`, 220.`}, {568.95`,
 19.`, 170.`}, {471.81`, 20.`, 70.`}, {537.35`, 22.`, 
210.`}, {514.07`, 21.`, 200.`}, {174.09`, 40.`, 300.`}, {1720.81`,
 32.`, 290.`}, {611.48`, 20.`, 70.`}, {251.19`, 24.`, 
150.`}, {97.97`, 38.`, 190.`}, {406.81`, 24.`, 240.`}, {265.4`, 
25.`, 100.`}, {1323.29`, 35.`, 250.`}, {196.65`, 36.`, 
210.`}, {1326.6`, 27.`, 280.`}, {1380.69`, 33.`, 230.`}};

ListPointPlot3D[t1]
$\endgroup$
  • $\begingroup$ Yeah every ones got the sheet fit, the plane is pretty rare to find $\endgroup$ – Jethro Devøn Nov 13 '14 at 15:51
  • $\begingroup$ Define fit = Fit[t1, {1, x, y}, {x, y}], then use Plot3D. $\endgroup$ – user484 Nov 13 '14 at 16:18
12
$\begingroup$

Here's a way to get the plane of best fit: subtract the centroid of the data, and then plot the plane generated by the first two left singular vectors of the singular value decomposition of the resulting data:

Y = # - Mean /@ # &[t1\[Transpose]]
{U, S, V} = SingularValueDecomposition[Y];
Graphics3D[{InfinitePlane[{{0, 0, 0}, U[[;; , 1]], U[[;; , 2]]}], 
  Point /@ (Y\[Transpose])}, BoxRatios -> {1, 1, 1}]

which produces

enter image description here

The RMS error is the third singular value:

S[[3, 3]]

which produces

25.8873
| improve this answer | |
$\endgroup$
  • $\begingroup$ +1, I think there may be a syntax error as Y is not defined. Copy and paste does not run. Should it be: {U, S, V} = SingularValueDecomposition[Y = # - Mean /@ # &[t1\[Transpose]]]; Graphics3D[{InfinitePlane[{{0, 0, 0}, U[[;; , 1]], U[[;; , 2]]}], Point /@ (Y\[Transpose])}, BoxRatios -> {1, 1, 1}] or alternatively plotting original points by translating plane. $\endgroup$ – ubpdqn Nov 15 '14 at 11:56
  • $\begingroup$ @ubpdqn: Dang, you're right, sorry 'bout that! Fixing it. $\endgroup$ – DumpsterDoofus Nov 15 '14 at 14:40
  • $\begingroup$ Why don't you plot the plane against the actual data t1? $\endgroup$ – Taiki Oct 15 '15 at 18:51
8
$\begingroup$
Show[
 ListPointPlot3D[t1, PlotStyle -> Red, Filling -> Bottom],
 Plot3D[
  Evaluate@Fit[t1, {1, x, y}, {x, y}]
  , {x, Min[t1[[All, 1]]], Max[t1[[All, 1]]]}
  , {y, Min[t1[[All, 2]]], Max[t1[[All, 2]]]}
  ]]

Mathematica graphics

| improve this answer | |
$\endgroup$
  • $\begingroup$ That's very beautiful graph, mathematica is very powerful thank you for your help, this problem is solved :) $\endgroup$ – Jethro Devøn Nov 13 '14 at 17:05
  • 1
    $\begingroup$ Am I correct in thinking that the plane of best fit is also the sum of least squares? $\endgroup$ – Jethro Devøn Nov 13 '14 at 17:06
  • 4
    $\begingroup$ This particular fit minimizes sums of squares of vertical discrepancies. For total least squares (sums of squares of shortest rather than vertical distances to plane), one would have to do this a bit differently. $\endgroup$ – Daniel Lichtblau Nov 13 '14 at 20:59
  • $\begingroup$ really?? uh oh, so is there a more effective function than fit?? $\endgroup$ – Jethro Devøn Nov 13 '14 at 21:22
  • 3
    $\begingroup$ @rhermans: Not quite, the plane of best fit is different from the plane of best vertical fit, although for point clouds which nearly outline a horizontal plane the difference will be fairly small. One way to get the true plane of best fit (and the associated RMS error) is to use principal value decomposition of the data (see my answer for details). $\endgroup$ – DumpsterDoofus Nov 14 '14 at 16:00
3
$\begingroup$

You can now use ResourceFunction["PlaneOfBestFit"] to find the best plane for a set of points. Using the example from the OP,

In[40]:= ResourceFunction["PlaneOfBestFit"]@t1

Out[40]= InfinitePlane[{{702.569, 29.95, 213.}, {703.567, 29.9501, 
   213.055}, {702.623, 29.8805, 212.004}}]

I recently wanted this functionality and my searches led me to DumpsterDoofus's answer as well as this math forum post. I made a standalone function to get the plane from a set of points, and Daniel Lichblau extended it to work with n-dimensional points.

The reason I wanted to do this was to visualize the plane of a ring in a molecule, to visualize sterochemistry:

m = Molecule["adenosine"];
aliphaticRing = 
  Flatten@Values@
    FindMoleculeSubstructure[m, 
     ResourceFunction["MoleculeRingPattern"][{"O", _, _, _, _}]];
points = QuantityMagnitude@m["AtomCoordinates"][[aliphaticRing]];
plane = ResourceFunction["PlaneOfBestFit"][points];
Show[MoleculePlot3D[m], Graphics3D@{Red, Opacity[.7], plane}]

The 3D plane helps to see how the adenine ring is on the opposite side of the ribose ring from the hydroxyl groups. enter image description here

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ You might like this OG reference: Pearson 1901 $\endgroup$ – Chris K Apr 23 at 1:58
  • $\begingroup$ Very nice - I can add a link to that on the reference page. $\endgroup$ – Jason B. Apr 23 at 2:24
  • $\begingroup$ @Chris and Jason, this would be the official link. ;) Jason, I'll also e-mail you later about this. $\endgroup$ – J. M.'s technical difficulties Apr 23 at 2:41

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.