[root] / trunk / math / mathcore / src / FitResult.cxx Repository:
ViewVC logotype

Diff of /trunk/math/mathcore/src/FitResult.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 25485, Mon Sep 22 07:52:52 2008 UTC revision 25486, Mon Sep 22 12:43:03 2008 UTC
# Line 14  Line 14 
14    
15  #include "Fit/FitConfig.h"  #include "Fit/FitConfig.h"
16    
17    #include "Fit/BinData.h"
18    
19  #include "Math/Minimizer.h"  #include "Math/Minimizer.h"
20    
21  #include "Math/IParamFunction.h"  #include "Math/IParamFunction.h"
22    #include "Math/OneDimFunctionAdapter.h"
23    
24  #include "Math/DistFunc.h"  #include "Math/DistFunc.h"
25    
26    #include "TMath.h"
27    #include "Math/RichardsonDerivator.h"
28    #include "Math/Error.h"
29    
30  #include <cassert>  #include <cassert>
31  #include <cmath>  #include <cmath>
32    #include <iostream>
33    #include <iomanip>
34    
35  namespace ROOT {  namespace ROOT {
36    
# Line 29  Line 38 
38    
39    
40  FitResult::FitResult() :  FitResult::FitResult() :
41     fValid(false), fNormalized(false), fVal(0), fEdm(0), fChi2(0), fNdf(0), fNCalls(0), fDataSize(0), fFitFunc(0)     fValid(false), fNormalized(false), fNFree(0), fNdf(0), fNCalls(0),
42       fStatus(-1), fVal(0), fEdm(0), fChi2(0), fFitFunc(0)
43  {  {
44     // Default constructor implementation.     // Default constructor implementation.
45  }  }
46    
47        FitResult::FitResult(ROOT::Math::Minimizer & min, const FitConfig & fconfig, const IModelFunction & func,  bool isValid,  unsigned int sizeOfData, const  ROOT::Math::IMultiGenFunction * chi2func, bool minosErr, unsigned int ncalls ) :        FitResult::FitResult(ROOT::Math::Minimizer & min, const FitConfig & fconfig, IModelFunction & func,  bool isValid,  unsigned int sizeOfData, const  ROOT::Math::IMultiGenFunction * chi2func, bool minosErr, unsigned int ncalls ) :
48     fValid(isValid),     fValid(isValid),
49     fNormalized(false),     fNormalized(false),
50       fNFree(min.NFree() ),
51       fNCalls(min.NCalls()),
52       fStatus(min.Status() ),
53     fVal (min.MinValue()),     fVal (min.MinValue()),
54     fEdm (min.Edm()),     fEdm (min.Edm()),
55     fNCalls(min.NCalls()),     fFitFunc(&func),
56     fParams(std::vector<double>(min.X(), min.X() + min.NDim() ) ),     fParams(std::vector<double>(min.X(), min.X() + min.NDim() ) )
    fFitFunc(&func)  
57  {  {
58     // Constructor from a minimizer, fill the data     // Constructor from a minimizer, fill the data. ModelFunction  is passed as non const
59       // since it will be managed by the FitResult
60    
61     if (sizeOfData > 0) fNdf = sizeOfData - min.NFree();     if (sizeOfData > 0) fNdf = sizeOfData - min.NFree();
62    
63       // set right parameters in function (in case minimizer did not do before)
64       // do also when fit is not valid
65       fFitFunc->SetParameters(&fParams.front());
66    
67     if (min.Errors() != 0)     if (min.Errors() != 0)
68        fErrors = std::vector<double>(min.Errors(), min.Errors() + min.NDim() ) ;        fErrors = std::vector<double>(min.Errors(), min.Errors() + min.NDim() ) ;
69    
70       // check for fixed or limited parameters
71       for (unsigned int ipar = 0; ipar < fParams.size(); ++ipar) {
72          const ParameterSettings & par = fconfig.ParSettings(ipar);
73          if (par.IsFixed() ) fFixedParams.push_back(ipar);
74          if (par.IsBound() ) fBoundParams.push_back(ipar);
75       }
76    
77     if (chi2func == 0)     if (chi2func == 0)
78        fChi2 = fVal;        fChi2 = fVal;
79     else {     else {
# Line 65  Line 89 
89  //    // fill error matrix  //    // fill error matrix
90     // cov matrix rank     // cov matrix rank
91     if (fValid) {     if (fValid) {
92    
93        unsigned int r = n * (  n + 1 )/2;        unsigned int r = n * (  n + 1 )/2;
94        fCovMatrix.reserve(r);        fCovMatrix.reserve(r);
95        for (unsigned int i = 0; i < n; ++i)        for (unsigned int i = 0; i < n; ++i)
# Line 72  Line 97 
97              fCovMatrix.push_back(min.CovMatrix(i,j) );              fCovMatrix.push_back(min.CovMatrix(i,j) );
98    
99        assert (fCovMatrix.size() == r );        assert (fCovMatrix.size() == r );
100     }  
101          // normalize errors if requested in configuration
102          if (fconfig.NormalizeErrors() ) NormalizeErrors();
103    
104     // minos errors     // minos errors
105     if (minosErr) {     if (minosErr) {
# Line 85  Line 112 
112        }        }
113     }     }
114    
115          // globalCC
116          fGlobalCC.reserve(n);
117          for (unsigned int i = 0; i < n; ++i) {
118             double globcc = min.GlobalCC(i);
119             if (globcc < 0) break; // it is not supported by that minimizer
120             fGlobalCC.push_back(globcc);
121          }
122    
123       }
124    
125     fMinimType = fconfig.MinimizerType();     fMinimType = fconfig.MinimizerType();
126    
127       // append algorithm name for minimizer that support it
128       if ( (fMinimType.find("Fumili") == std::string::npos) &&
129            (fMinimType.find("GSLMultiFit") == std::string::npos)
130          ) {
131     if (fconfig.MinimizerAlgoType() != "") fMinimType += " / " + fconfig.MinimizerAlgoType();     if (fconfig.MinimizerAlgoType() != "") fMinimType += " / " + fconfig.MinimizerAlgoType();
132  }  }
133    
134    }
135    
136    FitResult::FitResult(const FitResult &rhs) {
137       // Implementation of copy constructor
138       (*this) = rhs;
139    }
140    
141    FitResult & FitResult::operator = (const FitResult &rhs) {
142       // Implementation of assignment operator.
143       if (this == &rhs) return *this;  // time saving self-test
144    
145       // Manages the fitted function
146       if (fFitFunc) delete fFitFunc;
147       fFitFunc = 0;
148       if (rhs.fFitFunc != 0 ) {
149          fFitFunc = dynamic_cast<IModelFunction *>( (rhs.fFitFunc)->Clone() );
150          assert(fFitFunc != 0);
151       }
152    
153       // copy all other data members
154       fValid = rhs.fValid;
155       fNormalized = rhs.fNormalized;
156       fNFree = rhs.fNFree;
157       fNdf = rhs.fNdf;
158       fNCalls = rhs.fNCalls;
159       fStatus = rhs.fStatus;
160       fVal = rhs.fVal;
161       fEdm = rhs.fEdm;
162       fChi2 = rhs.fChi2;
163    
164       fFixedParams = rhs.fFixedParams;
165       fBoundParams = rhs.fBoundParams;
166       fParams = rhs.fParams;
167       fErrors = rhs.fErrors;
168       fCovMatrix = rhs.fCovMatrix;
169       fGlobalCC = rhs.fGlobalCC;
170       fMinosErrors = rhs.fMinosErrors;
171    
172       fMinimType = rhs.fMinimType;
173    
174       return *this;
175    
176    }
177    
178    
179  void FitResult::NormalizeErrors() {  void FitResult::NormalizeErrors() {
180     // normalize errors and covariance matrix according to chi2 value     // normalize errors and covariance matrix according to chi2 value
181     if (fNdf == 0 || fChi2 <= 0) return;     if (fNdf == 0 || fChi2 <= 0) return;
# Line 116  Line 203 
203     return -1; // case name is not found     return -1; // case name is not found
204  }  }
205    
206    bool FitResult::IsParameterBound(unsigned int ipar) const {
207       for (unsigned int i = 0; i < fBoundParams.size() ; ++i)
208          if ( fBoundParams[i] == ipar) return true;
209       return false;
210    }
211    
212    bool FitResult::IsParameterFixed(unsigned int ipar) const {
213       for (unsigned int i = 0; i < fFixedParams.size() ; ++i)
214          if ( fFixedParams[i] == ipar) return true;
215       return false;
216    }
217    
218  void FitResult::Print(std::ostream & os, bool doCovMatrix) const {  void FitResult::Print(std::ostream & os, bool doCovMatrix) const {
219     // print the result in the given stream     // print the result in the given stream
220       // need to add minos errors , globalCC, etc..
221     if (!fValid) {     if (!fValid) {
222        os << "\n****************************************\n";        os << "\n****************************************\n";
223        os << "            Invalid FitResult            ";        os << "            Invalid FitResult            ";
# Line 126  Line 226 
226     }     }
227    
228     os << "\n****************************************\n";     os << "\n****************************************\n";
229     os << "            FitResult                   \n\n";     //os << "            FitResult                   \n\n";
230     os << "Minimizer is " << fMinimType << std::endl;     os << "Minimizer is " << fMinimType << std::endl;
231     unsigned int npar = fParams.size();     unsigned int npar = fParams.size();
232     os << "Chi2/Likelihood  =\t" << fVal << std::endl;     const unsigned int nw = 25;
233     if (fVal != fChi2)     if (fVal != fChi2)
234     os << "Chi2             =\t" << fChi2<< std::endl;        os << std::setw(nw) << std::left << "Likelihood" << " =\t" << fVal << std::endl;
235     os << "NDf              =\t" << fNdf << std::endl;     os << std::setw(nw) << std::left <<  "Chi2"  << " =\t" << fChi2<< std::endl;
236     os << "Edm              =\t" << fEdm << std::endl;     os << std::setw(nw) << std::left << "NDf"    << " =\t" << fNdf << std::endl;
237     os << "NCalls           =\t" << fNCalls << std::endl;     if (fMinimType.find("Linear") == std::string::npos) {  // no need to print this for linear fits
238          os << std::setw(nw) << std::left << "Edm"    << " =\t" << fEdm << std::endl;
239          os << std::setw(nw) << std::left << "NCalls" << " =\t" << fNCalls << std::endl;
240       }
241     assert(fFitFunc != 0);     assert(fFitFunc != 0);
242     for (unsigned int i = 0; i < npar; ++i) {     for (unsigned int i = 0; i < npar; ++i) {
243        os << fFitFunc->ParameterName(i) << "\t\t =\t" << fParams[i] << " \t+/-\t" << fErrors[i] << std::endl;        os << std::setw(nw) << std::left << fFitFunc->ParameterName(i) << " =\t" << fParams[i];
244          if (IsParameterFixed(i) )
245             os << " \t(fixed)";
246          else {
247             if (fErrors.size() != 0)
248                os << " \t+/-\t" << fErrors[i];
249             if (IsParameterBound(i) )
250                os << " \t (limited)";
251          }
252          os << std::endl;
253     }     }
254    
255     if (doCovMatrix) PrintCovMatrix(os);     if (doCovMatrix) PrintCovMatrix(os);
# Line 146  Line 258 
258  void FitResult::PrintCovMatrix(std::ostream &os) const {  void FitResult::PrintCovMatrix(std::ostream &os) const {
259     // print the covariance and correlation matrix     // print the covariance and correlation matrix
260     if (!fValid) return;     if (!fValid) return;
261     os << "\n****************************************\n";     if (fCovMatrix.size() == 0) return;
262    //   os << "****************************************\n";
263     os << "\n            Covariance Matrix            \n\n";     os << "\n            Covariance Matrix            \n\n";
264     unsigned int npar = fParams.size();     unsigned int npar = fParams.size();
265     const int kPrec = 8;     const int kPrec = 5;
266     const int kWidth = 12;     const int kWidth = 10;
267       const int parw = 12;
268       const int matw = kWidth+4;
269       os << std::setw(parw) << " " << "\t";
270       for (unsigned int i = 0; i < npar; ++i)
271          if (!IsParameterFixed(i) ) {
272             os << std::setw(matw)  << fFitFunc->ParameterName(i) ;
273          }
274       os << std::endl;
275     for (unsigned int i = 0; i < npar; ++i) {     for (unsigned int i = 0; i < npar; ++i) {
276        for (unsigned int j = 0; j <= i; ++j) {        if (!IsParameterFixed(i) ) {
277           os.precision(kPrec); os.width(kWidth);  os << CovMatrix(i,j);           os << std::setw(parw) << std::left << fFitFunc->ParameterName(i) << "\t";
278             for (unsigned int j = 0; j < npar; ++j) {
279                if (!IsParameterFixed(j) ) {
280                   os.precision(kPrec); os.width(kWidth);  os << std::setw(matw) << CovMatrix(i,j);
281                }
282        }        }
283        os << std::endl;        os << std::endl;
284     }     }
285       }
286    //   os << "****************************************\n";
287     os << "\n            Correlation Matrix         \n\n";     os << "\n            Correlation Matrix         \n\n";
288       os << std::setw(parw) << " " << "\t";
289       for (unsigned int i = 0; i < npar; ++i)
290          if (!IsParameterFixed(i) ) {
291             os << std::setw(matw)  << fFitFunc->ParameterName(i) ;
292          }
293       os << std::endl;
294     for (unsigned int i = 0; i < npar; ++i) {     for (unsigned int i = 0; i < npar; ++i) {
295        for (unsigned int j = 0; j <= i; ++j) {        if (!IsParameterFixed(i) ) {
296           os.precision(kPrec); os.width(kWidth);  os << Correlation(i,j);           os << std::setw(parw) << std::left << fFitFunc->ParameterName(i) << "\t";
297             for (unsigned int j = 0; j < npar; ++j) {
298                if (!IsParameterFixed(j) ) {
299                   os.precision(kPrec); os.width(kWidth);  os << std::setw(matw) << Correlation(i,j);
300                }
301        }        }
302        os << std::endl;        os << std::endl;
303     }     }
304  }  }
305    }
306    
307    void FitResult::GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double * x, double * ci, double cl ) const {
308       // stride1 stride in coordinate  stride2 stride in dimension space
309       // i.e. i-th point in k-dimension is x[ stride1 * i + stride2 * k]
310       // compute the confidence interval of the fit on the given data points
311       // the dimension of the data points must match the dimension of the fit function
312       // confidence intervals are returned in array ci
313    
314       // use student quantile
315       //double t = - TMath::StudentQuantile((1.-cl)/2, f->GetNDF());
316       double t = TMath::StudentQuantile(0.5 + cl/2, fNdf);
317       double chidf = TMath::Sqrt(fChi2/fNdf);
318    
319       if (!fFitFunc) {
320          MATH_ERROR_MSG("FitResult::GetConfidenceIntervals","cannot compute Confidence Intervals without fitter function");
321          return;
322       }
323    
324       unsigned int ndim = fFitFunc->NDim();
325       unsigned int npar = fFitFunc->NPar();
326    
327       std::vector<double> xpoint(ndim);
328       std::vector<double> grad(npar);
329       std::vector<double> vsum(npar);
330    
331       // loop on the points
332       for (unsigned int ipoint = 0; ipoint < n; ++ipoint) {
333    
334          for (unsigned int kdim = 0; kdim < ndim; ++kdim) {
335             unsigned int i = ipoint * stride1 + kdim * stride2;
336             assert(i < ndim*n);
337             xpoint[kdim] = x[ipoint * stride1 + kdim * stride2];
338          }
339    
340          // calculate gradient of fitted function w.r.t the parameters
341    
342          // check first if fFitFunction provides parameter gradient or not
343    
344          // does not provide gradient
345          // t.b.d : skip calculation for fixed parameters
346          ROOT::Math::RichardsonDerivator d;
347          for (unsigned int ipar = 0; ipar < npar; ++ipar) {
348             ROOT::Math::OneDimParamFunctionAdapter<const ROOT::Math::IParamMultiFunction &> fadapter(*fFitFunc,&xpoint.front(),&fParams.front(),ipar);
349             d.SetFunction(fadapter);
350             grad[ipar] = d(fParams[ipar] ); // evaluate df/dp
351          }
352    
353          // multiply covariance matrix with gradient
354          vsum.assign(npar,0.0);
355          for (unsigned int ipar = 0; ipar < npar; ++ipar) {
356             for (unsigned int jpar = 0; jpar < npar; ++jpar) {
357                vsum[ipar] += CovMatrix(ipar,jpar) * grad[jpar];
358             }
359          }
360          // multiply gradient by vsum
361          double r2 = 0;
362          for (unsigned int ipar = 0; ipar < npar; ++ipar) {
363             r2 += grad[ipar] * vsum[ipar];
364          }
365          double r = std::sqrt(r2);
366          ci[ipoint] = r * t * chidf;
367       }
368    }
369    
370    void FitResult::GetConfidenceIntervals(const BinData & data, double * ci, double cl ) const {
371       // implement confidence intervals from a given bin data sets
372       // currently copy the data from Bindata.
373       // could implement otherwise directly
374       unsigned int ndim = data.NDim();
375       unsigned int np = data.NPoints();
376       std::vector<double> xdata( ndim * np );
377       for (unsigned int i = 0; i < np ; ++i) {
378          const double * x = data.Coords(i);
379          std::vector<double>::iterator itr = xdata.begin()+ ndim * i;
380          std::copy(x,x+ndim,itr);
381       }
382       // points are arraned as x0,y0,z0, ....xN,yN,zN  (stride1=ndim, stride2=1)
383       GetConfidenceIntervals(np,ndim,1,&xdata.front(),ci,cl);
384    }
385    
386     } // end namespace Fit     } // end namespace Fit
387    
Legend:
Removed from v.25485  
changed lines
  Added in v.25486
Subversion Admin
 ViewVC Help
Powered by ViewVC 1.0.9