Utilities/interface/Likelihood.h

0001 #ifndef PhysicsTools_Utilities_Likelihood_h
0002 #define PhysicsTools_Utilities_Likelihood_h
0003 #include "PhysicsTools/Utilities/interface/RootMinuitResultPrinter.h"
0004 #include "PhysicsTools/Utilities/interface/RootMinuitFuncEvaluator.h"
0005 #include <cmath>
0006 #include "TMath.h"
0007
0008 namespace fit {
0009   template <typename PDF, typename Tuple>
0010   struct LikelihoodEvaluator {
0011     static double evaluate(const PDF& pdf, const Tuple& tuple) { return pdf(tuple); }
0012   };
0013
0014   // to be generalized for tuples in N variables
0015   template <typename PDF>
0016   struct LikelihoodEvaluator<PDF, double> {
0017     static double evaluate(const PDF& pdf, const double val) { return pdf(val); }
0018   };
0019
0020   struct NoExtendedLikelihood {};
0021
0022   template <typename Sample, typename PDF, typename Yield = NoExtendedLikelihood>
0023   class Likelihood {
0024   public:
0025     Likelihood() {}
0026     Likelihood(const Sample& sample, PDF& pdf, Yield& yield) : pdf_(&pdf), yield_(&yield), sample_(sample) {}
0027     double operator()() const { return log(); }
0028     double log() const {
0029       double l = -(*yield_)();
0030       for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
0031         double p = Evaluator::evaluate(*pdf_, *i);
0032         l += std::log(p);
0033       }
0034       sampleSize_ = sample_.size();
0035       return l;
0036     }
0037     double logNFactorial() const { return std::log(TMath::Factorial(sampleSize_)); }
0038     double absoluteLog() const { return log() - logNFactorial(); }
0039     PDF& pdf() { return *pdf_; }
0040     const PDF& pdf() const { return *pdf_; }
0041     Yield& yield() { return *yield_; }
0042     const Yield& yield() const { return *yield_; }
0043     unsigned int sampleSize() const { return sampleSize_; }
0044
0045   private:
0046     typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
0047     PDF* pdf_;
0048     Yield* yield_;
0049     Sample sample_;
0050     mutable unsigned int sampleSize_ = 0u;
0051   };
0052
0053   template <typename Sample, typename PDF>
0054   class Likelihood<Sample, PDF, NoExtendedLikelihood> {
0055   public:
0056     Likelihood() {}
0057     Likelihood(const Sample& sample, PDF& pdf) : pdf_(&pdf), sample_(sample) {}
0058     double operator()() const { return log(); }
0059     double log() const {
0060       double l = 0;
0061       for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
0062         l += std::log(Evaluator::evaluate(*pdf_, *i));
0063       }
0064       return l;
0065     }
0066     PDF& pdf() { return *pdf_; }
0067     const PDF& pdf() const { return *pdf_; }
0068
0069   private:
0070     typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
0071     PDF* pdf_;
0072     Sample sample_;
0073   };
0074
0075   template <typename Sample, typename PDF, typename Yield>
0076   struct RootMinuitResultPrinter<Likelihood<Sample, PDF, Yield> > {
0077     static void print(double amin, unsigned int numberOfFreeParameters, const Likelihood<Sample, PDF, Yield>& f) {
0078       std::cout << "-2 log(maximum-likelihood) = " << amin << ", free parameters = " << numberOfFreeParameters
0079                 << std::endl;
0080     }
0081   };
0082
0083   template <typename Sample, typename PDF, typename Yield>
0084   struct RootMinuitFuncEvaluator<Likelihood<Sample, PDF, Yield> > {
0085     static double evaluate(const Likelihood<Sample, PDF, Yield>& f) { return -2 * f(); }
0086   };
0087
0088 }  // namespace fit
0089
0090 #endif