6feb61882aafd143f3c6bd63b901427a518c6022
[u/mrichter/AliRoot.git] / 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(0),
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(0),
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 == 0) {
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