* 本記事はスクリプトや操作関数に関する記載が多くなっています。あらかじめご了承ください。

最近時々相談を受ける、陰関数でのフィッティングについてご紹介します。陰関数でフィッティング、というのは、もう少し言葉を補うと、陰関数表示された関数でフィッティングをするという意味合いです。

通常の、つまり、陽関数表示された関数というのは、

y = f(x)

と、y= なんとか、という形、言い換えれば、関数を y について解いた表現形式を意味します。一方、陰関数表示というのは、

f(y,x) = 0

と表現される関数で、y= なんとか、という形では表現できない、あるいは表現していない関数という意味です。数学的な詳細は、専門書をご覧ください。

ここでは、IGOR Pro を使って、陰関数の代表格ともいえる、楕円の方程式でフィッティングしてみます。

▼ おさらい:楕円の方程式

陰関数の例として、ここではシンプルに原点を中心とする楕円の方程式を考えます。

ellipse.png

a > b なら、上下につぶれた楕円になるわけで、その焦点の座標は

focus.png

で与えられるわけです。さらに、θを媒介変数に使うと、

parameter.png

と、表現できるのでした。おー、懐かしいですね。


▼ 楕円の方程式を定義する

通常の陽関数表示のユーザ定義関数については、これまでにも触れていますように(ココ)、メニューか UI を使って定義することができますが、残念ながら、陰関数となる楕円の方程式は定義できません。次のような関数をプロシージャウィンドウに直接定義する必要があります。

Function FitEllipse(w, x, y) : FitFunc
    Wave w
    Variable x
    Variable y

    //Coefficients 2
    //w[0] = a
    //w[1] = b

    return (x/w[0])^2 + (y/w[1])^2 - 1

End

大切なのは、最後の return の行です。陽関数表示であれば、y = f(x) の右辺を定義してあげました。しかし陰関数表示では、y = f(x) になっていませんので、y - f(x) = 0 と左辺にまとめて、y - f(x) を return してあげます。

▼ フィッティングの実行

もう一つの問題は、フィッティングを解析メニューからは実行できないこと、です。IC 50 のケースのように、解析メニューから定義した関数を選択できません。コマンドラインで、

FuncFit /ODR=3 関数名, 係数ウェーブ, /X={入力ウェーブ} /XD={出力ウェーブ}

と操作する必要があります。/ODRフラグがポイントです。ODR は Orthogonal Distance Regression の略で、Y データだけでなく、X データの誤差も考慮するフィッティング方法です。/ODR フラグを使うと、こういった特殊なフィッティング方法が利用できるようになるわけです。

▼ 具体例

では、やってみましょう。コンカイモ ソウサハ ゼンブ コマンド デス。まず、データを作ってグラフに表示してみます。

Make/N=20 theta,ellipseY,ellipseX
//データ作成
theta = 2*pi*p/20
ellipseY = 2*cos(theta)
ellipseX=3*sin(theta)
//グラフ表示
Display ellipseY vs ellipseX
ModifyGraph mode=3,marker=8
ModifyGraph rgb=(0,0,65280)
ModifyGraph width={perUnit,30,bottom},height={perUnit,30,left}

20ポイントの ellipseX と ellipseY を作成し、theta を媒介変数に、データを合成しています。

data.png

これではきれいな楕円になってしまうので、gnoise() 関数を使って、ノイズを乗せます。

SetRandomSeed 0.5
ellipseY += gnoise(.2)
ellipseX += gnoise(.2)

random.png

SetRandomSeed を使うと、乱数を初期化することができ、どんな環境でも同じ結果を得ることが可能になります。

次に、出力ウェーブと係数ウェーブを作成します。

Duplicate ellipseY, ellipseYFit, ellipseXFit
Make/D/O ellipseCoefs={3,2}			// a, b

係数ウェーブを倍精度で定義することは、フィッティングポイントです。

最後に、フィッティングを実行しましょう。

FuncFit/ODR=3 FitEllipse, ellipseCoefs /X={ellipseX, ellipseY}/XD={ellipseXFit,ellipseYFit}

うまくいけば、次の出力が得られるはずです。

fitting.png

*FuncFit/ODR=3 FitEllipse, ellipseCoefs /X={ellipseX, ellipseY}/XD={ellipseXFit,ellipseYFit}
  適切に収束しました
  ellipseCoefs={3.053,1.8973}
  V_chisq= 0.622108;V_npnts= 20;V_numNaNs= 0;V_numINFs= 0;
  W_sigma={0.0829,0.0653}
  係数値 ±標準偏差
  	w_0	=3.053 ア 0.0829
  	w_1	=1.8973 ア 0.0653

媒介変数表示のところで、ellipleX の係数を 3、ellipseY の係数を 2 としていますので、良い結果が得られていると言えるでしょう。フィッティング結果をプロットすると、こんな感じになりますね。青が元のデータで、赤×がフィッティングで得られたデータです。

AppendToGraph ellipseYFit vs ellipseXFit
ModifyGraph mode=3,marker(ellipseYFit)=1

result.png

▼ まとめ(種明かし)

今回ご紹介した陰関数として楕円の方程式でフィッティングする例は、IGOR Pro のヘルプCurve Fitting.ihf の Example: Fit to an Ellipse の項をさらに簡単な例に再編集したものです。ご興味のある方は、そちらもご覧いただければと思います。

井上@技術部でした。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/467
内容に対しての関連性がみられないものは削除する場合があります

コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報