ｎ次曲線による最小二乗法につていは説明しましたが、円や楕円、その他の一般式についても最小二乗法による近似は可能です。
今回は円の最小二乗法を例にとって説明しますが、他の一般式についても要領は同じです。

まず初めに近似する一般式を作成します。

円の場合、円の中心座標（a,b）、円の半径 r とすると

(X – a)² + (Y – b)² = r² ・・・①

となります。

ｎ次曲線による最小二乗法ではＸ座標がX_iの時のＹ座標Y_iと近似曲線上のＹ座標f(X_i)との差の二乗が最小になるように計算していましたが、円の場合、Ｘ座標がX_iの時にＹ座標が存在しなかったり、Ｙ座標が２点存在してしまう場合もあります。

そこで、今回はＹ座標の差の２乗の総和を求めるのではなく、①式を＝０の式に変形し、
点の座標(X_i,Y_i)を代入し、その式の２乗の総和を求めます。
（ｎ次式近似でやった計算も結局は近似式＝０になるように変形していた訳ですが．．．）

①式を＝０になるように変形し、２乗の総和は

∑｛(X_i – a)² + (Y_i – b)² – r²｝² = 0　・・・②

となります。

ただし、このまま未知数a,b,rに関して②式を偏微分してもa,b,rは４次関数となるので、各未知数の偏微分=0の時が最小とは限りません。

そこで、②式を展開し、下の式のように置きます。

∑｛X_i² + Y_i² + AX_i + BY_i + C｝² = 0　・・・③

ただし

A = -2a ・・・④
B = -2b ・・・⑤
C =  a² + b² – r² ・・・⑥

③式に関してＡ，Ｂ，Ｃについて偏微分すると

∂/∂A = A∑X_i² + B∑X_iY_i + C∑X_i + ∑X_i³ + ∑X_iY_i² = 0 ・・・⑦
∂/∂B = A∑X_iY_i + B∑Y_i² + C∑Y_i + ∑X_i²Y_i + ∑Y_i³ = 0 ・・・⑧
∂/∂C = A∑X_i + B∑Y_i + C∑1  + ∑X_i² + ∑Y_i² = 0 ・・・⑨

⑦、⑧、⑨式を行列を用いて解くと

となり、逆行列を求めるとＡ，Ｂ，Ｃが求まります。

さらに④、⑤、⑥式より円の中心座標（a,b）、円の半径 rが求まります。

  \(a=-\frac{A}{2}\)

  \(b=-\frac{B}{2}\)

  \(r=\sqrt{a^{2}+b^{2}-C}\)

←使える数学へ戻る

関連記事

最小自乗法という人もいますが、私は最小二乗法派です。

最小二乗法とは？

点の集まりから、近似直線（曲線）などの近似式を解くための手法です。
単に直線や二次曲線、三次曲線などに限らず、最小二乗法の応用範囲は広くあります。

最小二乗法は、言葉の通りに２乗を最小にする方法なのですが、この２乗と言うのが、方程式上の理論値と近似するデータとの差の２乗の合計を最小にする方法となります。と言っても分かりづらいので、具体的な最小二乗法の例を示します。

最小二乗法の計算方法

今回は説明を簡単にするため、ｆ（X）= aX + b の直線でデータを近似する場合で説明します。

X = X_i のとき、直線上のY座標は

$$Y_{i}=f(X_{i})=aX_{i}+b$$

より、点と近似直線との差は

$$Y_{i}-f(X_{i})=Y_{i}-(aX_{i}+b)$$

となります。

この差の合計が一番小さい時にもっとも点の集まりを直線で近似できていることとなりますが、
そのままこの差の総和を計算すると差の値がプラスの時とマイナスの時があるため、例えば下図

のように、正の差と負の差が同じようにある場合、差の合計が小さくなってしまうため、差を２乗してから合計した値が最小になるようにします。

よって、差の２乗の合計は

$$\sum(Y_{i}-(aX_{i}+b))^{2}\\
=\sum(a^{2}X_i^2+2aX_{i}b+b^{2}-2aX_{i}Y_{i}-2bY_{i}+Y_i^2 )$$

・・・式①

となり、この総和の値が最小となるａとｂを求めれば近似直線を求められることとなります。

この式を良くみると、ａとｂの二次関数でできているため、この二次関数が最小となる値を求めればよい。
二次関数が最小となるのは微分した値 = 0 の場合であるから、式①をａとｂで偏微分した値が
０（ゼロ）となる値を求めればよい。

式①をaとbでそれぞれ偏微分すると

$$\frac{\partial}{\partial a}=2a\sum X_i^2+2b\sum X_{i}-2\sum X_{i}Y_{i}=0\\
\frac{\partial}{\partial b}=2a\sum X_i+2b\sum 1-2\sum Y_{i}=0$$

より

$$a\sum X_i^2+b\sum X_{i}-\sum X_{i}Y_{i}=0\\
a\sum X_i+b\sum 1-\sum Y_{i}=0$$

となり、未知数がaとbで式が２本あることから、aとbの連立法的式を解くことで未知数a,bを求めることができます。

連立方程式を行列で表すと

$$\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\frac{1}{\sum X_i^2\sum 1-\sum X_{i}\sum X_{i}}\left(\begin{array}{c}\sum 1 & -\sum x_i\\ -\sum x_i & \sum X_i^2\end{array}\right)\left(\begin{array}{c}\sum X_iY_i\\ \sum Y_i\end{array}\right)$$

となり、行列を解くと未知数ａ，ｂが求まり、近似直線を求めることができます。

ちなみに、二次式（aX² + bX + c）で近似すると

$$\left(\begin{array}{c}a\\ b \\ c\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

三次式（aX³ + bX² + cX + d）で近似すると

$$\left(\begin{array}{c}
a\\
b\\
c\\
d
\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^6 & \sum X_i^5 & \sum X_i^4 & \sum X_i^3\\
\sum X_i^5 & \sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^4 & \sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^3 & \sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^3Y_i\\
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

と規則的に変化するため、ｎ次式近似のプログラムも簡単に作ることができます。
ただし、行列の部分を二次元配列と考えると、近似する次数（ｎ）が上がって行くと、配列の先頭に値が追加されていくのは、扱いにくいので、ｎ次式近似のプログラムを作る場合は、近似式を

$$Y=a+bX+cX^2+dX^3 \cdot\cdot\cdot\cdot$$

として、行列の表現は

$$\left(\begin{array}{c}a\\ b\\ c\\d\\:
\end{array}\right)=
\left(\begin{array}{c}
\sum 1 & \sum X_i & \sum X_i^2 & \sum X_i^3 & .. \\
\sum X_i & \sum X_i^2 & \sum X_i^3 & \sum X_i^4 & .. \\
\sum X_i^2 & \sum X_i^3 & \sum X_i^4 & \sum X_i^5 & ..\\
\sum X_i^3 & \sum X_i^4 & \sum X_i^5 & \sum X_i^6 & ..\\
: & : & : & : & :
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum Y_i\\
\sum X_iY_i \\
\sum X_i^2Y_i\\
\sum X_i^3Y_i\\
:
\end{array}\right)$$

として、ｎ次式で近似した最小二乗法を求めます。

この計算をエクセルで行ったものはExcelによる最小二乗近似からダウンロードできます。
逆行列を求めるプログラムにはガウスの消去法などが用いられます。

最小二乗法での注意点

三次式（a+bX+cX²+dX³）の近似を例にとってみると、近似処理をする過程でXに関する計算式はX～X⁶まで計算されます。

ここで、例えばXの値の範囲が

０．０１～１０００

だったとすると、Xに関する計算では

０．０００００００００００1～１００００００００００００００００００

までの計算がされます。
逆行列の計算ではこれらの値に関して足し算や掛け算がなされる訳ですからコンピュータの
計算誤差（桁落ち、情報落ちなど？）が発生する場合があります。
これを防ぐために、最小二乗法のプログラムでは近似式をそのまま使うのではなく、Xの値から
Xの平均値（X_ｂ）を引いた計算式

$$Y=a+b(X-X_b)+c(X-X_b)^2+d(X-X_b)^3$$

で処理を行うのが一般的です。

最小二乗法の関連記事

←使える数学へ戻る

回転行列では原点周りに点を回転させますが、任意の点（Ｃ_ｘ、Ｃ_ｙ）周りに回転させたい場合にはどうするのか？

これまでの知識を少し応用することで、意外と簡単に求めることができます。

まず、回転する点を回転中心座標が原点と一致するように点を移動させます。

次に移動した点を原点周りに回転移動させます。

回転移動後、点を原点から元の回転中心位置へ移動させます。

これで、任意点（Ｃ_ｘ、Ｃ_ｙ）周りに点を回転移動させることができました。

この 原点へ移動→原点周りの回転→元へ戻す の一連の処理を行列であらわすと

行列部分を整理すると

となり、任意点（Ｃ_ｘ、Ｃ_ｙ）周りに点を回転移動させる行列を求めることができます。

今回は回転移動について説明しましたが、拡大縮小についても、任意点を基点に拡大縮小する場合についても同様です。

この考え方は、二次元の座標だけでなく、三次元の場合も使えます。。

参考

←使える数学へ戻る

二次元座標(X,Y座標)の場合のアフィン変換行列についてはこちらで説明しましたが、今回は三次元座標(X,Y,Z座標)のアフィン変換となります。

三次元座標の場合、まず座標軸の定義、回転方向の定義を明確に覚えます。

この座標は右手座標系と呼ばれます。
フレミングの法則のときのように右手で親指、人差し指、中指をそれぞれ
直交するようにします。
このとき親指から順に親指がＸ軸、人差し指がＹ軸、中指がＺ軸の方向と
なります。
回転方向は電流と磁界の向きと同じように電流が軸の向き、磁界が回転方向
に相当します。（右ねじの法則と同じです。）

回転行列

三次元の回転行列の前に二次元の回転行列のおさらいです。
二次元の回転行列は以下の通りとなります。

これをベースに三次元座標の場合では、回転する軸の正の方向から原点の方向を見たときに、Ｘ軸、Ｙ軸はそれぞれ何軸に相当するのか？を考えれば、二次元座標のＸやＹの変数の置き換えで導き出すことができます。
行列変換しない軸に関しては単位行列でそのまま残します。

【Ｘ軸周りの回転】

【Ｙ軸周りの回転】

【Ｚ軸周りの回転】

拡大縮小行列

点（x, y, z）を原点に関してＸ軸方向にＳ_Ｘ倍、Ｙ軸方向にＳ_Ｙ倍、Ｚ軸方向にＳ_Z倍する行列は

平行移動行列

点（x, y, z）をＸ軸方向にＴ_Ｘ、Ｙ軸方向にＴ_Ｙ、Ｚ軸方向にＴ_Zだけ移動する行列は

補足

三次元の座標変換に関して検索すると座標変換は下記のように

行ベクトルで表記される場合もあるのですが、変換行列の値が変わるので、
混同しないようご注意下さい。
この表現はマイクロソフトがお得意で、ＤｉｒｅｃｔＸ（Ｄｉｒｅｃｔ３Ｄ）や.NETのアフィン変換でしか使われないので、特に必要の無い場合は覚えない方が無難です。

関連記事

二次元座標のアフィン変換についてはこちら↓にまとめています。

←使える数学へ戻る

本記事は、私がアフィン変換を勉強し始めた当初の記事になります。
今では、3×3行列の同次座標行列と呼ばれる行列しか用いておらず、こちらの方が断然おススメなので、下記ページを参照ください。

以下は、2×2行列を使ったアフィン変換の説明です。

回転行列

点（x, y）を原点まわりに反時計方向にθ度回転する行列は

拡大縮小行列

点（x, y）を原点に関してＸ軸方向にＳ_Ｘ倍、Ｙ軸方向にＳ_Ｙ倍する行列は

平行移動行列

点（x, y）をＸ軸方向にＴ_Ｘ、Ｙ軸方向にＴ_Ｙだけ移動する行列は

ただし、平行移動だけ行列の足し算になると、扱いにくい場合があるので３×３行列を用いて以下のように表す場合もあります。
というより、こちらを使う方が便利です。（私はこちらしか使いません。）

【回転行列】

【拡大縮小行列】

【平行移動行列】

とすることで、すべての座標変換を行列の積で扱うことができます。

この行列を同次座標行列と言います。
詳細はこちらを参照ください。

参考まで．．．

個人的には回転行列を覚えるのは苦手で、ＳｉｎとＣｏｓが逆になっりマイナスのつける位置を間違ったりしていたのですが、次のように考えることで少しは覚えやすくなりました。

下図のように

点（１，０）をθ度回転すると（Ｃｏｓθ、Ｓｉｎθ）
点（０，１）をθ度回転すると（－Ｓｉｎθ、Ｃｏｓθ）

に移動することはすぐにわかります。

このことを行列で表現すると
点（１，０）が（Ｃｏｓθ、Ｓｉｎθ）になることから

点（０，１）が（－Ｓｉｎθ、Ｃｏｓθ）になることから

という事がわかります。
これを合わせて表現すると

となり、回転行列が求まります。
この計算を何回か繰り返すと、そのうち覚えると思います。

←使える数学へ戻る

行列を使って連立方程式を解く方法を紹介します。

【計算例】

３つの方程式

を行列であらわすと

となります。この行列を逆行列を使ってＸ、Ｙ、Ｚに関して解くと

となり連立方程式を行列で解くことができます。

…というのは普通すぎて面白くないので、３点からなる円の方程式を行列で解く方法を紹介します。

中心（a、b）、半径　ｒ　の円の方程式

(X – a)² + (Y – b)² = r²

を展開して

X² + Y² – 2aX – 2bY + a² + b² – r² =  0

となり、A = 2a、B = 2b、C = r² – a² – b² とおくと、上式は

AX + BY + C  =   X² + Y²

となります。

ここで、円上の３点（X₀、Y₀）、（X₁、Y₁）、（X₂、Y₂）を代入すると

AX₀ + BY₀ + C  =  X₀² + Ｙ₀²
AX₁ + BY₁ + C  =  X₁² + Ｙ₁²
AX₂ + BY₂ + C  =  X₂² + Ｙ₂²

の３本の式が成り立ち、これを行列で表現すると、

となり、Ａ、Ｂ、Ｃに関して行列を解くと

となり、A、Ｂ、Ｃが求まることから、A = 2a、B = 2b、C = r² – a² – b² より
中心と半径を求めることができます。

行列を解く部分は

ガウスの消去法

のページを参照下さい。

←使える数学へ戻る