<p>最小二乘法应该很容易适应一个平面。一个平面的方程是:ax+by+c=z。所以用所有数据建立这样的矩阵:</p>
<pre><code> x_0 y_0 1
A = x_1 y_1 1
...
x_n y_n 1
</code></pre>
<p>以及</p>
<pre><code> a
x = b
c
</code></pre>
<p>以及</p>
<pre><code> z_0
B = z_1
...
z_n
</code></pre>
<p>换句话说:Ax=B,现在求x的系数。但由于你有3个以上的点,系统被过度确定,所以你需要使用左伪逆。所以答案是:</p>
<pre><code>a
b = (A^T A)^-1 A^T B
c
</code></pre>
<p>下面是一些简单的Python代码和一个例子:</p>
<pre><code>import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET = 5
EXTENTS = 5
NOISE = 5
# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
zs.append(xs[i]*TARGET_X_SLOPE + \
ys[i]*TARGET_y_SLOPE + \
TARGET_OFFSET + np.random.normal(scale=NOISE))
# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')
# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
tmp_A.append([xs[i], ys[i], 1])
tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)
print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print residual
# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
for c in range(X.shape[1]):
Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
</code></pre>
<p>你的观点的解决方案是:</p>
<pre><code>0.143509 x + 0.057196 y + 1.129595 = z
</code></pre>
<p><a href="https://i.stack.imgur.com/kqXzP.gif" rel="noreferrer"><img src="https://i.stack.imgur.com/kqXzP.gif" alt="plane fit"/></a></p>