Python3で非線形回帰

はじめに

円形アーチの計測結果から，円の中心および半径を逆算する必要性が生じたのですが，これには Scipy の非線形回帰が使えそうなので，早速やってみることにしました．

回帰したい円の方程式は

$$
x^2+(y-y_0)^2=r^2
$$

の形であり，未知数は，円中心の $y$ 座標$y_0$ と 半径 $r$ です．

実行用スクリプト

非線形回帰用のPythonプログラムは py_arch.py とし，入力ファイル名は inp_0.csv とします．データ数は少ないので，スクリプト内で入力データファイルを作成しています．1列目が計測値 $x$，２列目が計測値 $y$ であり，データをExcelでくれ〜という人が出てくるだろうから csv 形式にしておきます．(これはプログラムのテスト用データで，実際のデータはもっと多いです)

cat << EOT > inp_0.csv
-4, 0.00
-3, 1.63
-2, 2.22
-1, 2.53
 0, 2.65
 1, 2.56
 2, 2.22
 3, 1.61
 4, 0.00
EOT

python3 py_arch.py

Pythonのプログラム

scipy.optimize.leastsq を用いて，円中心のy座標および半径を求めます．
非線形回帰計算の主要部分は以下の部分です．

def func(param,x,y):
    y0 = param[0]
    rr=  param[1]
    return x**2+(y-y0)**2-rr**2

param0 = [-1.0, 1.0]
result = optimize.leastsq(func,param0,args=(x,y))
y0,rr=result[0][:]

def func は評価される関数であり，イコールゼロ ($=0$) となるべき関数（数式）を書き込みます．ここでは $x^2+(y-y_0)^2-r^2$ です．

param0 には求めるべき未知数の初期値ベクトルを入力します．
ここでは，$(y_0)_{initial}=-1.0$ および $(r)_{initial}=1.0$ としています．

argsには，求めるべき未知数ベクトル以外の引数名を与えます．ここでは計測値である，$x$ と $y$ で，args=(x,y) としています．

以下に実際に用いたPythonプログラムを示します．ついでに matplotlib による作図も行っています．

py_arch.py

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

# Function for Least squares method with scipy.optimize
def func(param,x,y):
    y0 = param[0]
    rr=  param[1]
    return x**2+(y-y0)**2-rr**2


# Main routine
fnameR='inp_0.csv'
fnameF='fig_0.png'

data = np.loadtxt(fnameR, delimiter=',')
x = data[:,0]
y = data[:,1]

param0 = [-1.0, 1.0]
result = optimize.leastsq(func,param0,args=(x,y))
y0,rr=result[0][:]
print(y0,rr,rr+y0)

# Data for Plotting
xx=np.arange(-4, 4.1, 0.1)
yy=y0+np.sqrt(rr**2-xx**2)
str0='$x^2+(y-y_0)^2=r^2$'
str1='$y_0='+'{0:.3f}'.format(y0)+'$'
str2='$r='+'{0:.3f}'.format(rr)+'$'

# PLotting
plt.plot(x,y,'bx', label='Data')
plt.plot(xx,yy,'k-', label='fitted curve', linewidth=5, alpha=0.3)
plt.text(0,1.5, str0, ha = 'center', va = 'center', fontsize=16) 
plt.text(0,1.0, str1, ha = 'center', va = 'center', fontsize=16) 
plt.text(0,0.5, str2, ha = 'center', va = 'center', fontsize=16) 
plt.xlim(-5,5)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True, numpoints=1)
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
plt.savefig(fnameF)
plt.show()

出力結果

以上