logo
生物科学研究所 井口研究室
Laboratory of Biology, Okaya, Nagano, Japan
Home

円の最小2乗近似と描画にRのoptimxとdraw.circle関数

井口豊(生物科学研究所,長野県岡谷市)
最終更新:2023 年 12 月 29 日

ウェブ上を見ていて,意外に解説が少ないのが,円の最小2乗近似の方法であった。以下のページには, Excel を使った方法が述べられている。

しかし,これらは,「中心から各点への距離の2乗」と「円の半径の2乗」の差を最小化しようとしているため,実態に合わない円の近似となってしまう。

ここでは,実際に 8 個の点
(3, 0), (5, 0), (0, 3), (0, 5), (-3, 0), (-5, 0), (0, -3), (0, -5)
に対する円の近似を統計解析ソフト R を使い,まず直感的に考えてみよう。

x<- c(3, 5, 0, 0, -3, -5, 0, 0)
y<- c(0, 0, 3, 5, 0, 0, -3, -5)
plot(x, y, asp=1, col="magenta")

あとで近似円を描きやすいように,アスペクト比(asp)を 1 にしてある。

統計解析ソフト R による円近似用のデータプロット

この散布図ならば,中心 (0, 0),半径 4 の円が描けると推定される。

ここで重要なことは,通常最小二乗法(Ordinary Least Square: OLS)は,y方向の値しか最小化しないため,円の近似には妥当でない,ということである。この通常最小 2 乗法の計算については,実は,大学の研究者でも知らない人がいるので注意が必要である。この点に関しては,私のウェブ解説も参照: 回帰と相関,知っているようで知らない,その本質

R を使って円を描く場合は,パッケージ plotrix の中の draw.circle 関数やパッケージ shape の中の plotcircle 関数を使えば,円が簡単に描画できる。

ここでは,draw.circle 関数を使う。 plot 関数から,再度書くと,以下のようになる。

plot(x, y, asp=1, col="magenta")
library(plotrix)
draw.circle(0, 0, 4) # 中心 (0, 0),半径 4 の円

実に簡単,1行で円を描くスクリプトが書ける。

統計解析ソフト R による円の描画

その上で,中心から各点までの距離と半径の差を最小 2 乗化するために,パッケージ optimx の中の optimx 関数を利用する。この関数の利点は,オプション all.methods=T を指定すると,多くの近似法の結果を返してくれる点である。それを適合度が高い(残差が小さい)順に並べれば良い。初期値として,中心 (1, 1),半径 5 として計算したのが次のスクリプト。

f<- function(arg) {
  a<- arg[1]
  b<- arg[2]
  r<- arg[3]
  t<- (sqrt((x-a)^2 + (y-b)^2)-r)^2
  return(sum(t))
}

library(optimx)

res<- optimx(
  c(1, 1, 5), # 初期値,中心 (1, 1),半径 5
  f, control=list(all.methods=T) # 全ての近似法
)

summary(res, order=value) # 残差が小さい順に示す

必要部分だけ,結果を示す。

統計解析ソフト R による円の最小2乗近似結果

L-BFGS-B法がベストな方法で,中心 (0, 0),半径 4 の円で近似されるだろうと推定される。 初期値,中心 (0, 0),半径 4 として,再度計算させる。

f<- function(arg) {
  a<- arg[1]
  b<- arg[2]
  r<- arg[3]
  t<- (sqrt((x-a)^2 + (y-b)^2)-r)^2
  return(sum(t))
}

res<- optimx(
  c(0, 0, 4), # 初期値,中心 (0, 0),半径 4
  f, control=list(all.methods=T) # 全ての近似法
)

summary(res, order=value) # 残差が小さい順に示す

統計解析ソフト R による円の最小2乗近似結果

今度は,すっきりした結果が出た。

Home