回帰平面を見つけて点集合に描画する

Question

いくつかのデータ点に飛行機を合わせて描画したいと思います。

あなたは私が手動で飛行機を作成し、現時点では見ることができるように、次になり

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で可能ですね。エラー機能の種類はそれほど重要ではありません。私は最小の正方形（垂直の点 - 平面 - 距離）が良いと思う。もしあなたの誰かが私にそれをする方法を教えてくれれば涼しいでしょう。

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

不必要に尋ねて申し訳ありません。

Answer 2

もう1つの方法は、直線正方向最小二乗法です。 平面の方程式は、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 

：：だから、すべてのデータを使用してこのような行列を設定アックス= 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()