ROOT の TMinuit を使おうとして検索したけど、ネット上に分かりやすい例 (特に日本語) がなかったので、人に教わったことをメモしとく。
MINUIT
MINUIT は CERN で開発された最小値を探すプログラム。元々は Fortran で書かれていたけど ROOT にも移植されていて TMinuit というクラスになっている。TMinuit は TH1::Fit とかみたいにネット上にごろごろ例が転がってるわけじゃなくとっつきにくかったけど、大して難しいわけでもなかった。
簡単な例
二次関数
x^2+x
を最小にする $x$ を求めてみる。中学校で習ったように答えは -1/2. これを MINUIT で解くプログラムは、
#include <iostream>
#include <string>
using namespace std;
#include "TMinuit.h"
void chi2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
f = par[0]*par[0] + par[0] + 1;
}
int main (){
TMinuit *min = new TMinuit(1);
min->SetPrintLevel(1);
double par, parErr;
string parName = "test_para";
double stepSize = 0.1, minVal = -100.0, maxVal = 100.0;
min->DefineParameter(0, parName.c_str(), par, stepSize, minVal, maxVal);
min->SetFCN(chi2);
int migrad_stats = min->Migrad();
min->GetParameter(0, par, parErr);
cout << "Result: " <<par << " +/- " << parErr << endl;
cout << "Status of Migrad: " << migrad_stats << endl;
delete min;
return 0;
}
TMinuit *min = new TMinuit(1); の 1 は変数の数。min->SetFCN(chi2); で最小化したい関数をセットする。chi2 という関数名にしたけど、今回は全然 $\chi^2$ ではない。
http://www-glast.slac.stanford.edu/software/root/GRUG/docs/Feature/GRUGminuit.pdf によると、
TMinuit expects the function to minimize to be ALWAYS a static external function (no member functions).
だそうなので、extern void chi2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) というふうにした。あんまり意味は分かっていないけど、大事なのは、f = [関数] となるようにすることと、par がパラメターの配列であるということ。
さらに min->DefineParameter(0, parName.c_str(), par, stepSize, minVal, maxVal); でパラメター0の名前やステップサイズ、範囲などを決めている。
PrintLevel は -1, 0, 1 があるらしい。
min->Migrad(); で最小化の計算を実行できる。このとき返り値は fitting のステータスになっている。詳細はこれ。抜粋すると、
0: command executed normally
1: command is blank, ignored
2: command line unreadable, ignored
3: unknown command, ignored
4: abnormal termination (e.g., MIGRAD not converged)
5: command is a request to read PARAMETER definitions
6: 'SET INPUT' command
7: 'SET TITLE' command
8: 'SET COVAR' command
9: reserved
10: END command
11: EXIT or STOP command
12: RETURN command
結果を持ってくるには min->GetParameter(0, par, parErr); のようにする。
出力は以下のようになった。
**********
** 1 **SET PRINT 100
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 test_para 2.10781e-317 1.00000e-01 -1.00000e+02 1.00000e+02
**********
** 2 **MIGRAD
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-04
FCN=1 FROM MIGRAD STATUS=INITIATE 4 CALLS 5 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 test_para 0.00000e+00 1.00000e-01 1.00000e-03 1.00000e+02
NO ERROR MATRIX
FCN=0.75 FROM MIGRAD STATUS=PROGRESS 7 CALLS 8 TOTAL
EDM=1.74834e-11 STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 test_para -4.99998e-01 1.00000e-01 -5.00000e-03 4.18132e-04
RELATIVE CHANGE IN COV. MATRIX= 50.0 per cent
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
START COVARIANCE MATRIX CALCULATION.
EIGENVALUES OF SECOND-DERIVATIVE MATRIX:
1.0000e+00
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=0.75 FROM MIGRAD STATUS=CONVERGED 12 CALLS 13 TOTAL
EDM=8.74186e-12 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 test_para -4.99998e-01 9.99983e-01 6.45935e-06 4.18130e-04
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 1 ERR DEF=1
1.000e+00
Result: -0.499998 +/- 0.999983
Status of Migrad: 0
よく分かってないけど、答えが大体 -0.5 になった。