我想将平面拟合到 Python 中的一组点 (x, y, z)。如果误差是相对于 z 轴测量的，我找到了如何执行拟合的各种答案，但我想考虑正交方向的误差。我发现以下问题 (Best fit plane by minimizing orthogonal distances) 解决了相同的问题 - 但我不清楚如何在 Python 中实现它（可能使用 NumPy/SciPy）。有关数学推导的更多详细信息也可以在此处找到：http://www.ncorr.com/download/publications/eberlyleastsquares.pdf（第 2 节）。
【问题讨论】：
标签： python curve-fitting least-squares
我也不得不处理这种情况，起初数学符号可能会让人不知所措，但最终解决方案相当简单。
一旦你得到定义最佳拟合平面 Ax+By+Cz+D=0 的向量 (A,B,C) 的直觉这解释了坐标集的最小方差，那么解决方案就很简单了。
首先要做的是将你的坐标居中（这样D在你的平面方程中将是0）
coords -= coords.mean(axis=0)
那么您有 2 个选项来获取您感兴趣的向量：(1) 使用来自 sklearn 或 scipy 的 PCA 实现来获取解释最小方差的向量
pca = PCA(n_components=3)
pca.fit(coords)
# The last component/vector is the one with minimal variance, see PCA documentation
normal_vector = pca.components_[-1]
(2) 重新实现您链接的几何工具参考中描述的过程。
@njit
def get_best_fitting_plane_vector(coords):
# Calculate the covariance matrix of the coordinates
covariance_matrix = np.cov(coords, rowvar=False) # Variables = columns
# Calculate the eigenvalues & eigenvectors of the covariance matrix
e_val, e_vect = np.linalg.eig(covariance_matrix)
# The normal vector to the plane is the eigenvector associated to the minimum eigenvalue
min_eval = np.argmin(e_val)
normal_vector = e_vect[:, min_eval]
return normal_vector
在速度方面，重新实现的过程比使用PCA要快，如果使用numba会快很多（只需用@njit装饰函数即可）。
【讨论】：
您提供的第一个链接确实描述了正交距离拟合的算法，但相当简洁。如果有帮助，这里有一个更冗长的描述：
我想你有点（在你的情况下是 3d，但维度对算法没有影响） P[i], i=1..N 您想找到一个与您的点具有最小正交距离的（超）平面。
一个超平面可以用一个单位向量n和一个标量d来描述。平面上的点集是
{ P | n.P + d = 0 }
点P到平面的（正交）距离为
n.P + d
所以我们要找到 n 和 d 来最小化
Q(n,d) = Sum{ i | (n.P[i]+d)*(n.P[i]+d) } /N
（除以 N 不是必需的，对找到的 n 和 d 的值没有影响，但在我看来使代数更整洁）
首先要注意的是，如果我们知道 n，最小化 Q 的 d 将是
d = -n.Pbar where
Pbar = Sum{ i | P[i]}/N, the mean of the P[]
我们也可以使用 d 的这个值，这样，经过一点代数，问题就可以归结为最小化 Q^：
Q^(n) = Sum{ i | (n.P[i]-n.Pbar)*(n.P[i]-n.Pbar) } /N
= n' * C * n
where
C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N
Q^ 的形式告诉我们，使 Q^ 最小化的 n 值将是对应于最小特征值的 C 的特征向量。
所以（对不起我不能给出代码但我的python是可鄙的）：
一个/计算
Pbar = Sum{ i | P[i]}/N, the mean of the points
b/计算
C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N, the covariance matrix of the points
c/对角化C，并挑出一个最小特征值和对应的特征向量n
d/计算
d = -Pbar.n
然后n,d定义你想要的超平面。
【讨论】：
基于您的second refernce []
假设你有 n 个样本 (x,y,z)
我将这 3 个术语称为 M*A=V，并定义 column 数组
X=[ x_0, x_1 .. x_n ]'
Y=[ y_0, y_1 .. y_n ]'
Z=[ z_0, z_1 .. z_n ]'
定义（n×3）矩阵XY1=[X,Y,1n]：
[[x_0,y_0,1],
XY1= [x_1,y_1,1],
...
[x_n,y_n,1]]
矩阵M可以得到为
M = XY1' * XY1
其中撇号 (') 是转置运算符， (*) 是矩阵乘积。
而数组V是
V = XY1'*Z
通过moore-penrose pseoudoinverse可以得到最小二乘解：[(M'*M)^-1 * M']
~A = [(M'*M)^-1 * M'] * V
示例代码：
import numpy as np
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
#Input your values
A=3
B=2
C=1
#reserve memory
xy1=np.ones([n,3])
#Make random data, n ( x,y ) tuples.
n=30 #samples
xy1[:,:2]=np.random.rand(n,2)
#plane: A*x+B*y+C = z , the z coord is calculated from random x,y
z=xy1.dot (np.array([[A,B,C],]).transpose() )
#addnoise
xy1[:,:2]+=np.random.normal(scale=0.05,size=[n,2])
z+=np.random.normal(scale=0.05,size=[n,1])
#calculate M and V
M=xy1.transpose().dot(xy1)
V=xy1.transpose().dot(z)
#pseudoinverse:
Mp=np.linalg.inv(M.transpose().dot(M)).dot(M.transpose())
#Least-squares Solution
ABC= Mp.dot(V)
输出
In [24]: ABC
Out[24]:
array([[3.11395111],
[2.02909874],
[1.01340411]])
【讨论】：