TMinuitとは
欧州合同原子核研究機構(CERN)の物理学者Fred Jamesによって、最小化プログラムMinuitが開発された(このときはFortranで作られた)。F. JamesによるMinuitの解説はこちら。
このMinuitをC++を用いてROOTに実装したものがTMinuitである。
最小化プログラムは何らかの関数の最小点(厳密には極小点)を見つけるためのものであり、こちらの記事
『TMinuitの使い方』では2次関数の最小点を求める方法が紹介されている。
本記事ではTMinuitのより進んだ応用として、最小二乗法に基づいた、関数のデータへの当てはめ(フィッティング)をする方法を述べる。
最小二乗法
実験によって、次のようなn個のデータが得られたとする:
(x_1,y_1,\sigma_1),\ (x_2,y_2, \sigma_2),\ \cdots,\ (x_n,y_n,\sigma_n)
ここで$\sigma_i$は、$y_i$についての誤差であるとする。最小二乗法とは、このデータの分布$(x,y)$が、$y=f(x;\theta)$という関数(モデル)に従うとしたとき、関数のパラメータ$\theta$を
\chi^2 = \sum^n_{i}\left\{\frac{y_i-f(x_i;\theta)}{\sigma_i}\right\}^2
が最小になるように定める方法のことである。
最小二乗法を実行するために、関数$\chi^2$の最小値と、そのときのパラメータ$\theta$の値をTMinuitで計算してみよう。
サンプルデータ
ここでは、次のようなデータがあるとしよう:
| i | x | y | $\sigma$ |
|---|---|---|---|
| 1 | -5 | 25 | 2.5 |
| 2 | -4 | 16 | 1.6 |
| 3 | -3 | 9 | 0.9 |
| 4 | -2 | 4 | 0.4 |
| 5 | -1 | 1 | 0.1 |
| 6 | 0 | 0 | 0.1 |
| 7 | 1 | 1 | 0.1 |
| 8 | 2 | 4 | 0.4 |
| 9 | 3 | 9 | 0.9 |
| 10 | 4 | 16 | 1.6 |
| 11 | 5 | 25 | 2.5 |
$y$は$x$を2乗しただけで、それに適当なエラーをつけている。
これを、2次関数の一般的な形:
f(x) = [0]*(x-[1])^2 + [2]
でフィッテイングする。ここで、$[i]$はパラメータである。フィットがうまくいけば、
[0] = 1.0,\ \ [1]=0.0,\ \ [2]=0.0
に近い結果が得られるはずだ。
マクロコード
簡単のために、マクロ内にデータを埋め込んでいる。大きいデータを扱うならcsvファイルなどから読み込むよう改変してください。
//データを格納するための配列
//1列目:変数, 2列目:値, 3列目:error
double data[11][3] = {
-5.0, 25.0, 2.5,
-4.0, 16.0, 1.6,
-3.0, 9.0, 0.9,
-2.0, 4.0, 0.4,
-1.0, 1.0, 0.1,
0.0, 0.0, 0.1,
1.0, 1.0, 0.1,
2.0, 4.0, 0.4,
3.0, 9.0, 0.9,
4.0, 16.0, 1.6,
5.0, 25.0, 2.5
};
//最小化する関数
void chi2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
//カイ2乗値を入れる変数
//chi2()が呼び出されるたび、0に初期化する
double chisq=0;
for(int i=0; i<11; i++){
double func = 0;
func = par[0]*pow((data[i][0]-par[1]),2) + par[2];//2次関数のモデル：p0*(x-p1)**2 + p2
chisq += pow((data[i][1]-func)/data[i][2],2);//カイ2乗値の計算
}
f = chisq;
}
int data_fitting(){
//TMinuitクラスをインスタンス化
TMinuit *min = new TMinuit(3);//TMinuit(n) n:パラメータ数
min->SetPrintLevel(1);
min->SetFCN(chi2);//最小化する関数をchi2に設定する
int ierflg = 0;//関数の実行結果を教えてくれる変数(0=通常通りの動作)
//メンバ関数mnparmを用いてパラメータの初期化・微分ステップの設定を行う
//初期値を設定
double vstart[3];
vstart[0] = 1.1;
vstart[1] = 0.1;
vstart[2] = 0.1;
//ステップ幅
double step[3];
step[0] = 0.01;
step[1] = 0.01;
step[2] = 0.01;
//引数 : パラメータ番号, パラメータ名, 初期値, ステップ, パラメータの上限, 同下限
//ierflg=0 if no problems, >0 if MNPARM unable to implement definition
min->mnparm(0, "p0", vstart[0], step[0], 0, 0, ierflg);
min->mnparm(1, "p1", vstart[1], step[1], 0, 0, ierflg);
min->mnparm(2, "p2", vstart[2], step[2], 0, 0, ierflg);
//最小化に必要な設定は、メンバ関数mnexcmを通して行う
//mnexcmの引数 ： コマンド名、引数が入った配列、配列内の値の数、ierflg
double arglist[10];//mnexcmのための引数を入れておく配列(適当に10個分用意)
//1σを与えるカイ2乗値の増分を設定
arglist[0] = 3.53;//1パラメータ:1.0, 2パラメータ:2.30
min->mnexcm("SET ERR", arglist, 1, ierflg);
//最小化アルゴリズムをMIGRAD(DFP法)に設定し、実行
arglist[0] = 1000;//引数1:maxcalls
arglist[1] = 1;//引数2:tolerance
min->mnexcm("MIGRAD", arglist, 2, ierflg);
//fit結果を入れておくための変数
double par0, par1, par2;
double par0_err, par1_err, par2_err;
min->GetParameter(0,par0,par0_err);
min->GetParameter(1,par1,par1_err);
min->GetParameter(2,par2,par2_err);
cout<<"**************************************"<<endl;
cout<<"p0: "<< par0 <<" +/- "<< par0_err <<endl;
cout<<"p1: "<< par1 <<" +/- "<< par1_err <<endl;
cout<<"p2: "<< par2 <<" +/- "<< par2_err <<endl;
delete min;
return 0;
}
このフィッティングの結果は
FCN=4.701e-05 FROM MIGRAD STATUS=CONVERGED 58 CALLS 59 TOTAL
EDM=9.40099e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 1.00005e+00 6.95910e-02 2.90108e-05 1.12522e-01
2 p1 -1.92964e-04 5.49035e-02 2.67968e-05 -4.51941e-01
3 p2 3.27181e-05 1.23783e-01 5.16022e-05 4.80047e-02
ERR DEF= 3.53
となって、ERRORの範囲内で上で述べた予想に一致する。
FCNは$\chi^2$関数の最小値である。これが極めて小さい値になっているのは、元になったデータが恣意的に作られたもの($y_i=x_i^2$)であり、それに沿うようにモデルを2次関数とおいたことによる。