Finding a slope of 3D linear regression line

Question

I am trying to find a slope of a line in 3D space. The solution for plotting such line is given in this post

Here's the given code from the link above:

import numpy as np

pts = np.add.accumulate(np.random.random((10,3)))
x,y,z = pts.T

# this will find the slope and x-intercept of a plane
# parallel to the y-axis that best fits the data
A_xz = np.vstack((x, np.ones(len(x)))).T
m_xz, c_xz = np.linalg.lstsq(A_xz, z)[0]

# again for a plane parallel to the x-axis
A_yz = np.vstack((y, np.ones(len(y)))).T
m_yz, c_yz = np.linalg.lstsq(A_yz, z)[0]

# the intersection of those two planes and
# the function for the line would be:
# z = m_yz * y + c_yz
# z = m_xz * x + c_xz
# or:
def lin(z):
    x = (z - c_xz)/m_xz
    y = (z - c_yz)/m_yz
    return x,y

#verifying:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
zz = np.linspace(0,5)
xx,yy = lin(zz)
ax.scatter(x, y, z)
ax.plot(xx,yy,zz)
plt.savefig('test.png')
plt.show()

Math-wise I know how to find intersections of two planes and the slope of a given line but I am having trouble putting it in code.

How can I find the slope of resulting regression line using this solution?

Answer 1

The "slope" of a 3D line is generally taken to be slopes of the line "projected" onto the x, y and z planes. See the second answer to this question

If this is what you intended then it's easy enough to calculate these; this modified version of your code below does this into the sx, sy and sz variables:

import numpy as np
import matplotlib.pyplot as plt
from   mpl_toolkits.mplot3d import Axes3D
from   math import pow, sqrt

pts = np.add.accumulate(np.random.random((10,3)))

x, y, z = pts.T

# plane parallel to the y-axis
A_xz = np.vstack((x, np.ones(len(x)))).T
m_xz, c_xz = np.linalg.lstsq(A_xz, z, rcond=None)[0]

# plane parallel to the x-axis
A_yz = np.vstack((y, np.ones(len(y)))).T
m_yz, c_yz = np.linalg.lstsq(A_yz, z, rcond=None)[0]

# the intersection of those two planes and
# the function for the line would be:
# z = m_yz * y + c_yz
# z = m_xz * x + c_xz
# or:
def lin(z):
    x = (z - c_xz)/m_xz
    y = (z - c_yz)/m_yz
    return x,y


# get 2 points on the intersection line 
za = z[0]
zb = z[len(z) - 1]
xa, ya = lin(za)
xb, yb = lin(zb)

# get distance between points
len = sqrt(pow(xb - xa, 2) + pow(yb - ya, 2) + pow(zb - za, 2))

# get slopes (projections onto x, y and z planes)
sx = (xb - xa) / len  # x slope
sy = (yb - ya) / len  # y slope
sz = (zb - za) / len  # z slope

# integrity check - the sum of squares of slopes should equal 1.0
# print (pow(sx, 2) + pow(sy, 2) + pow(sz, 2))

fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x, slope: %.4f" %sx, color='blue')
ax.set_ylabel("y, slope: %.4f" %sy, color='blue')
ax.set_zlabel("z, slope: %.4f" %sz, color='blue')
ax.scatter(x, y, z)
ax.plot([xa], [ya], [za], markerfacecolor='k', markeredgecolor='k', marker = 'o')
ax.plot([xb], [yb], [zb], markerfacecolor='k', markeredgecolor='k', marker = 'o')
ax.plot([xa, xb], [ya, yb], [za, zb], color = 'r')

plt.show()

The output graph below shows the line in question, which is just drawn between the 2 extreme xyz points.

I hope this may help