IntegralError TF2

p_x · May 16, 2014, 11:02am

Dear all,

I am fitting a scatter plot using the template $ROOTSYS/tutorials/fit/fit2.C. The function is:

// --------
Double_t g2(Double_t x, Double_t par) {
Double_t r1 = Double_t((x[0]-par[1])/par[2]);
Double_t r2 = Double_t((x[1]-par[3])/par[4]);
return par[0]TMath::Exp(-0.5(r1r1+r2r2));
}
Double_t fun2(Double_t *x, Double_t *par) {
Double_t *p1 = &par[0];
Double_t *p2 = &par[5];
Double_t *p3 = &par[10];
Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3);
return result;
}
TF2 *f2 = new TF2(“f2”,fun2,-10,10,-10,10, npar);
// --------

I would like to obtain the Integral and Integralerror of the three indivual functions g2 using the covariance matrix which I obtain from the fit.

I have loooked through the forum but only found examples for TF1. I suppose the process is the same:

TFitResultPtr r = h2->Fit(f2, “S”);
TMatrixDSym cov = r->GetCovarianceMatrix();

and later use: f2->Integral() and f2->IntegralError() ?

But I seem not to get the right output when trying this, probably becuase my arguments into f2->Integral() and f2->IntegralError() are not correct.
Can somebody please help me out?

thanks in advance

moneta · May 17, 2014, 2:15pm

Hi,

You need to make a covariance matrix for the parameters of each function, by using a global covariance matrix obtained from the fit. Inn case of problem, please post your macro, and I can try to help you

Lorenzo

p_x · May 19, 2014, 7:05am

Dear Lorenzo,

thank you very much for your answer. I will try that.

Best regards

p_x · May 20, 2014, 9:38am

Dear Lorenzo,

I am back with a working example however I am not sure if it is correct.
First I obtain the full covariance matrix and create individual covariance matrices for the three subfunctions. Then three TF2 (fg1,fg2 and fg3) are assigned to the three g2 functions and they get the final parameters from the fit. The integral and Integralerror are calculated by dividing with the binsize. The output from this example is shown below the code.

The total integral agrees with the data points (~1e5). The errors seem little too small but I don’t know.

Do you think it is correct?
One question is why we only include the covariance matrix from each subfunction. What about the other elements between the subfunctions, should they not add to the total uncertainty or do they cancel out?
In my own fitting code I have fixed some parameters. Is there a complication by doing so in using the example below given if it is correct? I read for instance that < TVirtualFitter *fitter = TVirtualFitter::GetFitter(); > probable does not work with fixed parameters?

Thank you very much for your time and help. Best regards


#include "TF2.h"
#include "TH2.h"
#include "TMath.h"
#include "TVirtualFitter.h"
#include "TFitResultPtr.h"
#include "TMatrixDSym.h"
#include "TMatrixD.h"
#include "TMatrix.h"
#include "TFitResult.h"
#include "TCanvas.h"


// Fitting a 2-D histogram
// This tutorial illustrates :
//  - how to create a 2-d function
//  - fill a 2-d histogram randomly from this function
//  - fit the histogram
//  - display the fitted function on top of the histogram 
//
// This example can be executed via the interpreter or ACLIC
//   root > .x fit2.C
//   root > .x fit2.C++
//Author: Rene Brun
         
Double_t g2(Double_t *x, Double_t *par) {
   Double_t r1 = Double_t((x[0]-par[1])/par[2]);
   Double_t r2 = Double_t((x[1]-par[3])/par[4]);
   return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2));
}   
Double_t fun2(Double_t *x, Double_t *par) {
   Double_t *p1 = &par[0];
   Double_t *p2 = &par[5];
   Double_t *p3 = &par[10];
   Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3);
   return result;
}

void fit2() {
   const Int_t npar = 15;
   Double_t f2params[npar] = 
      {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7};
   TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar);
   f2->SetParameters(f2params);

   //Create an histogram and fill it randomly with f2
   TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10);
   Int_t nentries = 100000;
   h2->FillRandom("f2",nentries);
   //Fit h2 with original function f2
   Float_t ratio = 4*nentries/100000;
   f2params[ 0] *= ratio;
   f2params[ 5] *= ratio;
   f2params[10] *= ratio;
   f2->SetParameters(f2params);
 //  h2->Fit("f2");
TFitResultPtr r = h2->Fit(f2, "S");

   f2->Draw("cont1 same");

// ----- end of original example code --------

TMatrixDSym cov = r->GetCovarianceMatrix(); 	
r->Print("V"); 	

TVirtualFitter *fitter = TVirtualFitter::GetFitter();

TMatrixD covMatrix1(5,5);
TMatrixD covMatrix2(5,5);
TMatrixD covMatrix3(5,5);

for (Int_t i=0; i<5; i++){
    for (Int_t j=0; j<5; j++){
      covMatrix1[i][j] = fitter->GetCovarianceMatrixElement(i,j);
      covMatrix2[i][j] = fitter->GetCovarianceMatrixElement(i+5,j+5);
      covMatrix3[i][j] = fitter->GetCovarianceMatrixElement(i+10,j+10);
    }
}				

covMatrix1.Print();
covMatrix2.Print();
covMatrix3.Print();

TF2 *fg1 = new TF2("fg1",g2,-10,10,-10,10, 5);		
TF2 *fg2 = new TF2("fg2",g2,-10,10,-10,10, 5);		
TF2 *fg3 = new TF2("fg3",g2,-10,10,-10,10, 5);		 

fg1->SetParameters( f2->GetParameter(0),  f2->GetParameter(1), f2->GetParameter(2), f2->GetParameter(3), f2->GetParameter(4)) ;
fg2->SetParameters( f2->GetParameter(5),  f2->GetParameter(6), f2->GetParameter(7), f2->GetParameter(8), f2->GetParameter(9)) ;
fg3->SetParameters( f2->GetParameter(10),  f2->GetParameter(11), f2->GetParameter(12), f2->GetParameter(13), f2->GetParameter(14)) ;

float binsize = h2->GetXaxis()->GetBinWidth(1)*h2->GetYaxis()->GetBinWidth(1);
Double_t a[2] = {-10,-10};
Double_t b[2] = {10,10};

cout<<"\n"<<"Integral of entire fit f2: "<<f2->Integral(-10,10,-10,10)/binsize<<"\n" 
<<"Error integral of entire fit f2: "<<f2->IntegralError(2,a,b,f2->GetParameters(), cov.GetMatrixArray())/binsize<<"\n";

cout<<"\n"<<"Integral of subfunction fg1: "<<fg1->Integral(-10,10,-10,10)/binsize<<"\n" 
<<"Error integral of subfunction fg1: "<<fg1->IntegralError(2,a,b,fg1->GetParameters(), covMatrix1.GetMatrixArray())/binsize<<"\n";

cout<<"\n"<<"Integral of subfunction fg2: "<<fg2->Integral(-10,10,-10,10)/binsize<<"\n" 
<<"Error integral of subfunction fg2: "<<fg2->IntegralError(2,a,b,fg2->GetParameters(), covMatrix2.GetMatrixArray())/binsize<<"\n";

cout<<"\n"<<"Integral of subfunction fg3: "<<fg3->Integral(-10,10,-10,10)/binsize<<"\n" 
<<"Error integral of subfunction fg3: "<<fg3->IntegralError(2,a,b,fg3->GetParameters(), covMatrix3.GetMatrixArray())/binsize<<"\n";
}

The output:

$ root fit2.C++
root [0] 
Processing fit2.C++...
Info in <TUnixSystem::ACLiC>: creating shared library ./fit2_C.so
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
 FCN=1070.75 FROM MIGRAD    STATUS=CONVERGED     383 CALLS         384 TOTAL
                     EDM=6.11771e-08    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           3.93353e+02   2.01408e+00  -5.10889e-04  -5.57626e-05
   2  p1          -2.99826e+00   1.11673e-02  -7.12611e-06  -1.53188e-02
   3  p2           2.98248e+00   7.74670e-03  -4.36019e-07   1.68535e-02
   4  p3          -3.00160e+00   1.12693e-02  -8.26075e-06  -6.57076e-03
   5  p4           2.96923e+00   7.59681e-03   9.07255e-09   9.64234e-04
   6  p5           6.19687e+02   1.00395e+01   1.33223e-03   1.49900e-05
   7  p6           5.76875e-03   1.10874e-02  -1.01210e-06   1.28961e-02
   8  p7           8.01426e-01   9.52950e-03   4.04935e-06  -2.55908e-02
   9  p8          -1.33022e-03   1.23516e-02  -1.14423e-06   5.01667e-04
  10  p9           8.99011e-01   1.03696e-02  -1.86742e-06   5.99275e-04
  11  p10          1.53371e+02   4.12827e+00   1.14622e-03  -5.99469e-05
  12  p11          3.99460e+00   1.78039e-02   1.45363e-05   4.62300e-03
  13  p12          7.10306e-01   1.23898e-02   6.16665e-07   1.46045e-02
  14  p13          4.01838e+00   1.68437e-02   2.62466e-06  -6.82577e-03
  15  p14          6.83628e-01   1.22030e-02  -6.60457e-06   1.20282e-02

****************************************
Minimizer is Minuit / Migrad
Chi2                      =      1070.75
NDf                       =         1130
Edm                       =  6.11771e-08
NCalls                    =          384
p0                        =      393.353   +/-   2.01408     
p1                        =     -2.99826   +/-   0.0111673   
p2                        =      2.98248   +/-   0.0077467   
p3                        =      -3.0016   +/-   0.0112693   
p4                        =      2.96923   +/-   0.00759681  
p5                        =      619.687   +/-   10.0395     
p6                        =   0.00576875   +/-   0.0110874   
p7                        =     0.801426   +/-   0.0095295   
p8                        =  -0.00133022   +/-   0.0123516   
p9                        =     0.899011   +/-   0.0103696   
p10                       =      153.371   +/-   4.12827     
p11                       =       3.9946   +/-   0.0178039   
p12                       =     0.710306   +/-   0.0123898   
p13                       =      4.01838   +/-   0.0168437   
p14                       =     0.683628   +/-   0.012203    

Covariance Matrix:

            	          p0          p1          p2          p3          p4          p5          p6          p7          p8          p9         p10         p11         p12         p13         p14
p0          	      4.0565   0.0031661  -0.0082394   0.0027842  -0.0081168      2.2167   0.0028641  -0.0036563   0.0039919  -0.0049322     -1.7534 -0.00067286   0.0036741  -0.0003818   0.0031724
p1          	   0.0031661  0.00012471 -4.9989e-06  1.9897e-05 -5.3584e-06   0.0086387  7.5216e-06 -2.9335e-05  9.9454e-06 -2.0469e-05  -0.0014559  8.3964e-06 -4.3859e-06  1.6861e-06  3.3498e-06
p2          	  -0.0082394 -4.9989e-06  6.0011e-05  3.4692e-07  8.4019e-06  -0.0052454 -6.8294e-06  7.6388e-07 -5.2095e-06  5.6642e-06   0.0051031 -3.7182e-06 -1.0726e-06  5.7314e-06 -2.1386e-05
p3          	   0.0027842  1.9897e-05  3.4692e-07    0.000127 -7.1965e-06   0.0074377  1.2279e-05 -2.8866e-05  6.0175e-06 -2.0476e-05  -0.0016617  4.4035e-06  1.0273e-06  4.7362e-06 -1.2329e-06
p4          	  -0.0081168 -5.3584e-06  8.4019e-06 -7.1965e-06  5.7712e-05  -0.0029658 -3.1962e-06  1.3181e-06  -1.188e-05   5.371e-06   0.0054026  1.0492e-05 -2.4204e-05  1.5327e-06 -1.9143e-07
p5          	      2.2167   0.0086387  -0.0052454   0.0074377  -0.0029658      100.79   0.0063621   -0.048523   0.0095393   -0.057189      14.943  0.00024441    -0.02817 -0.00042756   -0.024943
p6          	   0.0028641  7.5216e-06 -6.8294e-06  1.2279e-05 -3.1962e-06   0.0063621  0.00012293 -1.7642e-05  8.0568e-06 -1.4193e-05 -8.3104e-05  2.6821e-06 -5.9384e-06  5.1839e-06  5.7325e-06
p7          	  -0.0036563 -2.9335e-05  7.6388e-07 -2.8866e-05  1.3181e-06   -0.048523 -1.7642e-05  9.0811e-05 -1.8324e-05  2.7857e-05   -0.007705 -2.8057e-07  1.9429e-05 -1.0566e-06  1.0269e-05
p8          	   0.0039919  9.9454e-06 -5.2095e-06  6.0175e-06  -1.188e-05   0.0095393  8.0568e-06 -1.8324e-05  0.00015256 -2.1903e-05  0.00049352  2.4928e-07  5.6046e-06  5.7185e-06 -7.2917e-06
p9          	  -0.0049322 -2.0469e-05  5.6642e-06 -2.0476e-05   5.371e-06   -0.057189 -1.4193e-05  2.7857e-05 -2.1903e-05  0.00010753   -0.011057 -5.1289e-06  1.8529e-05 -4.0345e-06  2.2539e-05
p10         	     -1.7534  -0.0014559   0.0051031  -0.0016617   0.0054026      14.943 -8.3104e-05   -0.007705  0.00049352   -0.011057      17.043   0.0031401   -0.018033  0.00029556   -0.019613
p11         	 -0.00067286  8.3964e-06 -3.7182e-06  4.4035e-06  1.0492e-05  0.00024441  2.6821e-06 -2.8057e-07  2.4928e-07 -5.1289e-06   0.0031401  0.00031698 -3.1303e-05  2.4583e-05 -1.1922e-06
p12         	   0.0036741 -4.3859e-06 -1.0726e-06  1.0273e-06 -2.4204e-05    -0.02817 -5.9384e-06  1.9429e-05  5.6046e-06  1.8529e-05   -0.018033 -3.1303e-05  0.00015351 -1.7482e-06 -2.2495e-05
p13         	  -0.0003818  1.6861e-06  5.7314e-06  4.7362e-06  1.5327e-06 -0.00042756  5.1839e-06 -1.0566e-06  5.7185e-06 -4.0345e-06  0.00029556  2.4583e-05 -1.7482e-06  0.00028371 -1.5889e-05
p14         	   0.0031724  3.3498e-06 -2.1386e-05 -1.2329e-06 -1.9143e-07   -0.024943  5.7325e-06  1.0269e-05 -7.2917e-06  2.2539e-05   -0.019613 -1.1922e-06 -2.2495e-05 -1.5889e-05  0.00014891

Correlation Matrix:

            	          p0          p1          p2          p3          p4          p5          p6          p7          p8          p9         p10         p11         p12         p13         p14
p0          	           1     0.14077    -0.52808     0.12267    -0.53049     0.10963     0.12826     -0.1905     0.16047    -0.23616    -0.21088   -0.018764     0.14723   -0.011254     0.12908
p1          	     0.14077           1   -0.057785      0.1581   -0.063162    0.077053    0.060748    -0.27566    0.072103    -0.17676   -0.031581    0.042231   -0.031699   0.0089641    0.024582
p2          	    -0.52808   -0.057785           1   0.0039739     0.14277   -0.067445   -0.079513    0.010348   -0.054445    0.070511     0.15957   -0.026959   -0.011176    0.043925    -0.22623
p3          	     0.12267      0.1581   0.0039739           1   -0.084061     0.06574     0.09827    -0.26879    0.043231    -0.17522   -0.035719    0.021948   0.0073573    0.024952  -0.0089654
p4          	    -0.53049   -0.063162     0.14277   -0.084061           1   -0.038886   -0.037947    0.018207    -0.12661    0.068181     0.17227    0.077575    -0.25715    0.011978   -0.002065
p5          	     0.10963    0.077053   -0.067445     0.06574   -0.038886           1    0.057156    -0.50719    0.076928    -0.54933     0.36053   0.0013674    -0.22647  -0.0025284     -0.2036
p6          	     0.12826    0.060748   -0.079513     0.09827   -0.037947    0.057156           1    -0.16697    0.058831    -0.12344  -0.0018156    0.013587   -0.043229    0.027758    0.042369
p7          	     -0.1905    -0.27566    0.010348    -0.26879    0.018207    -0.50719    -0.16697           1    -0.15568      0.2819    -0.19585  -0.0016537     0.16456  -0.0065826    0.088308
p8          	     0.16047    0.072103   -0.054445    0.043231    -0.12661    0.076928    0.058831    -0.15568           1    -0.17101   0.0096786   0.0011336    0.036623    0.027487   -0.048377
p9          	    -0.23616    -0.17676    0.070511    -0.17522    0.068181    -0.54933    -0.12344      0.2819    -0.17101           1    -0.25829   -0.027781     0.14422   -0.023099     0.17811
p10         	    -0.21088   -0.031581     0.15957   -0.035719     0.17227     0.36053  -0.0018156    -0.19585   0.0096786    -0.25829           1    0.042723    -0.35256   0.0042504    -0.38932
p11         	   -0.018764    0.042231   -0.026959    0.021948    0.077575   0.0013674    0.013587  -0.0016537   0.0011336   -0.027781    0.042723           1    -0.14191    0.081974  -0.0054876
p12         	     0.14723   -0.031699   -0.011176   0.0073573    -0.25715    -0.22647   -0.043229     0.16456    0.036623     0.14422    -0.35256    -0.14191           1  -0.0083768    -0.14878
p13         	   -0.011254   0.0089641    0.043925    0.024952    0.011978  -0.0025284    0.027758  -0.0065826    0.027487   -0.023099   0.0042504    0.081974  -0.0083768           1   -0.077304
p14         	     0.12908    0.024582    -0.22623  -0.0089654   -0.002065     -0.2036    0.042369    0.088308   -0.048377     0.17811    -0.38932  -0.0054876    -0.14878   -0.077304           1

5x5 matrix is as follows

     |      0    |      1    |      2    |      3    |      4    |
----------------------------------------------------------------------
   0 |      4.057    0.003166   -0.008239    0.002784   -0.008117 
   1 |   0.003166   0.0001247  -4.999e-06    1.99e-05  -5.358e-06 
   2 |  -0.008239  -4.999e-06   6.001e-05   3.469e-07   8.402e-06 
   3 |   0.002784    1.99e-05   3.469e-07    0.000127  -7.197e-06 
   4 |  -0.008117  -5.358e-06   8.402e-06  -7.197e-06   5.771e-05 


5x5 matrix is as follows

     |      0    |      1    |      2    |      3    |      4    |
----------------------------------------------------------------------
   0 |      100.8    0.006362    -0.04852    0.009539    -0.05719 
   1 |   0.006362   0.0001229  -1.764e-05   8.057e-06  -1.419e-05 
   2 |   -0.04852  -1.764e-05   9.081e-05  -1.832e-05   2.786e-05 
   3 |   0.009539   8.057e-06  -1.832e-05   0.0001526   -2.19e-05 
   4 |   -0.05719  -1.419e-05   2.786e-05   -2.19e-05   0.0001075 


5x5 matrix is as follows

     |      0    |      1    |      2    |      3    |      4    |
----------------------------------------------------------------------
   0 |      17.04     0.00314    -0.01803   0.0002956    -0.01961 
   1 |    0.00314    0.000317   -3.13e-05   2.458e-05  -1.192e-06 
   2 |   -0.01803   -3.13e-05   0.0001535  -1.748e-06  -2.249e-05 
   3 |  0.0002956   2.458e-05  -1.748e-06   0.0002837  -1.589e-05 
   4 |   -0.01961  -1.192e-06  -2.249e-05  -1.589e-05   0.0001489 


Integral of entire fit f2: 99013.4
Error integral of entire fit f2: 309.732

Integral of subfunction fg1: 85920.4
Error integral of subfunction fg1: 316.345

Integral of subfunction fg2: 11221.2
Error integral of subfunction fg2: 164.356

Integral of subfunction fg3: 1871.75
Error integral of subfunction fg3: 43.883
root [1]

moneta · May 22, 2014, 8:16pm

Hi,

Here are the answer:

The result of the macro seems correct to me
The other elements are related to parameters on which the function do not depend, so in principle it should be fine. However using the covariance matrix and the error propagation in IntegralError are only approximation. The correct thing to do is to fit directly those quantity you are interested and compute the uncertainty by profiling the other parameters (i.e. running Minos)
Yes, TVirtualFitter has a convention that removes the fixed parameter from the matrix. I would not use it since it will complicate the code.
Just get the covariance matrix from TFitResult. You can use FItResult::CovMatrix(i,j) to retrieve the single element.

Best Regards

Lorenzo

p_x · May 23, 2014, 7:57am

Dear Lorenzo,

thank you very much for the clear answers and for your time.

Best regards