找到并绘制回归平面到一组点

Posted

技术标签:

【中文标题】找到并绘制回归平面到一组点【英文标题】:Find and draw regression plane to a set of points 【发布时间】:2014-01-09 02:11:35 【问题描述】:

我想将平面拟合到一些数据点并绘制它。我目前的代码是这样的:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

points = [(1.1,2.1,8.1),
          (3.2,4.2,8.0),
          (5.3,1.3,8.2),
          (3.4,2.4,8.3),
          (1.5,4.5,8.0)]

xs, ys, zs = zip(*points)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(xs, ys, zs)

point  = np.array([0.0, 0.0, 8.1])
normal = np.array([0.0, 0.0, 1.0])
d = -point.dot(normal)
xx, yy = np.meshgrid([-5,10], [-5,10])
z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2]
ax.plot_surface(xx, yy, z, alpha=0.2, color=[0,1,0])

ax.set_xlim(-10,10)
ax.set_ylim(-10,10)
ax.set_zlim(  0,10)

plt.show()

结果如下:

正如您现在所看到的,我手动创建了平面。我该如何计算呢?我想scipy.optimize.minimize 是可能的。目前,这种错误函数对我来说并不重要。我认为最小二乘(垂直点平面距离)会很好。如果你们中的一个人能告诉我怎么做,那就太棒了。

【问题讨论】:

请参考Best fit plane algorithms why different results了解几种可能的方法 【参考方案1】:

哦,我突然想到了这个主意。这很容易。 :-)

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import scipy.optimize
import functools

def plane(x, y, params):
    a = params[0]
    b = params[1]
    c = params[2]
    z = a*x + b*y + c
    return z

def error(params, points):
    result = 0
    for (x,y,z) in points:
        plane_z = plane(x, y, params)
        diff = abs(plane_z - z)
        result += diff**2
    return result

def cross(a, b):
    return [a[1]*b[2] - a[2]*b[1],
            a[2]*b[0] - a[0]*b[2],
            a[0]*b[1] - a[1]*b[0]]

points = [(1.1,2.1,8.1),
          (3.2,4.2,8.0),
          (5.3,1.3,8.2),
          (3.4,2.4,8.3),
          (1.5,4.5,8.0)]

fun = functools.partial(error, points=points)
params0 = [0, 0, 0]
res = scipy.optimize.minimize(fun, params0)

a = res.x[0]
b = res.x[1]
c = res.x[2]

xs, ys, zs = zip(*points)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(xs, ys, zs)

point  = np.array([0.0, 0.0, c])
normal = np.array(cross([1,0,a], [0,1,b]))
d = -point.dot(normal)
xx, yy = np.meshgrid([-5,10], [-5,10])
z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2]
ax.plot_surface(xx, yy, z, alpha=0.2, color=[0,1,0])

ax.set_xlim(-10,10)
ax.set_ylim(-10,10)
ax.set_zlim(  0,10)

plt.show()

对不起,不必要地问了。

【讨论】:

【参考方案2】:

另一种方法是直接使用最小二乘法。 平面的方程是:ax + by + c = z。所以用你的所有数据设置这样的矩阵:

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1  

    a  
x = b  
    c

    z_0   
B = z_1   
    ...   
    z_n

换句话说:Ax = B。现在求解x,它们是你的系数。但是因为(我假设)你有超过 3 个点,系统是超定的,所以你需要使用左伪逆。所以答案是:

a 
b = (A^T A)^-1 A^T B
c

下面是一些简单的 Python 代码和示例:

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

【讨论】:

很好的答案,感谢分享!由于不推荐使用 np.matrix,因此我编辑了代码,使其适用于 np 数组【参考方案3】:

感谢@Ben 的分享!由于不推荐使用 np.matrix,因此我编辑了您的代码,使其适用于 np 数组

import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import inv

# Pass the function array of points, shape (3, X)
def plane_from_points(points):
    # Create this matrix correctly without transposing it later?
    A = np.array([
        points[0,:],
        points[1,:],
        np.ones(points.shape[1])
    ]).T
    b = np.array([points[2, :]]).T

    # fit = (A.T * A).I * A.T * b
    fit = np.dot(np.dot(inv(np.dot(A.T, A)), A.T), b)
    # errors = b - np.dot(A, fit)
    # residual = np.linalg.norm(errors)
    return fit

N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 3

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

points = np.array([xs, ys, zs])

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

fit = plane_from_points(points)
# 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()

【讨论】:

以上是关于找到并绘制回归平面到一组点的主要内容,如果未能解决你的问题,请参考以下文章

如何使用 Python 在一组点上绘制多边形(部分向内弯曲)边缘?

找到一个具有最大点数的圆 ((x,y,r));给定二维平面中的一组点(x,y)

将点列表排序为多边形

将平面拟合到 3D 中的一组点:scipy.optimize.minimize vs scipy.linalg.lstsq

如何在 Matplotlib 中绘制模糊点

如何在 scikit learn 中为多元回归绘制最佳拟合平面?