「回帰円」を求めてみよう

Last updated at 2025-01-18Posted at 2025-01-16

1. 目的

1-1. やりたいこと

2次元平面上の3つ以上の点に対して、尤もらしい円を求めたい。
この円を「回帰円」と呼ぶことにする。

1-2. 活用例

「渦巻き状の図形の、"半径"を評価する」など。
渦巻きは、「半径が増加する円」とみなすことができる。
逆に言えば、渦巻きの半径は一定ではなく、位置によって変化する。
渦巻きの、「ある地点での半径」を求めたいときに、周辺の3点を選び、その回帰円の半径を採用することが出来る。

渦巻きの例。灰色の破線は、点295,296,297,298,299,300の回帰円。半径は5028.01

2. 理論

円の方程式を線形回帰の形に持ち込み、最小二乗法によりパラメータを求めれば、回帰円を求めることが出来る。

2-1. 円の方程式

中心 $(x_e, y_e)$ で半径$r$ の円の方程式は
$$ (x-x_e)^2+(y-y_e)^2 = r^2 $$
である。

2-2. 線形回帰の形に変形

円の方程式を線形回帰の形にすると、
$$A\overrightarrow{c}-\overrightarrow{b}=\overrightarrow{\varepsilon}$$

となる。但し
$$ A =
\begin{pmatrix}
x_0 & y_0 & 1\\
x_1 & y_1 & 1\\
& \vdots \\
x_{n-1} & y_{n-1} & 1
\end{pmatrix},
\overrightarrow{b} =
\begin{pmatrix}
x_0^2+y_0^2\\
x_1^2+y_1^2\\
\vdots \\
x_{n-1}^2 + y_{n-1}^2
\end{pmatrix},
\overrightarrow{c}=
\begin{pmatrix}
2x_e \\
2y_e \\
-(x_e^2+y_e^2-r^2)
\end{pmatrix}
$$

であり、 $\overrightarrow{\varepsilon}$ は残差を表す。

解説

線形回帰の形とは、$A\overrightarrow{c}-\overrightarrow{b}=\overrightarrow{\varepsilon}$ のことである。
但し、

$A\overrightarrow{c}$ は理論値、 $\overrightarrow{b}$ は観測値、 $\overrightarrow{\varepsilon}$ は残差をそれぞれ表す、サンプルサイズ$n$だけ長さを持つ列ベクトルである。
$\overrightarrow{c}$はパラメータ群を表す、パラメータ数 $p$ だけ長さを持つ列ベクトルであり、
$A$ はデータを表す長さ $p$ の行ベクトルを $n$ 個積み重ねた $n$行$p$列の行列である。

この $A\overrightarrow{c}-\overrightarrow{b}=\overrightarrow{\varepsilon}$ の形に、円の方程式を無理矢理変形すればよい。

さて、円の方程式を変形すると

$$2x_ex+2y_ey-(x_e^2+y_e^2-r^2)-(x^2+y^2)=0$$
となる。ここで、

$$\overrightarrow{c}=
\begin{pmatrix}
2x_e \\
2y_e \\
-(x_e^2+y_e^2-r^2)
\end{pmatrix},
\overrightarrow{a}=
\begin{pmatrix}
x & y & 1
\end{pmatrix},
b=x^2+y^2
$$
とすると、円の方程式は

$$\overrightarrow{a}
\overrightarrow{c}-b=0$$
となる。

ここに、点 $(x_i, y_i)$ ($i$ は$0$から$n-1$までの整数) を$(x,y)$ に代入し、回帰円との残差が $\varepsilon_i$ であったとすると

$$\overrightarrow{a}_i
\overrightarrow{c}-b_i=\varepsilon_i$$
となる。但し $\overrightarrow{a}_i =
\begin{pmatrix}
x_i & y_i & 1
\end{pmatrix},
b_i = x_i^2+y_i^2
$である。

さて、行ベクトル$\overrightarrow{a}_i$もスカラ$b_i$もスカラ$\varepsilon_i$も「1行の行列」とみなせなくもないが、これを$n$行に拡張した行列$A$や$\overrightarrow{b}$や$\overrightarrow{\varepsilon}$を考える。

$$ A =
\begin{pmatrix}
x_0 & y_0 & 1\\
x_1 & y_1 & 1\\
& \vdots \\
x_{n-1} & y_{n-1} & 1
\end{pmatrix},
\overrightarrow{b} =
\begin{pmatrix}
x_0^2+y_0^2\\
x_1^2+y_1^2\\
\vdots \\
x_{n-1}^2 + y_{n-1}^2
\end{pmatrix},
\overrightarrow{\varepsilon} =
\begin{pmatrix}
\varepsilon_0\\
\varepsilon_1\\
\vdots \\
\varepsilon_{n-1}
\end{pmatrix}
$$

としたとき、
$$A\overrightarrow{c}-\overrightarrow{b}=\overrightarrow{\varepsilon}$$
が成り立つ。

(解説ここまで)

目的は、残差 $\overrightarrow{\varepsilon}$ の各成分を最小にすることなので、
そのような $\overrightarrow{c}$、つまり $A\overrightarrow{c}-\overrightarrow{b}$ を最も $\overrightarrow{0}$ に近づける $\overrightarrow{c}$ 、
つまり$\left|A\overrightarrow{c}-\overrightarrow{b}\right|$を最小にする$\overrightarrow{c}$ を探せばよい。

2-3. 正規方程式

$\left|A\overrightarrow{c}-\overrightarrow{b}\right|$を最小にするような$\overrightarrow{c}$を最小二乗法で見つけるには、次の正規方程式を解けばよい。
$$A^TA\overrightarrow{c}=A^T\overrightarrow{b}$$

導出

そもそも最小二乗法とは、残差の2乗和

$$E = \sum^{n-1}_{i=0}\varepsilon_i^2$$

を、各パラメータ、つまり$\overrightarrow{c}$の各成分で偏微分した結果が0となることを連立方程式にして解くものだった。

この $E$ は $\left|\overrightarrow{\varepsilon}\right|^2$ そのものであるから、
$$
\begin{array}{lll}
E=\left|\overrightarrow{\varepsilon}\right|^2 &=& \left|A\overrightarrow{c}-\overrightarrow{b}\right|^2 \\
&=& \left(A\overrightarrow{c}-\overrightarrow{b}\right)^T\left(A\overrightarrow{c}-\overrightarrow{b}\right)\\
&=& \left(
\left(A\overrightarrow{c}\right)^T-{\overrightarrow{b}}^T
\right)
\left(A\overrightarrow{c}-\overrightarrow{b}\right)\\
&=& \left(
\left(A\overrightarrow{c}\right)^T-{\overrightarrow{b}}^T
\right)A\overrightarrow{c}
- \left(
\left(A\overrightarrow{c}\right)^T-{\overrightarrow{b}}^T
\right)\overrightarrow{b}\\
&=& \left(A\overrightarrow{c}\right)^TA\overrightarrow{c}
-{\overrightarrow{b}}^TA\overrightarrow{c}
-\left(
\left(A\overrightarrow{c}\right)^T\overrightarrow{b}
- {\overrightarrow{b}}^T\overrightarrow{b}
\right)\\
&=&\left|A\overrightarrow{c}\right|^2- \overrightarrow{b}\bullet A\overrightarrow{c}
-\left(
A\overrightarrow{c}\bullet\overrightarrow{b} - \left|\overrightarrow{b}\right|^2
\right)\\
&=&\left|A\overrightarrow{c}\right|^2-2\overrightarrow{b}\bullet A\overrightarrow{c}+\left|\overrightarrow{b}\right|^2\\
&=&\sum^{n-1}_{i=0}\left(
\left(\sum^{p-1}_{j=0}A_{i,j}c_j
\right)^2
\right)
-2\sum^{n-1}_{i=0}\left(
b_i\sum^{p-1}_{j=0}A_{i,j}c_j
\right)
+\sum^{n-1}_{i=0}b_i^2
\end{array}
$$

である。これを $c_k$ で偏微分すると

$$
\begin{array}{lll}
\frac{\partial E}{\partial c_k}
&=& \sum^{n-1}_{i=0}\left(
\frac{\partial}{\partial c_k}\left(\sum^{p-1}_{j=0}A_{i,j}c_j
\right)^2
\right)
-2\sum^{n-1}_{i=0}\left(
b_i \frac{\partial}{\partial c_k}\sum^{p-1}_{j=0}A_{i,j}c_j
\right)
+\sum^{n-1}_{i=0}\frac{\partial}{\partial c_k}b_i^2\\
&=&\sum^{n-1}_{i=0}\left(
\frac{\partial}{\partial c_k}\left(\sum^{p-1}_{j=0}A_{i,j}c_j
\right)^2
\right)
-2\sum^{n-1}_{i=0}\left(
b_i A_{i,k}
\right)\\
&=&\sum^{n-1}_{i=0}\left(
2\left(\sum^{p-1}_{j=0}A_{i,j}c_j\right)
\frac{\partial}{\partial c_k}\sum^{p-1}_{j=0}A_{i,j}c_j
\right)
-2\sum^{n-1}_{i=0}\left(
b_i A_{i,k}
\right)\\
&=&2\sum^{n-1}_{i=0}\left(
\left(\sum^{p-1}_{j=0}A_{i,j}c_j\right)A_{i,k}
\right)
-2\sum^{n-1}_{i=0}\left(
b_i A_{i,k}
\right)\\
&=&2\sum^{n-1}_{i=0}\left(
A^T_{k,i}\left(\sum^{p-1}_{j=0}A_{i,j}c_j\right)
\right)
-2\sum^{n-1}_{i=0}\left(
A^T_{k,i} b_i
\right)\\
&=&2\sum^{n-1}_{i=0}\left(
A^T_{k,i}\left(A\overrightarrow{c}\right)_i
\right)
-2\sum^{n-1}_{i=0}\left(
A^T_{k,i} b_i
\right)\\
&=&2\left(A^TA\overrightarrow{c}\right)_k-2\left(A^T\overrightarrow{b}\right)_k
\end{array}
$$

これが $0$ になればいいので、

$$\left(A^TA\overrightarrow{c}\right)_k-\left(A^T\overrightarrow{b}\right)_k=0$$

これがすべての$k$で成り立つわけだから

$$A^TA\overrightarrow{c}-A^T\overrightarrow{b}=\overrightarrow{0}$$

あとはこれを移項するだけで正規方程式
$$A^TA\overrightarrow{c}=A^T\overrightarrow{b}$$
が示される。

(導出ここまで)

2-4. 正規方程式を解く

正規方程式$A^TA\overrightarrow{c}=A^T\overrightarrow{b}$はガウスの消去法(掃き出し法)で$\overrightarrow{c}=...$の形に解ける。
このようにして$\overrightarrow{c}$ が求められれば、 $x_e, y_e, r$ も求められる。

ガウスの消去法の解説

ガウスの消去法とは、 $A\overrightarrow{x}=\overrightarrow{y}$ の形の、長さ$n$の列ベクトル方程式を、 $\overrightarrow{x}=...$ の形に変換するための手法である。($A$ は $n$行$n$列)
この手法は「前進消去」と「後退代入」という2つのフェーズに分けられる。

前進消去

前進消去では、次のような変形を行う。

$$
\begin{pmatrix}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
A_{1,0} & A_{1,1} & A_{1,2} & \cdots & A_{1,n-1} \\
A_{2,0} & A_{2,1} & A_{2,2} & \cdots & A_{2,n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
A_{n-1,0} & A_{n-1,1} & A_{n-1,2} & \cdots & A_{n-1,n-1}
\end{pmatrix} \overrightarrow{x}=\begin{pmatrix}
y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1}
\end{pmatrix}
$$

$1$以上$n$未満の各整数 $j$ について、
$j$行目から$0$行目の $\frac{A_{j,0}}{A_{0,0}}$ 倍を引く

$$
\left(\begin{array}{lllll}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
0
& A_{1,1}-A_{0,1}\frac{A_{1,0}}{A_{0,0}}
& A_{1,2}-A_{0,2}\frac{A_{1,0}}{A_{0,0}}
& \cdots
& A_{1,n-1}-A_{0,n-1}\frac{A_{1,0}}{A_{0,0}} \\
0
& A_{2,1}-A_{0,1}\frac{A_{2,0}}{A_{0,0}}
& A_{2,2}-A_{0,2}\frac{A_{2,0}}{A_{0,0}}
& \cdots
& A_{2,n-1}-A_{0,n-1}\frac{A_{2,0}}{A_{0,0}} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0
& A_{n-1,1}-A_{0,1}\frac{A_{n-1,0}}{A_{0,0}}
& A_{n-1,2}-A_{0,2}\frac{A_{n-1,0}}{A_{0,0}}
& \cdots
& A_{n-1,n-1}-A_{0,n-1}\frac{A_{n-1,0}}{A_{0,0}}
\end{array}\right) \overrightarrow{x}=\left(\begin{array}{l}
y_0 \\ y_1-y_0\frac{A_{1,0}}{A_{0,0}} \\ y_2-y_0\frac{A_{2,0}}{A_{0,0}} \\ \vdots \\ y_{n-1}-y_0\frac{A_{n-1,0}}{A_{0,0}}
\end{array}\right)
$$

$1$以上$n$未満の整数$j$, $0$以上$n$未満の整数$k$について、
$A^{(1)}_{j,k}=A_{j,k}-A_{0,k}\frac{A_{j,0}}{A_{0,0}}$, $y^{(1)}_j=y_j-y_0\frac{A_{j,0}}{A_{0,0}}$ と定義すると、上記の式は

$$
\begin{pmatrix}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
0 & A^{(1)}_{1,1} & A^{(1)}_{1,2} & \cdots & A^{(1)}_{1,n-1} \\
0 & A^{(1)}_{2,1} & A^{(1)}_{2,2} & \cdots & A^{(1)}_{2,n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & A^{(1)}_{n-1,1} & A^{(1)}_{n-1,2} & \cdots & A^{(1)}_{n-1,n-1}
\end{pmatrix} \overrightarrow{x}=\begin{pmatrix}
y_0 \\ y^{(1)}_1 \\ y^{(1)}_2 \\ \vdots \\ y^{(1)}_{n-1}
\end{pmatrix}
$$

となる。

$2$以上$n$未満の各整数 $j$ について、
$j$行目から$1$行目の $\frac{A^{(1)}_{j,1}}{A^{(1)}_{1,1}}$ 倍を引く

$$
\left(\begin{array}{lllll}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
0 & A^{(1)}_{1,1} & A^{(1)}_{1,2} & \cdots & A^{(1)}_{1,n-1} \\
0
& 0
& A^{(1)}_{2,2}-A^{(1)}_{1,2}\frac{A^{(1)}_{2,1}}{A^{(1)}_{1,1}}
& \cdots
& A^{(1)}_{2,n-1}-A^{(1)}_{1,n-1}\frac{A^{(1)}_{2,1}}{A^{(1)}_{1,1}} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0
& 0
& A^{(1)}_{n-1,2}-A^{(1)}_{1,2}\frac{A^{(1)}_{n-1,1}}{A^{(1)}_{1,1}}
& \cdots
& A^{(1)}_{n-1,n-1}-A^{(1)}_{1,n-1}\frac{A^{(1)}_{n-1,1}}{A^{(1)}_{1,1}}
\end{array}\right) \overrightarrow{x}=\left(\begin{array}{l}
y_0 \\ y^{(1)}_1 \\ y^{(1)}_2-y^{(1)}_{1}\frac{A^{(1)}_{2,1}}{A^{(1)}_{1,1}} \\ \vdots \\ y^{(1)}_{n-1}-y^{(1)}_{1}\frac{A^{(1)}_{n-1,1}}{A^{(1)}_{1,1}}
\end{array}\right)
$$

$2$以上$n$未満の整数$j$, $0$以上$n$未満の整数$k$について、
$A^{(2)}_{j,k}=A^{(1)}_{j,k}-A^{(1)}_{1,k}\frac{A^{(1)}_{j,1}}{A^{(1)}_{1,1}}$, $y^{(2)}_j=y^{(1)}_j-y^{(1)}_1\frac{A^{(1)}_{j,1}}{A^{(1)}_{1,1}}$ と定義すると、上記の式は

$$
\begin{pmatrix}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
0 & A^{(1)}_{1,1} & A^{(1)}_{1,2} & \cdots & A^{(1)}_{1,n-1} \\
0 & 0 & A^{(2)}_{2,2} & \cdots & A^{(2)}_{2,n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & A^{(2)}_{n-1,2} & \cdots & A^{(2)}_{n-1,n-1}
\end{pmatrix} \overrightarrow{x}=\begin{pmatrix}
y_0 \\ y^{(1)}_1 \\ y^{(2)}_2 \\ \vdots \\ y^{(2)}_{n-1}
\end{pmatrix}
$$

となる。

一般に $0$ 以上 $n$ 未満の各整数 $i$ について、
$i+1$以上$n$未満の各整数 $j$ について、
$j$行目から$i$行目の $\frac{A^{(i)}_{j,i}}{A^{(i)}_{i,i}}$ 倍を引き、
$0$以上$n$未満の整数$k$について、
$A^{(i+1)}_{j,k}=A^{(i)}_{j,k}-A^{(i)}_{i,k}\frac{A^{(i)}_{j,i}}{A^{(i)}_{i,i}}$, $y^{(i+1)}_j=y^{(i)}_j-y^{(i)}_i\frac{A^{(i)}_{j,i}}{A^{(i)}_{i,i}}$ と定義すると、
式は最終的に

$$
\begin{pmatrix}
A_{0,0} & A_{0,1} & A_{0,2} & \cdots & A_{0,n-1} \\
0 & A^{(1)}_{1,1} & A^{(1)}_{1,2} & \cdots & A^{(1)}_{1,n-1} \\
0 & 0 & A^{(2)}_{2,2} & \cdots & A^{(2)}_{2,n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & A^{(n-1)}_{n-1,n-1}
\end{pmatrix} \overrightarrow{x}=\begin{pmatrix}
y_0 \\ y^{(1)}_1 \\ y^{(2)}_2 \\ \vdots \\ y^{(n-1)}_{n-1}
\end{pmatrix}
$$

となる。

後退代入

後退代入では、次のような計算を行う。

$$
\begin{pmatrix}
A_{0,0} & \cdots & A_{0,n-3} & A_{0,n-2} & A_{0,n-1} \\
\vdots & \ddots & \vdots & \vdots & \vdots \\
0 & \cdots & A^{(n-3)}_{n-3,n-3} & A^{(n-3)}_{n-3,n-2} & A^{(n-3)}_{n-3,n-1} \\
0 & \cdots & 0 & A^{(n-2)}_{n-2,n-2} & A^{(n-2)}_{n-2,n-1} \\
0 & \cdots & 0 & 0 & A^{(n-1)}_{n-1,n-1}
\end{pmatrix} \overrightarrow{x}=\begin{pmatrix}
y_0 \\ \vdots \\ y^{(n-3)}_{n-3} \\ y^{(n-2)}_{n-2} \\ y^{(n-1)}_{n-1}
\end{pmatrix}
$$

分かりにくいので、これを

$$
\begin{pmatrix}
A_{0,0}' & \cdots & A_{0,n-3}' & A_{0,n-2}' & A_{0,n-1}' \\
\vdots & \ddots & \vdots & \vdots & \vdots \\
0 & \cdots & A_{n-3,n-3}' & A_{n-3,n-2}' & A_{n-3,n-1}' \\
0 & \cdots & 0 & A_{n-2,n-2}' & A_{n-2,n-1}' \\
0 & \cdots & 0 & 0 & A_{n-1,n-1}'
\end{pmatrix}
\begin{pmatrix}
x_0 \\
\vdots \\
x_{n-3} \\
x_{n-2} \\
x_{n-1}
\end{pmatrix}
=\begin{pmatrix}
y_0' \\ \vdots \\ y_{n-3}' \\ y_{n-2}' \\ y_{n-1}'
\end{pmatrix}
$$

と書き直す。

$n-1$行目に注目すると、
$$
\begin{array}{rlcc}
0x_0+...+0x_{n-3}+0x_{n-2}+A_{n-1,n-1}'x_{n-1}& &=& y_{n-1}'\\
A_{n-1,n-1}'x_{n-1}&+\sum_{j=n}^{n-1}\left(A_{n-1,j}'x_{j}\right) &=& y_{n-1}'
\end{array}
$$
なので、
$$x_{n-1}=\frac{y_{n-1}'-\sum_{j=n}^{n-1}\left(A_{n-1,j}'x_{j}\right)}{A_{n-1,n-1}'}$$ である。

但し、　$\sum_{j=n}^{n-1}\left(A_{n-1,j}'x_{j}\right)$ は$j$の開始 $n$ が終了 $n-1$よりも大きいため、$0$である。

$n-2$行目に注目すると、
$$
\begin{array}{rlcc}
0x_0+...+0x_{n-3}+A_{n-2,n-2}'x_{n-2}&+A_{n-2,n-1}'x_{n-1} &=& y_{n-2}'\\
A_{n-2,n-2}'x_{n-2}&+\sum_{j=n-1}^{n-1}\left(A_{n-2,j}'x_{j}\right) &=& y_{n-2}'
\end{array}
$$
なので、
$$x_{n-2}=\frac{y_{n-2}'-\sum_{j=n-1}^{n-1}\left(A_{n-2,j}'x_{j}\right)}{A_{n-2,n-2}'}$$ である。

但し、　$\sum_{j=n-1}^{n-1}\left(A_{n-2,j}'x_{j}\right)$ に出てくる $x_j$ は $x_{n-1}$ のみであり、これは先ほど求めたので、結局 $x_{n-2}$ も求められる。

一般に $i$に$n-1, n-2, n-3, ..., 0$ をこの順に代入した場合について、
$i$行目に注目すると、
$$
\begin{array}{rlcc}
A_{i,i}'x_{i}&+\sum_{j=i+1}^{n-1}\left(A_{i,j}'x_{j}\right) &=& y_{i}'
\end{array}
$$
なので、
$$x_{i}=\frac{y_{i}'-\sum_{j=i+1}^{n-1}\left(A_{i,j}'x_{j}\right)}{A_{i,i}'}$$ である。

但し、$\sum_{j=i+1}^{n-1}\left(A_{i,j}'x_{j}\right)$ に出てくる $x_j$ はすべて事前に求められている。

以上で、 $\overrightarrow{x}$ (の各成分)が求められた。

(ガウスの消去法の解説ここまで)

3. コーディング (python)

fit_circle.py

def fit_circle(points:list[tuple[float, float]])->tuple[float, float, float]:
    A:list[list[float]] = []
    b:list[float] = []
    for (x, y) in points:
        A.append([x, y, 1])
        b.append(x**2 + y**2)

    # 行列の転置と積を計算
    AT:list[list[float]] = [[A[j][i] for j in range(len(A))] for i in range(3)]
    ATA:list[list[float]] = [
        [sum(AT[i][k] * A[k][j] for k in range(len(A))) for j in range(3)]
        for i in range(3)
    ]
    ATb:list[float] = [
        sum(AT[i][k] * b[k] for k in range(len(b))) for i in range(3)
    ]

    # ガウスの消去法で解を求める
    def gauss_eliminate(A:list[list[float]], y:list[float])->list[float]:
        n:int = len(y)

        # 前進消去
        for i in range(n):
            for j in range(i + 1, n):
                ratio:float = A[j][i] / A[i][i]
                for k in range(n):
                    A[j][k] -= ratio * A[i][k]
                y[j] -= ratio * y[i]

        # 後退代入
        x:list[float] = [0.0] * n
        for i in range(n-1, -1, -1):
            x[i] = (
                y[i] - sum(A[i][j] * x[j] for j in range(i + 1, n))
            ) / A[i][i]
        return x

    c:list[float] = gauss_eliminate(ATA, ATb)
    xe:float = c[0] / 2
    ye:float = c[1] / 2
    r:float = (xe**2 + ye**2 + c[2]) ** 0.5
    return (xe, ye, r)

例0、点 $(0,0), (1,1), (2,3)$ に対する回帰円を求めたいなら、次のようにする。

python fit_circle.py

>>> points = [(0,0), (1,1), (2,3)]
>>> fit_circle(points)
(-3.4999999999999774, 4.499999999999985, 5.700877125495664)

こうして、中心 $(-3.5, 4.5)$、半径 $5.7$ の回帰円が求められた。

例1、点 $(-3,-4), (0,1), (3,5), (7,9)$ に対する回帰円を求めたいなら、次のようにする。

python fit_circle.py

>>> points = [(-3,-4), (0,1), (3,5), (7,9)]
>>> fit_circle(points)
(36.09310344827625, -23.65862068965548, 43.74630271450496)

こうして、中心 $(36.1, -23.7)$、半径 $43.7$ の回帰円が求められた。

4. 終わりに

3章のコードはまだまだ完ぺきではない。
例えば、点が直線上に並ぶ場合は、回帰円の半径が理論上無限大となり、この計算を行うと ZeroDivisionError が発生する。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up