$x_{i},y_{i},z_{i}\in\mathbf{x},\mathbf{y},\mathbf{z}\in\mathbb{R}^{N}
,a,b,cx,y,z\in\mathbb{R}
,$ モデル関数$z=ax+by+c$
及びデータ$d_{i}=ax_{i}+by_{i}+c-z_{i}d$
をそれぞれ定義する。モデル関数とデータの誤差関数を
E \equiv \sum_{i=0}^{N}\left(ax_{i}+by_{i}+c-z_{i}\right)^{2}
とする。
\frac{\partial E}{\partial a} = \sum_{i=0}^{N}2x_{i}\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
\frac{\partial E}{\partial b} = \sum_{i=0}^{N}2y_{i}\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
\frac{\partial E}{\partial c} = \sum_{i=0}^{N}2\left(ax_{i}+by_{i}+c-z_{i}\right)=0\\
したがって、
a\sum_{i=0}^{N}x_{i}^{2}+b\sum_{i=0}^{N}x_{i}y_{i}+c\sum_{i=0}^{N}x_{i}-\sum_{i=0}^{N}x_{i}z_{i} = 0\\
a\sum_{i=0}^{N}x_{i}y_{i}+b\sum_{i=0}^{N}y_{i}^{2}+c\sum_{i=0}^{N}y_{i}-\sum_{i=0}^{N}y_{i}z_{i} = 0\\
a\sum_{i=0}^{N}x_{i}+b\sum_{i=0}^{N}y_{i}+Nc-\sum_{i=0}^{N}z_{i} = 0
上式を下記のように表記する。
aA+bC+cF-D = 0\\
aC+bB+cG-E = 0\\
aF+bG+cN-H = 0\\
ここで、
A\equiv\sum_{i=0}^{N}x_{i}^{2},\\
B\equiv\sum_{i=0}^{N}y_{i}^{2},\\
C\equiv\sum_{i=0}^{N}x_{i}y_{i},\\
D\equiv\sum_{i=0}^{N}x_{i}z_{i},\\
E\equiv\sum_{i=0}^{N}y_{i}z_{i},\\
F\equiv\sum_{i=0}^{N}x_{i},\\
G\equiv\sum_{i=0}^{N}y_{i},\\
H\equiv\sum_{i=0}^{N}z_{i}
。上式を行列表示し、行基本変形を用いて下記のように$a,b, c
$を求める。
\left[\begin{array}{cccc}
A & C & F & -D\\
C & B & G & -E\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & C/A & F/A & -D/A\\
C & B & G & -E\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & B-\frac{C^{2}}{A} & G-C\frac{F}{A} & -E+C\frac{D}{A}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{G-C\frac{F}{A}}{B-\frac{C^{2}}{A}} & -\frac{E-C\frac{D}{A}}{B-\frac{C^{2}}{A}}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
F & G & N & -H
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & G-F\frac{C}{A} & N-\frac{F^{2}}{A} & -H+F\frac{D}{A}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] = 0\\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\left(G-F\frac{C}{A}\right)\frac{AG-CF}{AB-C^{2}} & -H+F\frac{D}{A}+\left(G-F\frac{C}{A}\right)\frac{AE+CD}{AB-C^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] \\
\left[\begin{array}{cccc}
1 & \frac{C}{A} & \frac{F}{A} & -\frac{D}{A}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\frac{\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & -H+F\frac{D}{A}+\frac{\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] \\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & N-\frac{F^{2}}{A}-\frac{\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & -H+F\frac{D}{A}+\frac{\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] \\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & \frac{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}{A^{2}B-AC^{2}} & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE+CD\right)}{A^{2}B-AC^{2}}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] \\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & 1 & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right] \\
\left[\begin{array}{cccc}
1 & 0 & \frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}} & -\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\\
0 & 1 & \frac{AG-CF}{AB-C^{2}} & -\frac{AE-CD}{AB-C^{2}}\\
0 & 0 & 1 & \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)\left(AG-CF\right)}
\end{array}\right]\left[\begin{array}{c}
a\\
b\\
c\\
1
\end{array}\right]
上式より、
-c = \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)^{2}}\\
= \frac{\left(-AH+FD\right)\left(AB-C^{2}\right)+\left(AG-FC\right)\left(AE-CD\right)}{\left(AN-F^{2}\right)\left(AB-C^{2}\right)-\left(AG-FC\right)^{2}}\\
= \frac{IJ+KL}{MJ-K^{2}}\\
c = -\frac{IJ+KL}{MJ-K^{2}}\\
-b = c\frac{AG-CF}{AB-C^{2}}-\frac{AE-CD}{AB-C^{2}}\\
= c\frac{K}{J}-\frac{L}{J}\\
= \frac{cK-L}{J}\\
b = \frac{L-cK}{J}\\
-a = c\left(\frac{F}{A}-\frac{ACG-C^{2}F}{A^{2}B-AC^{2}}\right)+\left(-\frac{D}{A}+\frac{ACE-C^{2}D}{A^{2}B-AC^{2}}\right)\\
= c\left(\frac{F}{A}-\frac{C\left(AG-CF\right)}{A\left(AB-C^{2}\right)}\right)+\left(-\frac{D}{A}+\frac{C\left(AE-CD\right)}{A\left(AB-C^{2}\right)}\right)\\
= c\left(\frac{F}{A}-\frac{CK}{AJ}\right)+\left(-\frac{D}{A}+\frac{CL}{AJ}\right)
= c\left(\frac{FJ-CK}{AJ}\right)+\left(\frac{-DJ+CL}{AJ}\right)\\
= \frac{cFJ-cCK-DJ+CL}{AJ}
= \frac{\left(cF-D\right)J-C\left(cK-L\right)}{AJ}\\
a = \frac{C\left(cK-L\right)-\left(cF-D\right)J}{AJ}\\
\\
I \equiv \left(-AH+FD\right)\\
J \equiv \left(AB-C^{2}\right)\\
K \equiv \left(AG-CF\right)\\
L \equiv \left(AE-CD\right)\\
M \equiv \left(An-F^{2}\right)\\
下記実行結果。総和を要素が全て1のカーネルとの畳込みにしています。
from scipy.signal import *
from scipy.linalg import *
from numpy import *
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
x,y = mgrid[0:100:10,0:100:10] * 1.
#x = (x.flatten() + random.rand(x.size)*100).reshape(x.shape)
#y = (y.flatten() + random.rand(y.size)*100).reshape(y.shape)
a = 2.00
b = 2.00
c = 0.
z = a * x + b * y + c
#k = exp(-((x-50)**2+(y-50)**2)/800)
#z = -sqrt((x-50) ** 2 + (y-50) ** 2) + 100
#z = zeros_like(x)
h = ones((20,20))
n = x.size
A = fftconvolve((x**2),h,"same")
B = fftconvolve((y**2),h,"same")
C = fftconvolve((x*y),h,"same")
D = fftconvolve(x*z,h,"same")
E = fftconvolve(y*z,h,"same")
F = fftconvolve(x,h,"same")
G = fftconvolve(y,h,"same")
H = fftconvolve(z,h,"same")
I = (-A*H+F*D)
J = (A*B-C**2)
K = (A*G-C*F)
L = (A*E-C*D)
M = (A*n-F**2)
c = -(I*J+K*L)/(M*J-K**2)
b = (L-c*K)/J
a = (C*(c*K-L)-(c*F-D)*J)/(A*J)
# 法線ベクトルN[(a, b, d)]
norms = sqrt(a**2 + b ** 2 + 1)
a = a/norms
b = b/norms
d = -ones_like(b) / norms
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c="r", marker="o")
ax.quiver(x.flatten(), y.flatten(), z.flatten(), a.flatten(),b.flatten(),d.flatten(), length=10)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
処理結果(平面)
点を密にすると苔むしたみたいになります(円錐)。
