missing commits
[u/mrichter/AliRoot.git] / 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