平面の計算方法

質問

3つ以上の座標点(n個)から最小二乗法を用いて，平面の中心座標，XY・XZ・YZ方向の傾きなどを計算したいのですが，どのように計算したらよろしいでしょうか？

回答No.4

age_momo · ベストアンサー

質問者さんが誤差をどの方向に考えているかで変わってくると思います。

全ての点の平面からの距離の二乗和の最小を考えるなら皆さんが回答されている通りだと思います。
ただ、単純な最小2乗法は全ての点から直線までの距離が最小になるよう定めているわけ
ではありません。y方向に誤差を考えていて、その2乗和が最小になるy=a+bxのa,bを
求めています。回帰の説明などを見れば明らかですね。点から直線に伸ばす矢印(誤差の量)は
y軸に平行に書かれています。

http://szksrv.isc.chubu.ac.jp/lms/lms1.html

直線までの距離なら直線に対して垂直に書かなければなりません。
その場合はこの質問同様、もう少し複雑な式になります。

これと同じに考えるなら(つまりz軸方向への誤差の最小を考える)単純な重回帰に
なります。n個のデータ(xi,yi,zi)を使ってz=ax+by+cに回帰するには

aΣxi^2+bΣxiyi+cΣxi=Σxizi
aΣxiyi+bΣyi^2+cΣyi=Σyizi
aΣxi+bΣyi+cn=Σzi

の連立方程式を解いてください。なお、ここで例えばΣxi^2は

Σxi^2=x1^2+x2^2+x3^2+・・・・・＋xn^2

でデータが得られているなら定数です。

お礼

age_momoさん，ありがとうございます．

ｎ個のデータ(xi,yi,zi)を使ってz=ax+by+cに回帰させるとき，
>　aΣxi^2+bΣxiyi+cΣxi=Σxizi
>　aΣxiyi+bΣyi^2+cΣyi=Σyizi
>　aΣxi+bΣyi+cn=Σzi
この3つの式を行列式に置き換えて計算したいのですが，
計算式は，
| Σxi^2　Σxiyi　Σxi || a | | Σxizi 　|
| Σxiyi　Σyi^2　Σyi || b |=| Σyizi |
| Σxi　　Σyi　　n　 　 || c | | Σzi 　 |
で大丈夫でしょうか？

回答No.6

age_momo

>この3つの式を行列式に置き換えて計算したいのですが，
>計算式は，
>| Σxi^2　Σxiyi　Σxi || a | | Σxizi 　|
>| Σxiyi　Σyi^2　Σyi || b |=| Σyizi |
>| Σxi　　Σyi　　n　 　 || c | | Σzi 　 |
>で大丈夫でしょうか？

行列式というか行列そのもので
A=
| Σxi^2　Σxiyi　Σxi |
| Σxiyi　Σyi^2　Σyi |
| Σxi　　Σyi　　n　 |

の逆行列A^(-1)を計算して
A^(-1)| Σxizi 　|
| Σyizi |
| Σzi 　|

を計算すればいいです。

お礼

返事が遅くなってしまい申し訳ありません．

質問の趣旨がおかしくなってしまい，本来聞きたかったことが聞けませんでした．
聞きたかったことは，

aΣxi^2+bΣxiyi+cΣxi=Σxizi
aΣxiyi+bΣyi^2+cΣyi=Σyizi
aΣxi+bΣyi+cn=Σzi

↑の式から平面の中心座標(a,b,c)を求めるには，↑の式を行列に置き換えて，クラーメルの公式を使って求める方法でも(a,b,c)がでるかどうかとういことです．

本当に時間がかかってしまい，申し訳ありません．

補足

>A=
>| Σxi^2　Σxiyi　Σxi |
>| Σxiyi　Σyi^2　Σyi |
>| Σxi　　Σyi　　n　 |

>の逆行列A^(-1)を計算して
>A^(-1)| Σxizi 　|
>| Σyizi |
>| Σzi 　|

これらの式はクラーメルの公式だと思うのですが，
求まったa，b，cの値はそのまま中心座標になりますか？

試しに自分で考えた測定点で，エクセルを使って計算しているのですが答えが求まりません．
10　10　10
20　10　10
10　10　20
20　10　20

お時間がありましたら，教えてください．

回答No.5

inara

ANo.3です。
実は私も仕事で同じような計算をしたことがあります。具体的に言うと、数cm角の平面板が弾性体で支えられている試作品がいっぱいあって、その平面の傾斜分布を定量化するために、平面上の４点以上の点の高さをレーザ変位計で測定し、age_momoさんの連立方程式を解いて、最小２乗平面をそれぞれの試作品について求めました。

最終的には、法線ベクトルを方位角と天頂角の２変数に変換して、極座標グラフで分布（バラツキ）を表示させ、傾向の分析を行いましました。測定した平面は、傾斜しているといっても、そのバラツキは小さく、xyz 座標でいえば、法線ベクトルが z 軸から数度以内だったので、「測定点と平面のｚ座標の差の２乗」が最小となるようにして平面を決めていました（平面の法線ベクトルとz軸との角度が小さい場合は |z座標の差| ≒ 点と平面の距離 と考えていいので）。

しかし、求めたい平面がｚ軸にほぼ平行である場合を含め、向きが全く分からない場合は、点と平面の距離の２乗を最小とする方法が良いのではと判断して、ANo.3ではその方法を紹介しました。私は回帰分析の専門家ではないので数学的にどれが最良かというのはお答えできませんが、ある程度、平面の向きがわかっている場合には、計算が比較的簡単なので、age_momoさんの方法でも良いかと思います。

回答No.3

inara

考え方としてはANo.2の参考URLにあるとおりですが大変面倒な計算になります。
うまく収束するかどうか分かりませんが、Excelのソルバーを使って、残差２乗和が最小となるような、５個の未知パラメータを求めるのはどうでしょうか。

結果だけ書きますが、n個のデータ点と平面との距離の２乗の和 S が最小となるような次式の a, b, x0, y0, z0 を求めれば、平面が決定されます。

S = Σ [ i = 1, 2, ... n ] [ a*( xi - x0 ) + b*( yi - y0 ) + √{ 1 - ( a^2 + b^2 ) } *( zi - z0 ) ]^2

( xi, yi, zi ) はn個のデータ点（i = 1,2,.. n ）です。a、b は平面の法線ベクトル（平面に垂直なベクトル）のx成分とy成分です。（ x0, y0 , z0 ) この平面が通る点の座標です。法線ベクトルのz成分がないのは、a^2 + b^2 + c^2 = 1 と規格化したので、c = √{ 1 - ( a^2 + b^2 ) } と自動的に決まるからです。なぜこのような式になるかについては【詳細】に書きました。

未知パラメータ a, b, x0, y0, z0 のうち、x0, y0, z0 は、ANo.1さんのコメントにある「中心座標の計算はn個の座標を足してnで割れば」というのは、常に成り立つかどうか分かりません。ただ、初期値として、x0 = Σ(xi)、y0 = Σ(yi)、z0 = Σ(zi) とするのは良いと思います（Excelのソルバーは初期値の選択を誤ると収束しないので、もっともらしい値を初期値とすれば収束しやすくなります）。

【詳細】
法線ベクトル（平面に垂直なベクトル）が (a,b,c) で、点（x0,y0,z0) を通る平面の方程式は次式で表されます。
a*( x - x0 )+b*( y - y0 )+c*( z - z0 ) = 0 --- (1)
一方、n 個の点　( xi, yi, zi ) [ i = 1, 2, ... n ] と平面との距離を Li [ i = 1, 2, ... n ] とすれば
Li^2 = ｛a*( xi - x0 ) + b*( yi - y0 ) + c*( zi - z0 ) }^2 / ( a^2 + b^2 + c^2 ) ---- (2)
となります。したがって、n個の点と平面(1)との距離の２乗の総和を S とすれば、
S = Σ [ i = 1, 2, ... n ] ( Li^2 ) = Σ [ i = 1, 2, ... n ] [ ｛a*( xi - x0 ) + b*( yi - y0 ) + c*( zi - z0 ) }^2 / ( a^2 + b^2 + c^2 ) ] --- (3)
となります。これが最小 となるような、a, b, c, x0, y0, z0 を求めるというのが最終目的です。

S が最小 であるとき、∂S/∂a = ∂S/∂b = ∂S/∂c = ∂S/∂x0 = ∂S/∂y0 = ∂S/∂z0 = 0 が成り立ちます。Sを偏微分すると、次の式が得られます（和の範囲はいつも [ i = 1, 2, ... n ] なので省略します)。

a*( b^2+c^2 )*Σ( xi - x0 )^2 + b*Σ( xi - x0 )*( yi - y0 ) + c*Σ( zi - z0 )*( xi - x0 ) = 0 --- (4)
b*( c^2+a^2 )*Σ( yi - y0 )^2 + a*Σ( xi - x0 )*( yi - y0 ) + c*Σ( yi - y0 )*( zi - z0 ) = 0 --- (5)
c*( a^2+b^2 )*Σ( zi - z0 )^2 + b*Σ( yi - y0 )*( zi - z0 ) + a*Σ( zi - z0 )*( xi - x0 ) = 0 --- (6)
a*Σ( xi - x0 ) + b*Σ( xi - x0 ) + c*Σ( zi - z0 ) = 0 --- (7)

（）などを展開などして、この計算を続ければANo.2の参考URLにある式(4)に行く着くでしょうが、大変複雑な結果となります。
手っ取り早くには、式(3)の S が最小になるような、a, b, c, x0, y0, z0 を求めればいいわけです。

未知パラメータは６個あるように思えますが、実は５個です。a,b,c は比が分かればいいので、 a^2 + b^2 + c^2 = 1 と規格化してしまいます（とすれば |a|≦1、 |b|≦1、 |c|≦1 です）。すると、c^2 = 1 - ( a^2 + b^2 ) ≧ 0 ですから、a,bが決まれば c が決まりますので、c は求めなくていいのです。c の符号は±の２通りありますが、法線ベクトルというのは 、平面の表と裏を区別しないなら、(a,b,c) も (-a,-b,-c) も同じことですから、cの符号は＋か－かどちらに決めてしまえばいいのです。したがって c > 0 とすることにすれば、c = √{ 1 - ( a^2 + b^2 ) } です。すると式(3)は

S = Σ [ i = 1, 2, ... n ] [ a*( xi - x0 ) + b*( yi - y0 ) + √{ 1 - ( a^2 + b^2 ) } *( zi - z0 ) ]^2 ---- (2')

となります。

なお、ANo1さんの「中心座標の計算はn個の座標を足してnで割れば」というのは、式(7)で、Σ( xi - x0 ) =0、 Σ( xi - x0 ) = 0
Σ( zi - z0 ) = 0 としてx0, y0, z0 を求めるものですが、Σ() <> 0でも式(7)が成り立つ場合もあるかと思います。

回答No.2

noname#29127

検索して下記の文章を見つけました。

http://www.aichi-inst.jp/html/reports/repo2001/R01_11.pdf

この部分の2.1の部分で、最小二乗平面(Z=Ax + By + C)を求める
方法が記載されています。(4)式を解いてA,B,Cを求めるようです。

回答No.1

noname#29127

間違っているかもしれないので、参考意見として

中心座標の計算はn個の座標を足してnで割れば、中心が出ますね。

XY, XZ, YZ方向の傾きはXYの場合、X座標とY座標だけ用いて
最小二乗法で傾きが出ますね。他のXZ, YZでも同様に。

ただ、この傾きが実際どういうことに使えるのか不明です。

多くの点があり、その平均的な所を通る平面を求めたいのでは
ないでしょうか？
その場合はXY方向でなく、「XY軸上での傾き」などになるのでしょうね。

平均的な所を通る平面ですが、互いに垂直なベクトルが求まれば、
平面が定義できそうですが、それをどう求めるのか、ちょっと
分かっていません。

中心を求めた後に、中心の点から各座標へのベクトルを足し合わせる
といった感じでできそうな気がしますが。

お礼

popiahさん，ありがとうございます.
popiahさんが仰るとおり，多くの点があって，それらの点の平均的な所を通る平面を求めたいのです.

　最小二乗法で，平均的な所を通る平面を求める計算式などはやはりないのでしょうか？