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

Diff of /trunk/math/mathcore/src/FitUtil.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 21  Line 21 
21  #include "Math/IntegratorMultiDim.h"  #include "Math/IntegratorMultiDim.h"
22    
23  #include "Math/Error.h"  #include "Math/Error.h"
24    #include "Math/Util.h"  // for safe log(x)
25    
26  #include <limits>  #include <limits>
27  #include <cmath>  #include <cmath>
# Line 114  Line 115 
115              // construct from function and gradient dimension gdim              // construct from function and gradient dimension gdim
116              // gdim = npar for parameter gradient              // gdim = npar for parameter gradient
117              // gdim = ndim for coordinate gradients              // gdim = ndim for coordinate gradients
118                // construct (the param values will be passed later)
119              // construct for parameter gradient calculation (the param values will be passed later)              // one can choose between 2 points rule (1 extra evaluation) istrat=1
120              // dimension is npar              // or two point rule (2 extra evaluation)
121              SimpleGradientCalculator(const IModelFunction & func) :              // (found 2 points rule does not work correctly - minuit2FitBench fails)
122                 fEps(1.0E-4),              SimpleGradientCalculator(int gdim, const IModelFunction & func,double eps = 2.E-8, int istrat = 1) :
123                 fN(func.NPar() ),                 fEps(eps),
124                   fPrecision(1.E-8 ), // sqrt(epsilon)
125                   fStrategy(istrat),
126                   fN(gdim ),
127                 fFunc(func),                 fFunc(func),
128                 fVec(std::vector<double>(fN) )                 fVec(std::vector<double>(gdim) ) // this can be probably optimized
129              {}              {}
130    
             // construct for coordinate gradient calculator  
             // dimension is ndim  
             SimpleGradientCalculator(const IModelFunction & func, const double * par) :  
                fEps(1.0E-4),  
                fN(func.NDim() ),  
                fFunc(func),  
                fVec(std::vector<double>(fN) ),  
                fPar(std::vector<double>(par, par + func.NPar() ) )  
             {}  
131    
132              // calculate gradient at point (x,p) knnowing already value f0 (we gain a function eval.)              // calculate gradient at point (x,p) knnowing already value f0 (we gain a function eval.)
133              void ParameterGradient(const double * x, const double * p, double f0, double * g) {              void ParameterGradient(const double * x, const double * p, double f0, double * g) {
134                 // fVec are the cached parameter values                 // fVec are the cached parameter values
135                 std::copy(p, p+fN, fVec.begin());                 std::copy(p, p+fN, fVec.begin());
136                 for (unsigned int k = 0; k < fN; ++k) {                 for (unsigned int k = 0; k < fN; ++k) {
137                    fVec[k] += fEps;                    double p0 = p[k];
138                      double h = std::max( fEps* std::abs(p0), 8.0*fPrecision*(std::abs(p0) + fPrecision) );
139                      fVec[k] += h;
140                    // t.b.d : treat case of infinities                    // t.b.d : treat case of infinities
141                    //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )                    //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
142                      double f1 = fFunc(x, &fVec.front() );
143                      if (fStrategy > 1) {
144                         fVec[k] = p0 - h;
145                    double f2 = fFunc(x, &fVec.front() );                    double f2 = fFunc(x, &fVec.front() );
146  //         double f2 = invError * ( y - func(x,&p2.front()) );                       g[k] = 0.5 * ( f2 - f1 )/h;
147                    g[k] = ( f2 - f0 )/fEps;                    }
148                      else
149                         g[k] = ( f1 - f0 )/h;
150    
151                    fVec[k] = p[k]; // restore original p value                    fVec[k] = p[k]; // restore original p value
152                 }                 }
153              }              }
154    
155              // calculate gradient w.r coordinate values              // calculate gradient w.r coordinate values
156              void Gradient(const double * x, double f0, double * g) {              void Gradient(const double * x, const double * p, double f0, double * g) {
157                 // fVec are the cached coordinate values                 // fVec are the cached coordinate values
158                 std::copy(x, x+fN, fVec.begin());                 std::copy(x, x+fN, fVec.begin());
159                 for (unsigned int k = 0; k < fN; ++k) {                 for (unsigned int k = 0; k < fN; ++k) {
160                    fVec[k] += fEps;                    double x0 = x[k];
161                      double h = std::max( fEps* std::abs(x0), 8.0*fPrecision*(std::abs(x0) + fPrecision) );
162                      fVec[k] += h;
163                    // t.b.d : treat case of infinities                    // t.b.d : treat case of infinities
164                    //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )                    //if (fval > - std::numeric_limits<double>::max() && fval < std::numeric_limits<double>::max() )
165                    double f2 = fFunc( &fVec.front(), &fPar.front() );                    double f1 = fFunc( &fVec.front(), p );
166                    g[k] = ( f2 - f0 )/fEps;                    if (fStrategy > 1) {
167                         fVec[k] = x0 - h;
168                         double f2 = fFunc( &fVec.front(), p  );
169                         g[k] = 0.5 * ( f2 - f1 )/h;
170                      }
171                      else
172                         g[k] = ( f1 - f0 )/h;
173    
174                    fVec[k] = x[k]; // restore original x value                    fVec[k] = x[k]; // restore original x value
175                 }                 }
176              }              }
# Line 166  Line 178 
178           private:           private:
179    
180              double fEps;              double fEps;
181                double fPrecision;
182                int fStrategy; // strategy in calculation ( =1 use 2 point rule( 1 extra func) , = 2 use r point rule)
183              unsigned int fN; // gradient dimension              unsigned int fN; // gradient dimension
184              const IModelFunction & fFunc;              const IModelFunction & fFunc;
185              std::vector<double> fVec; // cached coordinates (or parameter values in case of gradientpar)              std::vector<double> fVec; // cached coordinates (or parameter values in case of gradientpar)
             std::vector<double> fPar; // cached parameter values  
186           };           };
187    
188    
# Line 230  Line 243 
243    
244     // get fit option and check case if using integral of bins     // get fit option and check case if using integral of bins
245     const DataOptions & fitOpt = data.Opt();     const DataOptions & fitOpt = data.Opt();
246     IntegralEvaluator igEval( func, p, fitOpt.fIntegral);     bool useBinIntegral = fitOpt.fIntegral;
247       IntegralEvaluator igEval( func, p, useBinIntegral);
248    
249    
250     for (unsigned int i = 0; i < n; ++ i) {     for (unsigned int i = 0; i < n; ++ i) {
251    
 //#define OLD  
 #ifdef OLD  
       const double * x = data.Coords(i);  
       double y = data.Value(i);  
       double invError = data.InvError(i);  
 #else  
252    
253        double y, invError;        double y, invError;
254        const double * x = 0;        // in case of no error in y invError=1 is returned
255        if (!fitOpt.fErrors1)        const double * x = data.GetPoint(i,y, invError);
          x = data.GetPoint(i,y, invError);  
       else {  
          x = data.GetPoint(i,y);  
          invError = 1.;  
       }  
   
 #endif  
256    
257    
258        double fval = 0;        double fval = 0;
259    
260        if (!fitOpt.fIntegral )        if (!useBinIntegral )
261           fval = func ( x, p );           fval = func ( x, p );
262        else {        else {
263           // calculate normalized integral (divided by bin volume)           // calculate normalized integral (divided by bin volume)
# Line 267  Line 268 
268    
269  #ifdef DEBUG  #ifdef DEBUG
270        std::cout << x[0] << "  " << y << "  " << 1./invError << " params : ";        std::cout << x[0] << "  " << y << "  " << 1./invError << " params : ";
271        for (int ipar = 0; ipar < func.NPar(); ++ipar)        for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
272           std::cout << p[ipar] << "\t";           std::cout << p[ipar] << "\t";
273        std::cout << "\tfval = " << fval << std::endl;        std::cout << "\tfval = " << fval << std::endl;
274  #endif  #endif
# Line 276  Line 277 
277        double tmp = ( y -fval )* invError;        double tmp = ( y -fval )* invError;
278        double resval = tmp * tmp;        double resval = tmp * tmp;
279    
280    
281          // avoid inifinity or nan in chi2 values due to wrong function values
282        if ( resval < std::numeric_limits<double>::max() )        if ( resval < std::numeric_limits<double>::max() )
283           chi2 += resval;           chi2 += resval;
284        else        else
285           nRejected++;           nRejected++;
286    
287    
288     }     }
289    
290     // reset the number of fitting data points     // reset the number of fitting data points
# Line 307  Line 311 
311    
312  #ifdef DEBUG  #ifdef DEBUG
313     std::cout << "\n\nFit data size = " << n << std::endl;     std::cout << "\n\nFit data size = " << n << std::endl;
314     std::cout << "evaluate chi2 using function " << &func << "  " << p << std::endl;     std::cout << "evaluate effective chi2 using function " << &func << "  " << p << std::endl;
315  #endif  #endif
316    
317     assert(data.HaveCoordErrors() );     assert(data.HaveCoordErrors() );
# Line 319  Line 323 
323     //func.SetParameters(p);     //func.SetParameters(p);
324    
325     unsigned int ndim = func.NDim();     unsigned int ndim = func.NDim();
326     SimpleGradientCalculator gradCalc(func, p );     SimpleGradientCalculator gradCalc(ndim, func );
327     std::vector<double> grad( ndim );     std::vector<double> grad( ndim );
328    
329    
# Line 354  Line 358 
358        if (j < ndim) {        if (j < ndim) {
359           for (unsigned int icoord = 0; icoord < ndim; ++icoord) {           for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
360              if (ex[icoord] != 0)              if (ex[icoord] != 0)
361                 gradCalc.Gradient(x, fval, &grad[0]);                 gradCalc.Gradient(x, p, fval, &grad[0]);
362              double edx = ex[icoord] * grad[icoord];              double edx = ex[icoord] * grad[icoord];
363              e2 += edx * edx;              e2 += edx * edx;
364           }           }
# Line 363  Line 367 
367        double resval = w2 * ( y - fval ) *  ( y - fval);        double resval = w2 * ( y - fval ) *  ( y - fval);
368    
369  #ifdef DEBUG  #ifdef DEBUG
370        std::cout << x[0] << "  " << y << "  " << e2 << " ey  " << ey << " params : ";        std::cout << x[0] << "  " << y << " ex " << e2 << " ey  " << ey << " params : ";
371        for (int ipar = 0; ipar < func.NPar(); ++ipar)        for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
372           std::cout << p[ipar] << "\t";           std::cout << p[ipar] << "\t";
373        std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;        std::cout << "\tfval = " << fval << "\tresval = " << resval << std::endl;
374  #endif  #endif
# Line 394  Line 398 
398  double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {  double FitUtil::EvaluateChi2Residual(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g) {
399     // evaluate the chi2 contribution (residual term) only for data with no coord-errors     // evaluate the chi2 contribution (residual term) only for data with no coord-errors
400     // This function is used in the specialized least square algorithms like FUMILI or L.M.     // This function is used in the specialized least square algorithms like FUMILI or L.M.
401     // if we have error on the coordinates method is not implemented     // if we have error on the coordinates the method is not yet implemented
402     //  integral option is also not yet implemented     //  integral option is also not yet implemented
403     //  one can use in that case normal chi2 method     //  one can use in that case normal chi2 method
404    
405     if (data.Opt().fCoordErrors )     if (data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
406        MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");        MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 residual");
407          return 0; // it will assert otherwise later in GetPoint
408       }
409    
410     if (data.Opt().fIntegral )     if (data.Opt().fIntegral )
411        MATH_WARN_MSG("FitUtil::EvaluateChi2Residual","Bin integral is not used in calculating Chi2 residual");        MATH_WARN_MSG("FitUtil::EvaluateChi2Residual","Bin integral is not used in calculating Chi2 residual");
# Line 432  Line 438 
438        gfunc->ParameterGradient(  x , p, g );        gfunc->ParameterGradient(  x , p, g );
439     }     }
440     else {     else {
441        SimpleGradientCalculator  gc(func);        SimpleGradientCalculator  gc( func.NPar(), func);
442        gc.ParameterGradient(x, p, fval, g);        gc.ParameterGradient(x, p, fval, g);
443     }     }
444     // mutiply by - 1 * weihgt     // mutiply by - 1 * weihgt
# Line 445  Line 451 
451    
452  }  }
453    
   
454  void FitUtil::EvaluateChi2Gradient(const IModelFunction & f, const BinData & data, const double * p, double * grad, unsigned int & nPoints) {  void FitUtil::EvaluateChi2Gradient(const IModelFunction & f, const BinData & data, const double * p, double * grad, unsigned int & nPoints) {
455     // evaluate gradient of chi2     // evaluate gradient of chi2
456     // this function is used when the model function knows how to calculate the derivative and we can     // this function is used when the model function knows how to calculate the derivative and we can
# Line 456  Line 461 
461     if (data.Opt().fIntegral )     if (data.Opt().fIntegral )
462        MATH_WARN_MSG("FitUtil::EvaluateChi2Residual","Bin integral is not used in calculating Chi2 gradient");        MATH_WARN_MSG("FitUtil::EvaluateChi2Residual","Bin integral is not used in calculating Chi2 gradient");
463    
464     if (data.Opt().fCoordErrors )     if ( data.GetErrorType() == BinData::kCoordError && data.Opt().fCoordErrors ) {
465        MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 gradient");        MATH_ERROR_MSG("FitUtil::EvaluateChi2Residual","Error on the coordinates are not used in calculating Chi2 gradient");            return; // it will assert otherwise later in GetPoint
466       }
467    
468     unsigned int nRejected = 0;     unsigned int nRejected = 0;
469    
# Line 475  Line 481 
481    
482     //int nRejected = 0;     //int nRejected = 0;
483     // set values of parameters     // set values of parameters
484     //func.SetParameters(p);  
485     unsigned int npar = func.NPar();     unsigned int npar = func.NPar();
486  //   assert (npar == NDim() );  // npar MUST be  Chi2 dimension  //   assert (npar == NDim() );  // npar MUST be  Chi2 dimension
487     std::vector<double> gradFunc( npar );     std::vector<double> gradFunc( npar );
# Line 488  Line 494 
494        double y, invError = 0;        double y, invError = 0;
495        const double * x = data.GetPoint(i,y, invError);        const double * x = data.GetPoint(i,y, invError);
496    
497        double fval = func ( x );        double fval = func ( x, p );
498        func.ParameterGradient(  x , p, &gradFunc[0] );        func.ParameterGradient(  x , p, &gradFunc[0] );
499    
500  #ifdef DEBUG  #ifdef DEBUG
501        std::cout << x[0] << "  " << y << "  " << 1./invError << " params : ";        std::cout << x[0] << "  " << y << "  " << 1./invError << " params : ";
502        for (int ipar = 0; ipar < npar; ++ipar)        for (unsigned int ipar = 0; ipar < npar; ++ipar)
503           std::cout << p[ipar] << "\t";           std::cout << p[ipar] << "\t";
504        std::cout << "\tfval = " << fval << std::endl;        std::cout << "\tfval = " << fval << std::endl;
505  #endif  #endif
# Line 529  Line 535 
535    
536  }  }
537    
   
538  //______________________________________________________________________________________________________  //______________________________________________________________________________________________________
539  //  //
540  //  Log Likelihood functions  //  Log Likelihood functions
# Line 537  Line 542 
542    
543  // utility function used by the likelihoods  // utility function used by the likelihoods
544    
 inline double EvalLogF(double fval) {  
    // evaluate the log with a protections against negative argument to the log  
    // smooth linear extrapolation below function values smaller than  epsilon  
    // (better than a simple cut-off)  
    const static double epsilon = 2.*std::numeric_limits<double>::min();  
    if(fval<= epsilon)  
       return fval/epsilon + std::log(epsilon) - 1;  
    else  
       return std::log(fval);  
 }  
   
545  // for LogLikelihood functions  // for LogLikelihood functions
546    
547  double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {  double FitUtil::EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int i, double * g) {
548     // evaluate the pdf contribution to the logl     // evaluate the pdf contribution to the logl
549       // return actually the log of the pdf and its derivatives
550    
551     //func.SetParameters(p);     //func.SetParameters(p);
552    
553     const double * x = data.Coords(i);     const double * x = data.Coords(i);
554     double fval = func ( x );     double fval = func ( x, p );
555     //return EvalLogF(fval);     double logPdf = ROOT::Math::Util::EvalLog(fval);
556     if (g == 0) return fval;     //return
557       if (g == 0) return logPdf;
558    
559     const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);     const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
560    
# Line 569  Line 565 
565     }     }
566     else {     else {
567        // estimate gradieant numerically with simple 2 point rule        // estimate gradieant numerically with simple 2 point rule
568        SimpleGradientCalculator gc(func);        // should probably calculate gradient of log(pdf) is more stable numerically
569          SimpleGradientCalculator gc(func.NPar(), func);
570        gc.ParameterGradient(x, p, fval, g );        gc.ParameterGradient(x, p, fval, g );
571     }     }
572       // divide gradient by function value since returning the logs
573       for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar) {
574          g[ipar] /= fval; // this should be checked against infinities
575       }
576    
577  #ifdef DEBUG  #ifdef DEBUG
578     std::cout << x[i] << "\t";     std::cout << x[i] << "\t";
579     std::cout << "\tpar = [ " << func.NPar() << " ] =  ";     std::cout << "\tpar = [ " << func.NPar() << " ] =  ";
580     for (int ipar = 0; ipar < func.NPar(); ++ipar)     for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
581        std::cout << p[ipar] << "\t";        std::cout << p[ipar] << "\t";
582     std::cout << "\tfval = " << fval;     std::cout << "\tfval = " << fval;
583     std::cout << "\tgrad = [ ";     std::cout << "\tgrad = [ ";
584     for (int ipar = 0; ipar < func.NPar(); ++ipar)     for (unsigned int ipar = 0; ipar < func.NPar(); ++ipar)
585        std::cout << g[ipar] << "\t";        std::cout << g[ipar] << "\t";
586     std::cout << " ] "   << std::endl;     std::cout << " ] "   << std::endl;
587  #endif  #endif
588    
589    
590     return fval;     return logPdf;
591  }  }
592    
593  double FitUtil::EvaluateLogL(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int &nPoints) {  double FitUtil::EvaluateLogL(const IModelFunction & func, const UnBinData & data, const double * p, unsigned int &nPoints) {
# Line 619  Line 620 
620           nRejected++; // reject points with negative pdf (cannot exist)           nRejected++; // reject points with negative pdf (cannot exist)
621        }        }
622        else        else
623           logl += EvalLogF( fval);           logl += ROOT::Math::Util::EvalLog( fval);
624    
625     }     }
626    
# Line 665  Line 666 
666  // for binned log likelihood functions  // for binned log likelihood functions
667  //------------------------------------------------------------------------------------------------  //------------------------------------------------------------------------------------------------
668  double FitUtil::EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g ) {  double FitUtil::EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * p, unsigned int i, double * g ) {
669     // evaluate the pdf contribution to the logl     // evaluate the pdf contribution to the logl (return actually log of pdf)
670     // t.b.d. implement integral option     // t.b.d. implement integral option
671    
672    
# Line 679  Line 680 
680    
681     double fval = func ( x , p);     double fval = func ( x , p);
682    
683     // remove constant term depending on N     // logPdf for Poisson: ignore constant term depending on N
684     double logPdf =   y * EvalLogF( fval) - fval;     double logPdf =   y * ROOT::Math::Util::EvalLog( fval) - fval;
685     // need to return the pdf contribution (not the log)     // need to return the pdf contribution (not the log)
686    
687     double pdfval =  std::exp(logPdf);     //double pdfval =  std::exp(logPdf);
688    
689     if (g == 0) return pdfval;    //if (g == 0) return pdfval;
690       if (g == 0) return logPdf;
691    
692     unsigned int npar = func.NPar();     unsigned int npar = func.NPar();
693     const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);     const IGradModelFunction * gfunc = dynamic_cast<const IGradModelFunction *>( &func);
# Line 696  Line 698 
698        gfunc->ParameterGradient(  x , p, g );        gfunc->ParameterGradient(  x , p, g );
699        // correct g[] do be derivative of poisson term        // correct g[] do be derivative of poisson term
700        for (unsigned int k = 0; k < npar; ++k) {        for (unsigned int k = 0; k < npar; ++k) {
701           g[k] *= ( y/fval - 1.) * pdfval;           g[k] *= ( y/fval - 1.) ;//* pdfval;
702        }        }
703     }     }
704     else {     else {
705        // estimate gradient numerically with simple 2 point rule  //       // estimate gradient numerically with simple 2 point rule
706        static const double kEps = 1.0E-4;  //       static const double kEps = 1.0E-6;
707        std::vector<double> p2(p,p+npar);  //       std::vector<double> p2(p,p+npar);
708        for (unsigned int k = 0; k < npar; ++k) {  //       for (unsigned int k = 0; k < npar; ++k) {
709           p2[k] += kEps;  //          p2[k] += kEps*p2[k];
710           // t.b.d : treat case of infinities  //          // t.b.d : treat case of infinities
711           // make derivatives of the log (more stables ) and correct afterwards  //          // make derivatives of the log (more stables ) and correct afterwards
712           double fval2 = func(x,&p2.front() );  //          double fval2 = func(x,&p2.front() );
713           double logPdf2 = y * EvalLogF( fval2) - fval2;  //          double logPdf2 = y * ROOT::Math::Util::EvalLog( fval2) - fval2;
714           double dlogPdf = (logPdf2 - logPdf)/kEps;  //          double dlogPdf = (logPdf2 - logPdf)/kEps;
715           g[k] = dlogPdf * pdfval;  //          g[k] = dlogPdf;// * pdfval;
716           p2[k] = p[k]; // restore original p value  //          p2[k] = p[k]; // restore original p value
717    
718          SimpleGradientCalculator  gc(func.NPar(), func);
719          gc.ParameterGradient(x, p, fval, g);
720          // correct g[] do be derivative of poisson term
721          for (unsigned int k = 0; k < npar; ++k) {
722             g[k] *= ( y/fval - 1.) ;
723        }        }
724    
725     }     }
726    
727     return pdfval;  #ifdef DEBUG
728       std::cout << "x = " << x[0] << " logPdf = " << logPdf << " grad";
729       for (unsigned int ipar = 0; ipar < npar; ++ipar)
730          std::cout << g[ipar] << "\t";
731       std::cout << std::endl;
732    #endif
733    
734    //   return pdfval;
735       return logPdf;
736  }  }
737    
738  double FitUtil::EvaluatePoissonLogL(const IModelFunction & func, const BinData & data, const double * p, unsigned int &nPoints) {  double FitUtil::EvaluatePoissonLogL(const IModelFunction & func, const BinData & data, const double * p, unsigned int &nPoints) {
# Line 746  Line 762 
762        }        }
763    
764    
765        loglike +=  fval - y * EvalLogF( fval);        loglike +=  fval - y * ROOT::Math::Util::EvalLog( fval);
766    
767    
768     }     }
# Line 755  Line 771 
771     if (nRejected != 0)  nPoints = n - nRejected;     if (nRejected != 0)  nPoints = n - nRejected;
772    
773  #ifdef DEBUG  #ifdef DEBUG
774     std::cout << "Logl = " << logl << " np = " << nPoints << std::endl;     std::cout << "Loglikelihood  = " << loglike << " np = " << nPoints << std::endl;
775  #endif  #endif
776    
777     return loglike;     return loglike;
Legend:
Removed from v.25485  
changed lines
  Added in v.25486
Subversion Admin
 ViewVC Help
Powered by ViewVC 1.0.9