円の方程式を最小二乗法で求める

工具顕微鏡で測定した測定点の座標から、エクセルにて円の方程式を最小二乗法で求める方法をお教え下さい。
　過去の質問から、「楕円」についてのご回答があり、参照させていただき、自分なりに応用（xa+yb+c=-(x^2+y^2)として）してみたのですが、測定点の座標から得られる行列P、及び、その転置行列P'との積の逆行列とP'と()内の式の右辺から得られる行列との積を計算することができません。（３点のデータでは計算できましたが、Pを入力し、P'を求め、…と云った段階的な計算方法を採りました。）
　宜しくお願いします。

No.3 ベストアンサー 回答者： stomachman 回答日時： 2003/02/12 19:14

【１】Excelの基本的機能を使って計算を行うworksheetを作る方法を説明します。

●A列にp(n)（x座標の測定値)、B列にq(n)（y座標の測定値)を並べます。説明のためにn=1～10としておきましょう。

●C列とD列にそれぞれ2p(n)と2q(n)を入れ、またE列には1を入れ、F列に(p(n)^2 + q(n)^2)を入れます。このためにC列1行目に「=2*A1」、D列1行目に「=2*B1」、E列1行目に「1」、F列1行目に「=A1^2+B1^2」と入力します。そしてC列1行目からF列10行目までを選択して、＜下方向へコピー＞をやる。

●H列1行目からK列3行目に正規方程式の係数表を作ります。
まずH列1行目に「=MMULT(TRANSPOSE(C1:E10),C1:F10)」と入力します。でもこれだけだと、エラーになりますよね。
H列1行目からK列3行目を選択すると数式バーに今入力した式「=MMULT(TRANSPOSE(C1:E10),C1:F10)」が表示されます。この式にカーソルを合わせて一度クリックし、それからCTRLキーを押しながらEnterキーを押します。（Macintoshならリンゴマークキーを押しながらenter。）
すると数式バーが「{=MMULT(TRANSPOSE(C1:E10),C1:F10)}」という表示に変わり、行列が表示されます。

●H列5行目からJ列7行目に、H列1行目からJ列3行目の3x3の行列に対する逆行列を作ります。このために、H列5行目に「=MINVERSE(H1:J3)」と入力します。そして、H列5行目からJ列7行目を選択し、数式バーに表示された「=MINVERSE(H1:J3)」をクリックして、それからCTRLキーを押しながらEnterキーを押します。（Macintoshならリンゴマークキーを押しながらenter。）
* なおここで、「H列5行目からJ列7行目を選択」する際には、まずH列5行目のセルをクリックし、次にshiftキーを押しながらJ列7行目のセルをクリックします。（先にH列5行目のセルをクリックしないと式が表示されません。）

●M列1行目からM列3行目に、a,b,cの値を計算します。M列1行目のセルに「=MMULT(H5:J7,K1:K3)」と入力し、M列1行目からM列3行目を選択し、数式バーに表示された「=MMULT(H5:J7,K1:K3)」をクリックして、それからCTRLキーを押しながらEnterキーを押します。（Macintoshならリンゴマークキーを押しながらenter。）

●M列5行目に、rを計算します。M列5行目のセルに「=SQRT(M3+M1^2+M2^2)」と入力するだけです。

【２】相関係数

二つのデータ系列の間で定義されるものですから、この問題の場合にはそのままでは馴染みません。必要なのは「フィッティングした円がどのぐらい、データと旨く合っているか」でしょう。そこで、各データ(p(n),q(n))と推定した中心(a,b)との距離と、推定した半径rとのずれを残差と考えるのが適当かと思います。すなわち
ε(n) = √((p(n)-a)^2+(q(n)-b)^2) - r
です。

●N列1行目から10行目にこの表を作ってみましょう。N列1行目に「=sqrt((A1-$M$1)^2+(B1-$M$2)^2)-$M$5」と入力し、N列1行目からN列10行目までを選択して、＜下方向へコピー＞をやる。

これがどのぐらいのばらつきであるかを見るために、平均（=average(N1:N10)）と標準偏差(=stdev(N1:N10)）を計算してみると良いでしょう。絶対値が最大のものを計算するにはセルに「=max(abs(N1:N10))」と入力して、数式バーに表示された式をクリックして、それからCTRLキーを押しながらEnterキーを押します。（Macintoshならリンゴマークキーを押しながらenter。）

この回答へのお礼

非常に解り易い説明をして戴き、誠にありがとうございます。
今後とも宜しくお願い申し上げます。

お礼日時：2003/02/15 17:29

No.2 回答者： stomachman 回答日時： 2003/02/11 15:05

直交座標で表された点の集合　(p(n),q(n))　に円をフィッティングするんですね。

円を表すモデルは
(p(n) - a)^2 + (q(n) - b)^2 = r^2
です。これを展開して
p(n)^2 - 2ap(n) + a^2 + q(n)^2 - 2bq(n) + b^2 - r^2 = 0
整理すると
(p(n)^2 + q(n)^2) = 2ap(n) + 2bq(n) - (a^2 +b^2 - r^2)
ここで、
c = -(a^2 +b^2 - r^2)
とおけば
(p(n)^2 + q(n)^2) = a(2p(n)) + b(2q(n)) + c
という、a,b,cに関する一次方程式になります。これがn=1,2,3....,N (N=測定点の個数）
だけ並んでいる連立一次方程式です。
こいつを最小二乗の意味で解く方法は、既におわかりなんでしょうか？（Excelの使い方の話に過ぎないんですが）

で、a,b,cが分かり、a,b,cからrが計算できます。

この回答への補足

ご回答ありがとうございます．
連立1次方程式をExcelで解く方法が良く分からずに苦労しております。
また、相関係数の導き方もお教え願えればとお願いする次第であります。

補足日時：2003/02/12 17:58

No.1 回答者： 0shiete 回答日時： 2003/02/10 19:41

最小二乗法ではありませんが、次のようにすれば、頭をまったく使わず

目的を達成することができると思います。

まず、[ツール]->[アドイン]でソルバーを入れる
（あらかじめ、入っていれば、この操作は必要ありません。
ソルバーは非線形最適化問題を解くための道具です）

考え方としては、近似円の式と測定点との誤差の２乗
の総和をソルバー機能をつかって最小化させます。

白紙のシートに、次のようにセルに入力します
カンマごとにセルを変えて入力してください。

（xi,yi）は測定点の座標、
(Cx,Cy)は近似円の中心座標（初期値を適当にめぼしをつけて入れてください）
Rは近似円の半径（初期値を適当にめぼしをつけて入れてください）


-------------------------------------------------
x1,y1,Cx,Cy,R,=((a1-c1)^2+(b1-d1)^2-e1^2)^2
x2,y2,=c1,=d1,=e1,=((a2-c2)^2+(b2-d2)^2-e2^2)^2
x3,y3,=c2,=d2,=e2,=((a3-c3)^2+(b3-d3)^2-e3^2)^2
x3,y3,=c3,=d3,=e3,=((a4-c4)^2+(b4-d4)^2-e4^2)^2
......
xn,yn,=c(n-1),=d(n-1),=e1,=((an-cn)^2+(bn-dn)^2-e1^2)^2
=sum(f1:fn)
--------------------------------------------------


これでソルバーをつかって「=sum(f1:fn)」の入れてあるセルを
目的とするセルとして選択し、c1,d1,e1のセルを変化させるセルとして
選択し、最小化すれば、Cx,Cy,Rがわかります。