TMinuitでフィッティングをする

#TMinuitとは
欧州合同原子核研究機構(CERN)の物理学者Fred Jamesによって、最小化プログラムMinuitが開発された(このときはFortranで作られた)。F. JamesによるMinuitの解説は[こちら](https://root.cern.ch/sites/d35c7d8c.web.cern.ch/files/minuit.pdf)。
このMinuitをC++を用いてROOTに実装したものがTMinuitである。

最小化プログラムは何らかの関数の最小点(厳密には極小点)を見つけるためのものであり、こちらの記事
『[TMinuitの使い方](https://qiita.com/pn11/items/36af3e3067f5d11f18eb)』では2次関数の最小点を求める方法が紹介されている。

本記事ではTMinuitのより進んだ応用として、最小二乗法に基づいた、関数のデータへの当てはめ(フィッティング)をする方法を述べる。

#最小二乗法
実験によって、次のようなn個のデータが得られたとする:

```math
(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$を

```math
\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次関数の一般的な形:

```math
f(x) = [0]*(x-[1])^2 + [2]
```
でフィッテイングする。ここで、$[i]$はパラメータである。フィットがうまくいけば、

```math
[0] = 1.0,\ \ [1]=0.0,\ \ [2]=0.0
```
に近い結果が得られるはずだ。

#マクロコード
簡単のために、マクロ内にデータを埋め込んでいる。大きいデータを扱うならcsvファイルなどから読み込むよう改変してください。

```c++:data_fitting.C
//データを格納するための配列
//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次関数とおいたことによる。

#補足

* 上のマクロコードでは、最小化アルゴリズムにMIGRADと呼ばれるものを用いた。このアルゴリズムは、[準ニュートン法](https://ja.wikipedia.org/wiki/%E6%BA%96%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95)の一種:[DFP法](https://en.wikipedia.org/wiki/Davidon%E2%80%93Fletcher%E2%80%93Powell_formula)を用いている。

* SET ERRの部分```
    arglist[0] = 3.53;//1パラメータ:1.0, 2パラメータ:2.30```
については、[Minuitの解説](https://root.cern.ch/sites/d35c7d8c.web.cern.ch/files/minuit.pdf)のp46あたりを参照。