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 になった。

Reference

Example