python 最小二乘法三维坐标拟合平面_【MQ笔记】超简单的最小二乘法拟合平面（Python）...

这篇笔记中，我主要通过解决“由离散点拟合平面”这个小问题，学习了超定方程最小二乘解的求解方法。在这里我整理了两种求解思路用以交流。

直接求解超定方程。

我们知道，对于一个平面，其方程可以用

来表示。由离散点拟合平面，实际上就是求解超定方程：

上述方程可以用

来表示。由于A是一个

的矩阵，因此我们先在等号两边分别乘以 A的转置矩阵

，使系数矩阵变为

的方阵，之后，通过乘以系数矩阵的逆矩阵求解，也就是说，

。

该方法的Python代码如下：

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# 创建函数，用于生成不同属于一个平面的100个离散点

def not_all_in_plane(a, b, c):

x = np.random.uniform(-10, 10, size=100)

y = np.random.uniform(-10, 10, size=100)

z = (a * x + b * y + c) + np.random.normal(-1, 1, size=100)

return x, y, z

# 调用函数，生成离散点

x, y, z = not_all_in_plane(2, 5, 6)

#创建系数矩阵A

a = 0

A = np.ones((100, 3))

for i in range(0, 100):

A[i, 0] = x[a]

A[i, 1] = y[a]

a = a + 1

#print(A)

#创建矩阵b

b = np.zeros((100, 1))

a = 0

for i in range(0, 100):

b[i, 0] = z[a]

a = a + 1

#print(b)

#通过X=(AT*A)-1*AT*b直接求解

A_T = A.T

A1 = np.dot(A_T,A)

A2 = np.linalg.inv(A1)

A3 = np.dot(A2,A_T)

X= np.dot(A3, b)

print('平面拟合结果为：z = %.3f * x + %.3f * y + %.3f'%(X[0,0],X[1,0],X[2,0]))

#计算方差

R=0

for i in range(0,100):

R=R+(X[0, 0] * x[i] + X[1, 0] * y[i] + X[2, 0] - z[i])**2

print ('方差为：%.*f'%(3,R))

# 展示图像

fig1 = plt.figure()

ax1 = fig1.add_subplot(111, projection='3d')

ax1.set_xlabel("x")

ax1.set_ylabel("y")

ax1.set_zlabel("z")

ax1.scatter(x,y,z,c='r',marker='o')

x_p = np.linspace(-10, 10, 100)

y_p = np.linspace(-10, 10, 100)

x_p, y_p = np.meshgrid(x_p, y_p)

z_p = X[0, 0] * x_p + X[1, 0] * y_p + X[2, 0]

ax1.plot_wireframe(x_p, y_p, z_p, rstride=10, cstride=10)

plt.show()

运行结果为：

利用最小二乘法公式求导

通过离散点拟合平面，也就是说，要找到一个平面(

)，使这平面到各个点的“距离”最近，根据最小二乘法，

，也就是说我们要求得一组a,b,c,使得对于已有的离散点来说，S的值最小。

要使得S的值最小，则有，也就是说。整理，得，将其改为矩阵方程,得：。

求解该恰定方程即可得到a，b，c。上述方程也可以用

表示，该方程可以通过两边同时乘以系数矩阵的逆矩阵求得，即

。

该方法的Python代码如下：

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# 创建函数，用于生成不同属于一个平面的100个离散点

def not_all_in_plane(a, b, c):

x = np.random.uniform(-10, 10, size=100)

y = np.random.uniform(-10, 10, size=100)

z = (a * x + b * y + c) + np.random.normal(-1,1,size=100)

return x, y, z

# 调用函数，生成离散点

x2, y2, z2 = not_all_in_plane(2, 5, 6)

#创建系数矩阵A

A=np.zeros((3,3))

for i in range(0,100):

A[0,0]=A[0,0]+x2[i]**2

A[0,1]=A[0,1]+x2[i]*y2[i]

A[0,2]=A[0,2]+x2[i]

A[1,0]=A[0,1]

A[1,1]=A[1,1]+y2[i]**2

A[1,2]=A[1,2]+y2[i]

A[2, 0] = A[0,2]

A[2, 1] = A[1, 2]

A[2, 2] = 100

#print(A)

#创建b

b = np.zeros((3,1))

for i in range(0,100):

b[0,0]=b[0,0]+x2[i]*z2[i]

b[1,0]=b[1,0]+y2[i]*z2[i]

b[2,0]=b[2,0]+z2[i]

#print(b)

#求解X

A_inv=np.linalg.inv(A)

X = np.dot(A_inv, b)

print('平面拟合结果为：z = %.3f * x + %.3f * y + %.3f'%(X[0,0],X[1,0],X[2,0]))

#计算方差

R=0

for i in range(0,100):

R=R+(X[0, 0] * x2[i] + X[1, 0] * y2[i] + X[2, 0] - z2[i])**2

print ('方差为：%.*f'%(3,R))

# 展示图像

fig1 = plt.figure()

ax1 = fig1.add_subplot(111, projection='3d')

ax1.set_xlabel("x")

ax1.set_ylabel("y")

ax1.set_zlabel("z")

ax1.scatter(x2,y2,z2,c='r',marker='o')

x_p = np.linspace(-10, 10, 100)

y_p = np.linspace(-10, 10, 100)

x_p, y_p = np.meshgrid(x_p, y_p)

z_p = X[0, 0] * x_p + X[1, 0] * y_p + X[2, 0]

ax1.plot_wireframe(x_p, y_p, z_p, rstride=10, cstride=10)

plt.show()

运行结果为：