STAT/AliTMinuitToolkit.cxx

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