
工具顕微鏡で測定した測定点の座標から、エクセルにて円の方程式を最小二乗法で求める方法をお教え下さい。
過去の質問から、「楕円」についてのご回答があり、参照させていただき、自分なりに応用(xa+yb+c=-(x^2+y^2)として)してみたのですが、測定点の座標から得られる行列P、及び、その転置行列P'との積の逆行列とP'と()内の式の右辺から得られる行列との積を計算することができません。(3点のデータでは計算できましたが、Pを入力し、P'を求め、…と云った段階的な計算方法を採りました。)
宜しくお願いします。
No.3ベストアンサー
- 回答日時:
【1】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)」と入力するだけです。
【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。)
No.2
- 回答日時:
直交座標で表された点の集合 (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で解く方法が良く分からずに苦労しております。
また、相関係数の導き方もお教え願えればとお願いする次第であります。
No.1
- 回答日時:
最小二乗法ではありませんが、次のようにすれば、頭をまったく使わず
目的を達成することができると思います。
まず、[ツール]->[アドイン]でソルバーを入れる
(あらかじめ、入っていれば、この操作は必要ありません。
ソルバーは非線形最適化問題を解くための道具です)
考え方としては、近似円の式と測定点との誤差の2乗
の総和をソルバー機能をつかって最小化させます。
白紙のシートに、次のようにセルに入力します
カンマごとにセルを変えて入力してください。
(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がわかります。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
このQ&Aを見た人はこんなQ&Aも見ています
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
数学の思考プロセスを理解する...
-
正規分布は一見、円と何も関係...
-
大学数学 測度の問題について
-
この余りが1、余りが3という...
-
高校数学 ベクトルの計算
-
式の展開
-
2m=8はわかるのですが、2n=6...
-
4は素数じゃないですよね? こ...
-
(0,1)=[0,1]?
-
【問題】 2次関数 f(x)=x^2−2ax...
-
【問題】 f(x) = x^2 - 4a x + ...
-
コピーしたい本のページ数
-
三角形の面積は、底辺✕高さ÷2 ...
-
暗号を解除(復号)できたとい...
-
計算で優劣でますかね? コップ...
-
決定性有限オートマトン
-
積分で絶対値が中にあるときっ...
-
n!=m^2-1
-
この180➗204の計算の仕方教えて...
-
数ⅱ等式の証明について。 条件...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
至急 a²b+a-b-1 の因数分解...
-
limn→∞、10∧n=0?
-
コピーしたい本のページ数
-
ルービックキューブと群論
-
この問題、解き方は理解したの...
-
三角形の面積は、底辺✕高さ÷2 ...
-
高校数学について
-
上が✖で下が〇になる理由が、何...
-
3つの無理数a,b,cでf(x)=x^3+ax...
-
文字置き 必要条件・十分条件に...
-
(0,1)=[0,1]?
-
数学の問題点を尋ねることがで...
-
写真は2変数関数の合成微分の公...
-
【問題】 f(x) = x^2 - 4a x + ...
-
1/(s(s^2+2s+5))を部分分数分解...
-
https://youtube.com/shorts/Kw...
-
青の吹き出しの何をどう考えれ...
-
数学の質問:関数の書き方
-
数ⅱ等式の証明について。 条件...
-
ランダウの記号のとある演算
おすすめ情報