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

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

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3922806/31e29c98-9605-2bcb-f889-75c46042d635.png)
渦巻きの例。灰色の破線は、点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}$ は残差を表す。

<details><summary>解説</summary>

線形回帰の形とは、$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}$$
が成り立つ。

(解説ここまで)

</details>

目的は、残差 $\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}$$

<details><summary>導出</summary>
そもそも最小二乗法とは、残差の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}$$
が示される。

(導出ここまで)
</details>


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

<details><summary>ガウスの消去法の解説</summary>

ガウスの消去法とは、 $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}$ (の各成分)が求められた。

(ガウスの消去法の解説ここまで)
</details>


# 3. コーディング (python)

```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$ の回帰円が求められた。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3922806/516fd99d-ec09-07a6-1352-d23867fb53e2.png)

例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$ の回帰円が求められた。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3922806/e677540c-9458-1e88-a64a-c39bf4dd0018.png)


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