測定点の楕円近似

質問

傾斜した円筒を水平に測定して得られたi個の測定点(X[i],Y[i])を楕円の公式にフィッティングし、その楕円の中心座標(X0,Y0),楕円のなす角度θを求めたいと考えております。
とあるサイト上で楕円の式
f(x,y)={(x-x0)cosθ+(y-y0)sinθ}^2/a^2+{(x-x0)sinθ-(y-y0)cosθ}^2/b^2=1
ただし、a=長半径　b=短半径
として、多変量の非線形最小二乗法にて当てはめる事により5つの変数を求めることが出来るとあり、そのプログラムも公開されていたのですが、私、数学ならびにプログラムに対して全くの素人の為、色々とサイトや文献を調べたのですがどのように行えば良いのか理解する事が出来ませんでした。
出来ればExcelシートで計算出来ればと思っております。
考え方、公式等だけでも構いませんので、どうかご教授願えませんでしょうか？
皆様方、どうか宜しくお願いいたします。

回答No.2

stomachman · ベストアンサー

おっとと
> 『q = inv(P) r ,q = inv(P' P) (P' r) 』を解く
のではありません。もう解いてあります。

このご質問の場合、測定点は５より多いでしょうから、
inv(P' P) (P' r)
を計算すれば答が出るんです。Excelの関数で書けば
=mmult(minverse(mmult(transpose(P),P),mmult(transpose(P),r))
ということです。

例えば10個測定点があるとしましょうか。
PをA列1行目からE列10行目までの範囲に入れることにしましょう。
j番目の測定値のxの値をxj, yの値をyjとするとき、
j行目のA列にxj^2, B列にyj^2, C列にxj*yj, D列にxj,E列にyj
を並べたものです。10行5列の行列ということになります。
具体的にはD列とE列にデータを入れて、
A列1行目に =C1^2
B列1行目に =D1^2
C列1行目に =C1*D1
を入力し、
A列1行目からC列10行目までの範囲を選択しておいて、「下方向へコピー」のコマンドを使えばできあがりです。

rというのは1を縦にずらりと10個並べたもの。
これをG列1行目からG列10行目までの範囲に置くことにしましょう。

そして、答のベクトルqをJ列1行目からJ列5行目に置くことにしましょうか。
答は5個の数値が縦に並んだもので、
1個目がa, 2個目がb, ..., 5個目がeの値です。

そうすると、J列1行目からJ列5行目の範囲を選んだ状態で、J列1行目のセルに
=mmult(minverse(mmult(transpose(A1:E10),A1:E10),mmult(transpose(A1:E10),G1:G10))
と入力します。で、数式入力の欄にカーソルを入れたまま、controlキーを押しながらEnterを押す。(Macならリンゴマークのキーを押しながらEnterを押す。)すると、数式入力の欄には
{=mmult(minverse(mmult(transpose(A1:E10),A1:E10),mmult(transpose(A1:E10),G1:G10))}
と表示され、答がJ列1行目からJ列5行目に表示されます。

Excelにおいて、行列全体にPだのrだのと言った名前をつける方法も知っていると便利です。Excelのマニュアルをご参照あれ。

お礼

お礼を申し上げるのが遅くなり、申し訳ありませんでした。ご丁寧なご回答本当にありがとうございました。

回答No.1

stomachman

しばしば非線形の問題として扱うようですが、実は線形の簡単な問題として解くことができます。

一般に、楕円や双曲線は
a(x^2)+b(y^2)+cxy+dx+ey=1
(a,b,c,d,eは実数）を満たす<x,y>の集合として表されます。
楕円あるいは双曲線の曲線上の座標<x,y>をn回測定したものを<x[j],y[j]> (j=1,2,...,n)とし、その値に誤差がないとするならば、５個の測定点があれば５元連立一次方程式
a(x[1]^2)+b(y[1]^2)+cx[1]y[1]+dx[1]+ey[1]=1
　：
a(x[5]^2)+b(y[5]^2)+cx[5]y[5]+dx[5]+ey[5]=1
を解けばa,b,c,d,eが決まります。これを行列で表せば
P[j,1] = x[j]^2
P[j,2] = y[j]^2
P[j,3] = x[j]y[j]
P[j,4] = x[j]
P[j,5] = y[j]
とし、
r[j] = 1
そして未知数のベクトルを
q[1] = a
q[2] = b
q[3] = c
q[4] = d
q[5] = e
とすれば
P q = r
を解くことに相当します。つまり
q = inv(P) r
ここにinv(P)はPの逆行列です。

測定点が５個より多い場合(n>5)にも j=1,2,...,nについてP,rを作り
q = inv(P' P) (P' r)
を計算すればq（すなわちa,b,c,d,e）が求められます。ここにP'はPの転置行列です。
nが大きいほど（つまり測定点が多いほど）測定値に少し誤差があってもその影響を受けにくくなります。

Excelでは 行列の転置はtranspose(), 逆行列はminverse()、行列（或いはベクトル）の積はmmult(,)という関数で計算できます。

a,b,c,d,eから、楕円の長軸・短軸の長さと向きおよび中心の座標を求めるのは、別の問題と考えれば良いでしょう。

お礼

早速のご回答誠にありがとうございます。
ご丁寧に解法をご説明していただいたのですが、何分レベルが高く、私、数学的知識がほとんど無いため、『q = inv(P) r ,q = inv(P' P) (P' r) 』を解くというお話から理解する事が出来ませんでした。本内容について理解すべく勉強いたします。
ありがとうございました。