使用 Scipy 最小二乘法的 3D 平面拟合产生完全错误的拟合答案

【问题标题】：3D plane fitting using Scipy least squares producing completely wrong fit使用 Scipy 最小二乘法的 3D 平面拟合产生完全错误的拟合
【发布时间】：2019-09-06 18:20:22
【问题描述】：

我从网格表面提取了一组节点的坐标，并将它们放在一个数组中：

[[-2.5      4.       0.     ]
 [-6.5      0.       0.     ]
 [-6.5      0.      20.     ]
 ...
 [-3.5      3.      10.5    ]
 [-3.16667  3.33333 10.5    ]
 [-2.83333  3.66667 10.5    ]]

这个数据集中的点来自一个平面，我想从中获得一个方程。我的 python 代码基于来自 gist 的代码。我正在再次绘制面部以验证计算出的平面是否正确/

我的代码如下（其中数据是上面显示的数组）：

def least_sq(self, data, order=1):

        # regular grid covering the domain of the data
        X,Y = np.meshgrid(np.arange(-3.0, 3.0, 0.5), np.arange(-3.0, 3.0, 0.5))
        XX = X.flatten()
        YY = Y.flatten()

        #1: linear, 2: quadratic
        if order == 1:
            # best-fit linear plane
            A = np.c_[data[:,0], data[:,1], np.ones(data.shape[0])]
            C,_,_,_ = scipy.linalg.lstsq(A, data[:,2])    # coefficients
            print(C)
            # evaluate it on grid
            #Z = C[0]*X + C[1]*Y + C[2]

            #or expressed using matrix/vector product
            Z = np.dot(np.c_[XX, YY, np.ones(XX.shape)], C).reshape(X.shape)

        elif order == 2:
            # best-fit quadratic curve
            A = np.c_[np.ones(data.shape[0]), data[:,:2], np.prod(data[:,:2], axis=1), data[:,:2]**2]
            C,_,_,_ = scipy.linalg.lstsq(A, data[:,2])
            print(C)
            # evaluate it on a grid
            Z = np.dot(np.c_[np.ones(XX.shape), XX, YY, XX*YY, XX**2, YY**2], C).reshape(X.shape)

        # plot points and fitted surface
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.8)
        ax.scatter(data[:,0], data[:,1], data[:,2], c='r', s=50)
        plt.xlabel('X')
        plt.ylabel('Y')
        ax.set_zlabel('Z')
        ax.axis('equal')
        ax.axis('tight')
        plt.show()

我唯一改变的是提供的数据。但是结果不正确，您可以在此处看到：

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

对于平面飞机/lineair：输入您的代码“order = 1” 在 if 行之前

*Z = C[0]*X + C[1]Y + C[2] 未注释！ #Z = np.dot(np.c_[XX, YY, np.ones(XX.shape)], C).reshape(X.shape)

对于样方平面：order = 2

顺便说一句：你用来装配的飞机让我想起了安斯科姆四重奏的情况。

【讨论】：