Optimization Toolboxのlsqnonlin関数とMATLABのfminsearch関数、バックスラッシュ演算(行列の左除算「\」)を使用する方法があります。
1.lsqnonlin関数によるフィッティング(Optimization Toolbox)
円の中心座標を(cx,cy)、半径をrとし、円の方程式 (x-cx)^2+(y-cy)^2=r^2 にフィットするcx、cy、rを最小二乗法により求めることができます。
% サンプルデータの作成
t = (0:pi/100:2*pi)';
x = 10*rand(1) + 4*cos(t) + rand(size(t));
y = 10*rand(1) + 4*sin(t) + rand(size(t));
% 初期値を指定 a0=[xの初期値, yの初期値, 半径の初期値]
a0 = [mean(x),mean(y),max(x)-mean(x)];
% 円の方程式にフィッティングする係数a(1)、a(2)、a(3)を求める
f = @(a) (x-a(1)).^2 + (y-a(2)).^2 - a(3).^2;
% lsqnonlin関数によるフィッティング
af1 = lsqnonlin(f, a0);
% 円の方程式
eq = ['(x-',num2str(af1(1)),')^2+(y-',num2str(af1(2)),')^2=',num2str(af1(3)),'^2'];
% 円の描画
figure
ezplot(eq,[min([af1(1)-af1(3),af1(2)-af1(3)]),max([af1(1)+af1(3),af1(2)+af1(3)])])
hold on
% 中心点と元データの描画
plot(x,y,'*r',af1(1),af1(2),'.b')
xlim([af1(1)-af1(3)-af1(3)/10, af1(1)+af1(3)+af1(3)/10])
ylim([af1(2)-af1(3)-af1(3)/10, af1(2)+af1(3)+af1(3)/10])
axis equal
hold off
2.fminsearch関数によるフィッティング(MATLAB)
Optimization Toolboxがない場合は、目的関数（円の方程式）を二乗和の形式にしてfminsearch関数を使用することで、1.と同様に最小二乗問題として取り扱うことができます。
% サンプルデータの作成
t = (0:pi/100:2*pi)';
x = 10*rand(1) + 4*cos(t) + rand(size(t));
y = 10*rand(1) + 4*sin(t) + rand(size(t));
% 初期値を指定 a0=[xの初期値, yの初期値, 半径の初期値]
a0 = [mean(x),mean(y),max(x)-mean(x)];
% 円の方程式にフィッティングする係数a(1)、a(2)、a(3)を求める
f = @(a) norm((x-a(1)).^2 + (y-a(2)).^2- a(3).^2);
% fminsearch関数によるフィッティング
af2 = fminsearch(f, a0);
% 円の方程式
eq = ['(x-',num2str(af2(1)),')^2+(y-',num2str(af2(2)),')^2=',num2str(af2(3)),'^2'];
% 円の描画
figure
ezplot(eq,[min([af2(1)-af2(3),af2(2)-af2(3)]),max([af2(1)+af2(3),af2(2)+af2(3)])])
hold on
% 中心点と元データの描画
plot(x,y,'*r',af2(1),af2(2),'.b')
xlim([af2(1)-af2(3)-af2(3)/10, af2(1)+af2(3)+af2(3)/10])
ylim([af2(2)-af2(3)-af2(3)/10, af2(2)+af2(3)+af2(3)/10])
axis equal
hold off
3.バックスラッシュ演算によるフィッティング(MATLAB)
任意データのラジアン情報が既知の場合、つまり、Xデータ：cx+r*cos(t)、Yデータ：cy+r*sin(t) とした時の t の値が既知の場合、これらの値にフィットするcx、cy、rを最小二乗法により求めることができます。バックスラッシュ演算時には、次の連立方程式を解いています。
% サンプルデータの作成
t = [0:pi/100:2*pi]';
x = 10*rand(1) + 4*cos(t) + rand(size(t));
y = 10*rand(1) + 4*sin(t) + rand(size(t));
% 円のXデータ：a(1)+a(3)*cos(t)、Yデータ：a(2)+a(3)*sin(t)に
% フィッティングする係数a(1)、a(2)、a(3)を求める
tsize = size(t);
W = [ones(tsize), zeros(tsize), cos(t); zeros(tsize), ones(tsize), sin(t)];
% バックスラッシュ演算によるフィッティング
af3 = W \ [x;y];
figure
% 元データの描画
plot(x,y,'*r');
hold on
% 中心点の描画
plot(af3(1), af3(2), '.')
% 円の描画
plot(af3(1) + af3(3)*cos(t), af3(2) + af3(3)*sin(t))
axis equal
hold off
本件に関連するヘルプドキュメントが下記 URL よりご覧いただけます。
●バックスラッシュ演算
