次の左グラフを見てください。こういうばらついたデータ(サンプルデータはこちら)を、右のグラフのように複数の直線(=折れ線)でフィッティングしたいと思うことありませんか?

ig_piecewise_1.png

今回は、領域ごとに異なる関数でフィッティングする方法(Piecewise Fitting)をご紹介します。


▼ Piecewise Fitting

1つのデータセットをいくつかの領域に分けて、それぞれ異なる関数、あるいは同じ関数でも異なるパラメータでフィッティングする方法を、Piecewise Fitting と呼びます。データレンジをあらかじめいくつかの領域に分けて、その領域ごとに個別にフィッティングする方法との大きな違いは、

  • フィッティングが1回で済むこと
  • 関数が連続になること*1
  • 領域の定義もフィッティングで得られること

です。

例えば、こんな例を考えてみましょう。

"温度をコントロールしながら、ある反応を観察しました。この反応は、温度を上昇させてもしばらくは大きな変化が見られず、傾きの小さな直線で表現できました(Linearな変化)。しかし、ある温度 TC を超えると、途端に劇的な変化(Exponentialな変化)がみられるようになりました。"

2つの領域のそれぞれのカーブ(Linear と Exponential)のパラメータと、この温度 TC とを同時に定めるような解析に、今回ご紹介する Piecewise Fitting は適しています*2

ig_piecewise_2.png

フィッティングですから、領域ごとの関数形はあらかじめ与える必要がありますが、このフィッティングでポイントとなるのは、x の値によって関数を切り替えることです。x が x1 < x < x2 の時は関数1、x2 を超えると関数2 に切り替える処理、いわば関数定義に if ~ else のような条件分岐する処理が可能であれば、フィッティングのできるグラフソフトで Piecewise Fitting は実現できるはずです。

さっそくグラフソフトで試してみたいところですが、冒頭のグラフを4つの直線でフィッティングできるように工夫するために、内分点についておさらいしておきます。


▼ [おさらい] 内分点

平面上の2点 A(xA, yA)と B(xB, yB) を m:n に内分する点 P(xP, yP)の座標は次の図の式1 と 2 で与えられます。

ig_piecewise_3.png

さらに、m と n はそれぞれ、xA、xB、xP を使うと、

m = xP - xA (3)
n = xB - xP (4)       (ただし、xA < xP < xB

と表現できますので、yP は領域を与える 2点A、B の座標をパラメータとする関数で定義することが可能になります。

yP = {(xB - xP)*yA + (xP - xA)*yB}/(xB - xA) (5)

非常に使いやすい形に整理できました。あと一歩です。今回のデータは、領域(線分)が4つありますので、領域の端点Q は 5個あります。端点 Q1~Q5 の x、y 座標を関数のパラメータとしますので、パラメータは全部で10個です。このパラメータを使って、4つの領域の関数 i~iv はそれぞれ次のように与えることができます。

ig_piecewise_4.png

x1 < x < x2     y = {(x2 - x)*y1 + (x - x1)*y2}/(x2 - x1) (i)
x2 < x < x3     y = {(x3 - x)*y2 + (x - x2)*y3}/(x3 - x2) (ii)
x3 < x < x4     y = {(x4 - x)*y3 + (x - x3)*y4}/(x4 - x3) (iii)
x4 < x < x5     y = {(x5 - x)*y4 + (x - x4)*y5}/(x5 - x4) (iv)

さらに、x データのレンジの外側にはフィッティングできませんので、x1 は x データの最小値(今回のデータでは1.5)、x5 は x データの最大値(今回のデータでは150)に固定することができます。この定義を、グラフソフトごとに作ってあげると、きっとうまくいくはず。ではまず、IGOR Pro で試してみます。


▼ IGOR Pro の場合

IGOR Pro のユーザ定義関数の使い方については、別の記事(ミカエリスメンテン式EC50)をご覧いただくことにして、ポイントとなる、ユーザ定義関数の中で if 文の利用について考えます。

IGOR Pro のユーザ定義関数は、ウィンドウメニュー>プロシージャウィンドウ で表示されるデフォルトとプロシージャウィンドウに関数として定義されます。この関数の中で if 文を使うだけですから、

Function Segment4(w, x) : FitFunc
	WAVE w
	Variable x

	//x1 = w[0] = hold, y1 = w[1]
	//x2 = w[2], y2 = w[3]
	//x3 = w[4], y3 = w[5]
	//x4 = w[6], y4 = w[7]
	//x5 = w[8] = hold, y5 = w[9]	
	
	Variable yy
	
	If (x < w[2])
		yy = ((w[2]-x)*w[1] + (x-w[0])*w[3])/(w[2]-w[0])
	ElseIf (x < w[4])
		yy = ((w[4]-x)*w[3] + (x-w[2])*w[5])/(w[4]-w[2])	
	ElseIf (x < w[6])
		yy = ((w[6]-x)*w[5] + (x-w[4])*w[7])/(w[6]-w[4])	
	Else
		yy = ((w[8]-x)*w[7] + (x-w[6])*w[9])/(w[8]-w[6])	
	EndIf
	
	Return yy
		
End

と、定義してあげると OK です。10個のパラメータを係数ウェーブ w[0]~w[9] に順番に割り当てます。このうち、x データのレンジから x1 と x5 はあらかじめ固定(hold)できますので、これはフィッティングを実行する際に、次のように制約条件(/H="1000000010")として指定する必要があります。

FuncFit/H="1000000010"/NTHR=0 Segment4 coef waveY /X=waveX /D

実行結果は次のようになります。

ig_piecewise_5.png

うまくいきました。


▼ KaleidaGraph

KaleidaGraph のユーザ定義関数の使い方は、こちらをご覧いただくことにして、定義する関数自体は次のようになります。

if (x<m2) ((m2-x)*m1 + (x-1.5)*m3)/(m2-1.5)
    else (if (x<m4) ((m4-x)*m3 + (x-m2)*m5)/(m4-m2) 
       else (if (x<m6) ((m6-x)*m5 + (x-m4)*m7)/(m6-m4)
          else ((150-x)*m7 + (x-m6)*m8)/(150-m6)));
            m1=12;m2=40;m3=20;m4=70;m5=90;m6=120;m7=70;m8=30

式中の 1.5 と 150 という数字は、x データのレンジに対応しています。KaleidaGraph で利用できるパラメータの数は9個(m1~m9)に限られているため、固定されることが明白な x データの両端はあらかじめ関数に記入してしまうことで、パラメータを8個まで減らしています。セミコロン ; で区切って書かれている値はフィッティングの初期値です。

ig_piecewise_6.png


▼ SigmaPlot

最後に取り上げる SigmaPlot には、なんと、2~5個の領域の直線からなる関数があらかじめ定義されているので、関数を定義する必要はありません!Regression Wizard (F5) を開いて、Equation Category のリストの末尾にある Piecewise を選択してください。関数のリスト(Equation Name)に2~5 segment linear という関数が表示されますのでこれを利用してください。

ig_piecewise_7.png

実行結果は次のようになります。

ig_piecewise_8.png


▼ まとめ

今回は IGOR Pro だけでなく、3つのソフトで Piecewise Fitting をやってみました。直線に限ると、SigmaPlot が簡単ですが、直線以外の関数による Piecewise Fitting の場合には、ソフトによらずユーザ定義関数としての定義が必要になります。今回紹介した Piecewise Fitting 以外にも、Global Fitting とか、Constrained Fitting(制約条件付きのフィッティング)とか、ちょっと変わったフィッティング方法がありますので、今後紹介していきたいと思います。

井上@技術部でした。

*1 もちろん、ここでいう連続とは、右からと左からの極限が一致するとかいうような、数学的に関数が連続という意味ではなくて、単純に、線が途切れず結ばれるという程度の意味です。
*2 Tc が物理的にどんな意味があるのかということは、カーブフィッティングではなく何か熱力学的な根拠を検討する必要があります。カーブフィッティングは、あるデータを与えた関数で一番うまく表現するパラメータを教えてくれる手法であって、物理的な現象を解き明かすための手段ではありません。

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

トラックバック

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

HULINKS サイトの新着情報