STAT/AliTMinuitToolkit.cxx

   1
   2 #include "AliTMinuitToolkit.h"
   3 #include <TNamed.h>
   4 #include <TVirtualFitter.h>
   5 #include <TH1F.h>
   6 #include <TH2F.h>
   7 #include <TF1.h>
   8 #include <TFormula.h>
   9 #include <TVector.h>
  10 #include <TMatrix.h>
  11 #include <TMath.h>
  12 #include <TString.h>
  13 #include <TROOT.h>
  14 #include <TCanvas.h>
  15 #include <TRandom.h>
  16
  17
  18 //--------------------------------------------------------------------------------------
  19 //
  20 // The AliTMinuitToolkit serves as an easy to use interface for the TMinuit
  21 // package:
  22 //
  23 // - It allows to fit a curve to one and two dimensional histograms
  24 //   (TH2F::Fit() only allows to fit a hyperplane).
  25 // - Or n points can be specified directly via a n x 2 matrix.
  26 // - An option for robust fitting with non-linear functions is implemented.
  27 //
  28 // A small example illustrating the usage of AliTMinuitToolkit is given in the function
  29 // "AliTMinuitToolkit::Test()".
  30 //
  31 //
  32 // 1. Setting the formula:
  33 //
  34 //  The formula is simply set via "void SetFitFunction(TFormula * formula)".
  35 //
  36 //
  37 // 2. Adding the data points
  38 //
  39 //  - In order to fit a histogram, use "void FitHistogram(TH1F * his)" or
  40 //    "void FitHistogram(TH2F * his)". The fitter is started automatically
  41 //  - Alternatively, the direct specification of the points is possible via
  42 //    "void SetPoints(TMatrixD * points)". Note, that the each point
  43 //    corresponds to one row in the matrix. The fitter is then started with
  44 //    the command "void Fit()". The weight of each point can be specified
  45 //    with an n-dimensional vector using "void SetWeights(TVectorD * weights)".
  46 //
  47 //
  48 // 3. Accessing the fit results
  49 //
  50 //  The N parameters of the formula are stored in a N-dimensional vector which
  51 //  is returned by "TVectorD * GetParameters()". In a similar way the covariance
  52 //  matrix of the fit is returned via "TMatrixD * GetCovarianceMatrix()" which
  53 //  is of the type N x N.
  54 //
  55 //
  56 // 4. Non-linear robust fitting:
  57 //
  58 //  Even a few outliers can lead to wrong results of a least-squares fitting
  59 //  procedure. In this case the use of robust(resistant) methods can be
  60 //  helpful, but a stronger dependence on starting values or convergence to
  61 //  local minima can occur.
  62 //
  63 //  The robust option becomes active if EnableRobust(true, sigma) is called. It is
  64 //  very much recommended that a normalization value (scale variable) corresponding
  65 //  to an expected deviation (sigma) is specified via
  66 //  "EnableRobust(Bool_t b, Double_t sigma)".
  67 //
  68 //  Performing the fit without knowledge of sigma is also possible if only
  69 //  "EnableRobust(true)" is activated, but this is NOT RECOMMENDED.
  70 //
  71 //  The method is based on another estimator instead of chi^2. For small deviations
  72 //  the function behaves like x^2 and for larger deviations like |x| - the so
  73 //  called Huber estimator:
  74 //
  75 //   h(x) = x^2                              , for x < 2.5*sigma
  76 //   h(x) = 2*(2.5*sigma)*x - (2.5*sigma)^2  , for x > 2.5*sigma
  77 //
  78 //  If a weighting function is specified in addition, a second run with the ordinary
  79 //  metric is started, but before entering the iteration every point is weighted
  80 //  according to its distance to the outcoming function of the first run. The weighting
  81 //  function w(x) must be defined on the intervall x in [0,1]. w(0) then
  82 //  corresponds to the weight of the closest point and w(1) to the point with the
  83 //  largest distance.
  84 //
  85 //  Some standard weighting functions are predefined in
  86 //  "SetWeightFunction(Char_t * name, Float_t param1, Float_t param2 = 0)":
  87 //   - "BOX" equals to 1 if x < param1 and to 0 if x > param1.
  88 //   - "EXPONENTIAL" corresponds to "Math::Exp(-TMath::Log(param1)*x)"
  89 //   - "ERRORFUNCTION" corresponds to "TMath::Erfc((x-param1)/param2)"
  90 //
  91 //
  92 //  REFERENCE for non-linear robust fitting:
  93 //  Ekblom H. and Madsen K. (1988), Alogrithms for non-linear Huber estimation,
  94 //  BIT Numerical Mathematics 29 (1989) 60-76.
  95 //  internet: http://www.springerlink.com/content/m277218542988344/
  96 //
  97 //
  98 // 5. examples:
  99 //
 100 //  A small example illustrating the working principles of AliTMinuitToolkit is given
 101 //  in the function "AliTMinuitToolkit::Test()".
 102 //
 103 //
 104 //
 105 // Comments and questions are always welcome: A.Kalweit@gsi.de
 106 //--------------------------------------------------------------------------------------
 107
 108
 109 ClassImp(AliTMinuitToolkit)
 110
 111 AliTMinuitToolkit::AliTMinuitToolkit() :
 112    TNamed(),
 113    fFormula(0),
 114    fWeightFunction(0),
 115    fFitAlgorithm(0),
 116    fPoints(0),
 117    fWeights(0),
 118    fParam(0),
 119    fParamLimits(0),
 120    fCovar(0),
 121    fChi2(0),
 122    fMaxCalls(0),
 123    fPrecision(0),
 124    fUseRobust(0),
 125    fExpectedSigma(0)
 126 {
 127  //
 128  // standard constructor
 129  //
 130  fMaxCalls = 500;
 131  fPrecision = 1;
 132  fUseRobust = false;
 133  fExpectedSigma = 0;
 134 }
 135
 136
 137 AliTMinuitToolkit::AliTMinuitToolkit(const AliTMinuitToolkit&) :
 138    TNamed(),
 139    fFormula(0),
 140    fWeightFunction(0),
 141    fFitAlgorithm(0),
 142    fPoints(0),
 143    fWeights(0),
 144    fParam(0),
 145    fParamLimits(0),
 146    fCovar(0),
 147    fChi2(0),
 148    fMaxCalls(0),
 149    fPrecision(0),
 150    fUseRobust(0),
 151    fExpectedSigma(0)
 152 {
 153
 154
 155 }
 156
 157
 158 AliTMinuitToolkit& AliTMinuitToolkit::operator=(const AliTMinuitToolkit&) {
 159
 160  return *this;
 161 }
 162
 163
 164
 165 AliTMinuitToolkit::~AliTMinuitToolkit(){
 166   //
 167   // destructor
 168   //
 169   delete fPoints;
 170   delete fWeights;
 171   delete fWeightFunction;
 172   delete fParamLimits;
 173   delete fFormula;
 174   delete fParam;
 175   delete fCovar;
 176   delete fChi2;
 177 }
 178
 179 void AliTMinuitToolkit::FitHistogram(TH1F * his) {
 180  //
 181  // Fit a one dimensional histogram
 182  //
 183  fPoints = new TMatrixD(his->GetNbinsX(), 2);
 184
 185  for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
 186   Double_t x = his->GetXaxis()->GetBinCenter(ibin+1);
 187   Double_t y = his->GetBinContent(ibin+1);
 188
 189   (*fPoints)(ibin, 0) = x;
 190   (*fPoints)(ibin, 1) = y;
 191  }
 192
 193  Fit();
 194 }
 195
 196
 197 void AliTMinuitToolkit::FitHistogram(TH2F * his) {
 198  //
 199  // Fit a curve to a two dimensional histogram
 200  //
 201  fPoints = new TMatrixD((Long64_t)his->GetEntries(), 2);
 202  Long64_t entry = 0;
 203
 204  for(Int_t ibin=0; ibin < his->GetNbinsX(); ibin++) {
 205   Double_t x = his->GetXaxis()->GetBinCenter(ibin);
 206   for(Int_t jbin=0; jbin < his->GetNbinsY(); jbin++) {
 207    Long64_t n = his->GetBin(ibin, jbin);
 208    Double_t y = his->GetYaxis()->GetBinCenter(jbin);
 209    for(Int_t ientries=0; ientries < his->GetBinContent(n); ientries++) {
 210     (*fPoints)(entry,0) = x;
 211     (*fPoints)(entry,1) = y;
 212     entry++;
 213    }
 214
 215   }
 216  }
 217
 218  Fit();
 219 }
 220
 221
 222 void AliTMinuitToolkit::SetWeightFunction(Char_t * name, Float_t param1, Float_t param2) {
 223  //
 224  // Set the weight function which must be defined on the interval [0,1].
 225  //
 226  TString FuncType(name);
 227  FuncType.ToUpper();
 228
 229  if (FuncType == "EXPONENTIAL") fWeightFunction = new TFormula("exp", Form("TMath::Exp(-TMath::Log(%f)*x)", param1));
 230  if (FuncType == "BOX") fWeightFunction = new TFormula("box", Form("TMath::Erfc((x-%f)/0.0001)", param1));
 231  if (FuncType == "ERRORFUNCTION") fWeightFunction = new TFormula("err", Form("TMath::Erfc((x-%f)/%f)", param1, param2));
 232
 233 }
 234
 235
 236 void AliTMinuitToolkit::FitterFCN(int &npar, double *dummy, double &fchisq, double *gin, int iflag){
 237   //
 238   // internal function which gives the specified function to the TMinuit function
 239   //
 240
 241   // suppress warnings for unused variables:
 242   dummy = dummy;
 243   iflag = iflag;
 244   npar = npar;
 245   //
 246   AliTMinuitToolkit * fitter = (AliTMinuitToolkit*)TVirtualFitter::GetFitter()->GetObjectFit();
 247   fchisq = 0;
 248   Int_t nvar       = fitter->GetPoints()->GetNcols()-1;
 249   Int_t npoints    = fitter->GetPoints()->GetNrows();
 250
 251   // calculate mean deviation for normalization or use user-defined sigma
 252   Double_t dev = 0.;
 253   if (fitter->GetExpectedSigma() == 0 && fitter->GetStatus() == true) {
 254    for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 255     Double_t x[100];
 256      for (Int_t ivar=0; ivar<nvar; ivar++){
 257       x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
 258      }
 259     Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
 260     Double_t delta = (*fitter->GetPoints())(ipoint, nvar) - funx;
 261     dev += TMath::Sqrt(TMath::Abs(delta));
 262    }
 263    dev = dev/npoints;
 264   } else {
 265    dev = fitter->GetExpectedSigma();
 266   }
 267   // calculate chisquare
 268   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 269     Double_t x[100];
 270     for (Int_t ivar=0; ivar<nvar; ivar++){
 271       x[ivar] = (*fitter->GetPoints())(ipoint, ivar);
 272     }
 273     Float_t funx = fitter->GetFormula()->EvalPar(x,gin);
 274     Double_t delta = TMath::Abs((*fitter->GetPoints())(ipoint, nvar) - funx);
 275     if (fitter->GetStatus() == true) {
 276      delta = delta/dev; // normalization
 277      if (delta <= 2.5) fchisq+= delta*delta; // new metric: Huber-k-estimator
 278      if (delta > 2.5) fchisq+= 2*(2.5)*delta - (2.5*2.5);
 279     } else {
 280      Double_t weight = (*fitter->GetWeights())(ipoint);
 281      fchisq+= delta*delta*weight; //old metric
 282     }
 283   }
 284  }
 285
 286
 287 void AliTMinuitToolkit::Fit() {
 288   //
 289   // internal function that calls the fitter
 290   //
 291   Int_t nparam  = fParam->GetNrows();
 292   Int_t npoints = fPoints->GetNrows();
 293   Int_t nvar    = fPoints->GetNcols()-1;
 294
 295   // set all paramter limits to infinity as default
 296   if (fParamLimits == 0) {
 297    fParamLimits = new TMatrixD(nparam ,2);
 298    for (Int_t iparam=0; iparam<nparam; iparam++){
 299     (*fParamLimits)(iparam, 0) = 0;
 300     (*fParamLimits)(iparam, 1) = 0;
 301    }
 302   }
 303
 304   // set all weights to 1 as default
 305   Bool_t weightFlag = false;
 306   if (fWeightFunction == 0) {
 307    fWeightFunction = new TFormula("constant", "1");
 308   } else {
 309    weightFlag = true;
 310   }
 311
 312   // migrad fit algorithm as default
 313   if (fFitAlgorithm == 0) {
 314    fFitAlgorithm = "migrad";
 315   }
 316
 317   // assign weights
 318   if (fWeights == 0) {
 319    fWeights = new TVectorD(npoints);
 320    for (Int_t ipoint=0; ipoint<npoints; ipoint++) (*fWeights)(ipoint) = 1;
 321   }
 322
 323   // set up the fitter
 324   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, nparam);
 325   minuit->SetObjectFit(this);
 326   minuit->SetFCN(AliTMinuitToolkit::FitterFCN);
 327
 328   // initialize paramters (step size???)
 329   for (Int_t iparam=0; iparam<nparam; iparam++){
 330    minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
 331   }
 332
 333   //
 334   Double_t argList[2];
 335   argList[0] = fMaxCalls; //maximal number of calls
 336   argList[1] = fPrecision; //tolerance normalized to 0.001
 337   if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
 338   if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
 339   // two additional arguments can be specified ExecuteCommand("migrad", argList, 2) - use 0,0 for default
 340
 341   // fill parameter vector
 342   for (Int_t ivar=0; ivar<nparam; ivar++){
 343    (*fParam)(ivar) = minuit->GetParameter(ivar);
 344    fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
 345   }
 346
 347   // if a weight function is specified -> enter 2nd run with weights
 348   if (weightFlag == true && fUseRobust == true) {
 349    // sort points for weighting
 350    Double_t *sortList = new Double_t[npoints];
 351    Int_t  *indexList = new Int_t[npoints];
 352    for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 353     Double_t funx = fFormula->Eval((*fPoints)(ipoint, 0));
 354     Double_t delta = TMath::Abs((*fPoints)[ipoint][nvar] - funx);
 355     sortList[ipoint] = delta;
 356    }
 357    TMath::Sort(npoints, sortList, indexList, false);
 358    for (Int_t ip=0; ip<npoints; ip++){
 359     Double_t t = ip/(Double_t)npoints;
 360     (*fWeights)(indexList[ip]) = fWeightFunction->Eval(t);
 361    }
 362
 363    // set up the fitter
 364    fUseRobust = false;
 365    for (Int_t iparam=0; iparam<nparam; iparam++){
 366     minuit->SetParameter(iparam, Form("p[%d]",iparam), (*fParam)(iparam), (*fParam)(iparam)/10, (*fParamLimits)(iparam, 0), (*fParamLimits)(iparam, 1));
 367    }
 368    // start fitting
 369    if (fMaxCalls == 500 && fPrecision == 1) minuit->ExecuteCommand(fFitAlgorithm, 0, 0);
 370    if (fMaxCalls != 500 || fPrecision != 1) minuit->ExecuteCommand(fFitAlgorithm, argList, 2);
 371    fUseRobust = true;
 372
 373    delete sortList;
 374    delete indexList;
 375   }
 376
 377   // fill parameter vector
 378   for (Int_t ivar=0; ivar<nparam; ivar++){
 379    (*fParam)(ivar) = minuit->GetParameter(ivar);
 380    fFormula->SetParameter(ivar, minuit->GetParameter(ivar));
 381   }
 382
 383   // fill covariance matrix
 384   fCovar = new TMatrixD(nparam, nparam);
 385   for(Int_t i=0; i < nparam; i++) {
 386    for(Int_t j=0; j < nparam; j++) {
 387     (*fCovar)(i,j) = minuit->GetCovarianceMatrixElement(i,j);
 388    }
 389   }
 390
 391   if (weightFlag == false) fWeightFunction = 0;
 392 }
 393
 394
 395
 396 void AliTMinuitToolkit::Test() {
 397  //
 398  // This test function shows the basic working principles of this class
 399  // and illustrates how a robust fit can improve the results
 400  //
 401
 402  // 1. provide some example histogram
 403  TH1F * hist = new TH1F("test", "with (red) and without (black) robust option", 20,0,4);
 404  TRandom * rand = new TRandom();
 405  for (Int_t i = 0; i < 10000; i++) {
 406   hist->Fill(rand->Exp(1));
 407   if (i < 1000) hist->Fill(3); //"outliers"
 408   if (i < 1070) hist->Fill(3.5);
 409   if (i < 670) hist->Fill(2);
 410   if (i < 770) hist->Fill(1.5);//"outliers"
 411   if (i < 740) hist->Fill(1);
 412  }
 413  TCanvas * canv = new TCanvas();
 414  canv->cd(1);
 415  hist->Draw();
 416
 417  // 2. example fit without robust option
 418  AliTMinuitToolkit * tool = new AliTMinuitToolkit();
 419  TFormula *FormExp = new TFormula("formExp", "[0]*TMath::Exp(-[1]*x)");
 420  tool->SetFitFunction(FormExp);
 421  TVectorD *vec1 = new TVectorD(2); // Set initial values
 422  (*vec1)(0) = 1800;
 423  (*vec1)(1) = 1;
 424  tool->SetInitialParam(vec1);
 425  tool->FitHistogram(hist);
 426
 427  // draw fit function
 428  TF1 *func = new TF1("test", "[0]*TMath::Exp(-[1]*x)", 0, 6);
 429  func->SetParameters((*tool->GetParameters())(0), (*tool->GetParameters())(1));
 430  func->Draw("same");
 431
 432  // 3 . robust fit
 433  TVectorD *vec2 = new TVectorD(2);
 434  (*vec2)(0) = 1800;
 435  (*vec2)(1) = 1;
 436  tool->SetInitialParam(vec2);
 437  tool->EnableRobust(true, 10);
 438  tool->SetWeightFunction("box", 0.75);
 439  tool->FitHistogram(hist);
 440  TF1 *func2 = new TF1("test2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
 441  func2->SetParameter(0, (*tool->GetParameters())(0));
 442  func2->SetParameter(1, (*tool->GetParameters())(1));
 443  func2->SetLineColor(kRed);
 444  func2->Draw("same");
 445
 446 }
 447