robjohn is giving me a hand with this, but in case anybody else knows...
I need to do a least-squares regression for linearity on a set of coordinates in 3space. If the dataset is linear, I need to see if it is close to vertical or horizontal. How could I do this?
Many thanks in advance
Joe Stavitsky
On
The core of the orthogonal regression is expressible by a "principal component" analysis as implemented in many statistical programs. You do a "PCA"-rotation on your centered dataset and use the first two "factors" (=first two principal components) as x and y-axes and the z-axis is then the measure of the vertical distance of the datapoints to that plane. The eigenvalue associated to that third component is then equivalent to the sum-of-squares of the vertical deviances and is minimized by the PC-rotation.
Here is a listing of a calculation (using my matrix-calculator "MatMate" and its own language, sorry, don't have Math'ica or Maple...). Also the occuring matrices are documented. The bracketed numbers are only statement numbers in the listing.
[14] set cldezweite=5 // decimal digits for printing
[15] n=8 // use 8 data points
[16] set randomstart=41
[17] udata = randomu(3,n) ' // get 3 columns of n rows random coordinates
------------------------------------------
udata :
x-axis y-axis z-axis
------------------------------------------
28.65775 4.77896 86.38809
8.65334 62.09805 20.01654
88.59632 78.99541 56.30901
33.55478 18.60703 90.15799
65.22586 81.45162 63.61901
87.68102 39.27773 4.50016
45.83511 48.76525 70.40108
67.04584 85.29792 40.20077
------------------------------------------
[18] center= meansp(udata) // determine the center as means columnwise
------------------------------------------
center :
53.15625 52.40900 53.94908
------------------------------------------
[19] abwdata= udata - center // generate centered data
------------------------------------------
abwdata :
-24.49850 -47.63004 32.43901
-44.50292 9.68905 -33.93254
35.44007 26.58641 2.35993
-19.60147 -33.80197 36.20891
12.06961 29.04263 9.66992
34.52477 -13.13127 -49.44892
-7.32114 -3.64375 16.45200
13.88958 32.88893 -13.74831
------------------------------------------
[20] chk = sqsumsp(abwdata)/n // check the square-sums of each column
------------------------------------------
chk :
725.63213 790.34643 814.84356 // we have s² = 814.84 in the z-axis
------------------------------------------
Get rotationmatrix for principal components rotation (columnwise)
[21] t = gettrans(abwdata,"pca")
; rotate centered data and also the center
[22] optabwdata=abwdata*t
[23] optcenter = center*t
------------------------------------------
optabwdata :
x-axis y-axis z-axis
------------------------------------------
-61.37125 -0.48959 -12.42644
2.16570 50.81196 25.28358
33.93348 -28.56924 -0.85030
-52.31193 -9.23432 -4.00819
18.95788 -20.78838 17.06163
38.42755 20.42223 -43.77032
-15.60503 -7.60057 6.02151
35.80360 -4.55208 12.68852
optcenter :
29.71799 -84.50733 21.40434
------------------------------------------
Translate the rotated data to the (rotated) center
[24] optdata = optabwdata + optcenter
------------------------------------------
optdata :
-31.65326 -84.99691 8.97790
31.88369 -33.69537 46.68792
63.65147 -113.07656 20.55405
-22.59394 -93.74165 17.39615
48.67587 -105.29571 38.46597
68.14554 -64.08510 -22.36598
14.11296 -92.10790 27.42586
65.52159 -89.05941 34.09286
------------------------------------------
Check squaresums in the three new axes. The z-axis(3'rd column) "is minimal" and contains a version of the least-squares-error of the related orthogonal regression
[25] chk1=sqsumsp(abwsp(optdata))/n
------------------------------------------
chk1 :
1377.57898 551.41041 401.83274 // we have s² = 401.83 in the z-axis
// this should be the least possible
// using rotations and translations only
------------------------------------------
; list the rotation-matrix for the columnwise rotation
------------------------------------------
t :
0.52648 -0.59715 -0.60517
0.62564 -0.20985 0.75136
-0.57567 -0.77419 0.26311
------------------------------------------
It is then possible to use the inverse(which is simply the transpose) of t for a function for the plane (or more correctly for the line as requested in your original question) .
For the line we had then for coordinates (x,y,z)' in the framework of the given points by a free real parameter x
(x,y,z)' = (x,0,0)* t' + center
= (x,0,0)* ( 0.52648 0.62564 -0.57567) + center
-0.59715 -0.20985 -0.77419)
-0.60517 0.75136 0.26311)
= (x*0.52648 ,x*0.62564, x* -0.57567) + center
For the plane we had then for coordinates (x,y,z)' in the framework of the given points by two free real parameters (x,y)
(x,y,z)' = (x,y,0)* t' + center
= (x,y,0)* ( 0.52648 0.62564 -0.57567) + center
-0.59715 -0.20985 -0.77419)
-0.60517 0.75136 0.26311)
= (x*0.52648 -y*0.59715,x*0.62564-y*0.20985, -x*0.57567-y*0.77419)
+ center
Here the formulae might be scaled by an arbitrary factor<>0 because (x,y) are arbitrary reals and we might factor out any number to make the decimal expansion of the coefficients at x any y in the equation more smooth.
(you can download MatMate and reproduce the computation)
On
I posted this this morning in chat to Joe, but since the question is here, too, I will post it here.
Uniquely, up to the sign of $u$, specify a line as $p+ut$ where $|u|=1$, $p\cdot u=0$, and $t\in\mathbb{R}$. The square of the distance of $x$ from this line is $$ |x-p|^2-((x-p)\cdot u)^2 = |x-p|^2-(x\cdot u)^2\tag{1} $$ Thus, for a given $\{x_k\}$, we wish to minimize $$ \Delta=\frac1n\sum_{k=1}^n\left[|x_k-p|^2-(x_k\cdot u)^2\right]\tag{2} $$ over all $p$ and $u$ so that $$ |u|=1\text{ and }p\cdot u=0\tag{3} $$ Thus, we need that for any variation of $u$ and $p$ that maintain $(3)$, the variation of $(2)$ must be $0$. That is, whenever $$ u\cdot\delta u=0\quad\text{and}\quad p\cdot\delta u+u\cdot\delta p=0\tag{4} $$ we want the variation of $\Delta$ to vanish $$ 0=\sum_{k=1}^n\left[(x_k-p)\cdot\delta p+u\cdot x_kx_k\cdot\delta u\right]\tag{5} $$ Define $$ \bar{x}=\frac1n\sum_{k=1}^nx_k\quad\text{and}\quad X=\frac1n\sum_{k=1}^nx_kx_k^T\tag{6} $$ and $(5)$ becomes $$ 0=(\bar{x}-p)\cdot\delta p+(Xu)\cdot\delta u\tag{7} $$ Considering $(4)$, if we hold $u$ constant and vary $p$, we get that for any $\delta p$ so that $u\cdot\delta p=0$ $$ 0=(\bar{x}-p)\cdot\delta p\tag{8} $$ Because $(8)$ holds for any $\delta p$ which is perpendicular to $u$, we get that $\bar{x}-p\,||\,u$, which leads to $$ \bar{x}-p=(\bar{x}-p)\cdot uu=\bar{x}\cdot uu\tag{9} $$ Combining $(4)$ and $(9)$ yields $$ \bar{x}\cdot\delta u+u\cdot\delta p=0\tag{10} $$ Plugging $(9)$ into $(7)$ and applying $(10)$ yields $$ \begin{align} 0 &=\bar{x}\cdot uu\cdot\delta p+(Xu)\cdot\delta u\\ &=-\bar{x}\cdot u\bar{x}\cdot \delta u+(Xu)\cdot\delta u\\ &=\left(\left(X-\bar{x}\bar{x}^T\right)u\right)\cdot\delta u\tag{11} \end{align} $$ Since $(11)$ holds for any $\delta u$ perpendicular to $u$, we must have that $$ \left(X-\bar{x}\bar{x}^T\right)u\,||\,u\tag{12} $$ That is, $u$ is an eigenvector of $X-\bar{x}\bar{x}^T$.
Plugging $(9)$ into $(2)$ yields $\Delta$ as a function of $u$, $$ \begin{align} \Delta &=\frac1n\sum_{k=1}^n\left[|x_k-p|^2-(x_k\cdot u)^2\right]\\ &=\frac1n\sum_{k=1}^n|(x_k-\bar{x})+\bar{x}\cdot uu|^2-u^TXu\\ &=\frac1n\sum_{k=1}^n\left(|x_k-\bar{x}|^2+2(x_k-\bar{x})\cdot u\bar{x}\cdot u+(\bar{x}\cdot u)^2\right)-u^TXu\\ &=\operatorname{Var}[\{x_k\}]-u^T\left(X-\bar{x}\bar{x}^T\right)u\tag{13} \end{align} $$
and since $$ X-\bar{x}\bar{x}^T=\frac1n\sum_{k=1}^n(x_k-\bar{x})(x_k-\bar{x})^T\tag{14} $$ we see that $X-\bar{x}\bar{x}^T$ is positive semi-definite. Therefore, the maximum eigenvalue corresponds to the minimum of $\Delta$.
Once we have found $u$, a unit eigenvector of $X-\bar{x}\bar{x}^T$ with the greatest eigenvalue, we can use $(9)$ to get $p=\bar{x}-\bar{x}\cdot uu$.
Typically, vertical would say all the $(x,y)$ coordinates are the same and horizontal would say all the $z$ coordinates are the same. So you could just look at the standard deviations of the coordinates of all the points to assess vertical or horizontal. That doesn't check if the points lie on an arbitrary line.
The discussion you cite was indeed to give a plane-that was the hypothesis. If you find a relation like $ax+by+cz=k$ you get a plane, as one equation reduces the dimension of the space by one. If you believe your points lie on a line in 3-space, you have two options. A standard linear regression (where one coordinate is fixed and you minimize the sum squared error in the other direction) may well work for you, and you can just do two regressions, one $x$ vs $y$ and another $x$ vs $z$. If the equations you get are $y=m_yx+b_y, z=m_zxb_z$, the line is then $(0,b_y,b_z)+k(1,m_y,m_z)$. If you want orthogonal distance regression, I believe (but didn't follow the derivation to confirm) the link you have is applicable, but you want the line to be in the direction of the maximum singular value, still through the centroid.