Установка линии в 3D

Ответ 1

Если вы пытаетесь предсказать одно значение из двух других, то вы должны использовать lstsq с аргументом a как ваши независимые переменные (плюс столбец 1 для оценки перехвата) и b как ваш зависимая переменная.

Если, с другой стороны, вы просто хотите получить наилучшую подходящую строку для данных, то есть линию, которая, если вы проецируете данные на нее, минимизирует квадрат расстояния между реальной точкой и ее проекцией, тогда то, что вы хотите, является первым основным компонентом.

Один из способов определить это линия, вектор направления которой является собственным вектором ковариационной матрицы, соответствующей самому большому собственному значению, которое проходит через среднее значение ваших данных. Тем не менее, eig(cov(data)) - действительно плохой способ вычислить его, поскольку он делает много ненужных вычислений и копирования и потенциально менее точен, чем при использовании svd. См. Ниже:

import numpy as np

# Generate some data that lies along a line

x = np.mgrid[-2:5:120j]
y = np.mgrid[1:9:120j]
z = np.mgrid[-5:3:120j]

data = np.concatenate((x[:, np.newaxis], 
                       y[:, np.newaxis], 
                       z[:, np.newaxis]), 
                      axis=1)

# Perturb with some Gaussian noise
data += np.random.normal(size=data.shape) * 0.4

# Calculate the mean of the points, i.e. the 'center' of the cloud
datamean = data.mean(axis=0)

# Do an SVD on the mean-centered data.
uu, dd, vv = np.linalg.svd(data - datamean)

# Now vv[0] contains the first principal component, i.e. the direction
# vector of the 'best fit' line in the least squares sense.

# Now generate some points along this best fit line, for plotting.

# I use -7, 7 since the spread of the data is roughly 14
# and we want it to have mean 0 (like the points we did
# the svd on). Also, it a straight line, so we only need 2 points.
linepts = vv[0] * np.mgrid[-7:7:2j][:, np.newaxis]

# shift by the mean to get the line in the right place
linepts += datamean

# Verify that everything looks right.

import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as m3d

ax = m3d.Axes3D(plt.figure())
ax.scatter3D(*data.T)
ax.plot3D(*linepts.T)
plt.show()

Вот что это выглядит:

Ответ 2

Если ваши данные достаточно хорошо себя ведут, то этого должно быть достаточно, чтобы найти сумму наименьших квадратов расстояний компонентов. Тогда вы можете найти линейную регрессию с z, не зависящую от x, а затем снова независимую от y.

Следуя документации:

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()

Если вы хотите свести к минимуму фактические ортогональные расстояния от линии (ортогональной линии) до точек в 3-пространстве (что я не уверен, даже упоминается как линейная регрессия). Затем я бы построил функцию, которая вычисляет RSS и использует функцию минимизации scipy.optimize для ее решения.