Likelihood.h

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_Likelihood_h
#define PhysicsTools_Utilities_Likelihood_h
#include "PhysicsTools/Utilities/interface/RootMinuitResultPrinter.h"
#include "PhysicsTools/Utilities/interface/RootMinuitFuncEvaluator.h"
#include <cmath>
#include "TMath.h"

namespace fit {
  template <typename PDF, typename Tuple>
  struct LikelihoodEvaluator {
    static double evaluate(const PDF& pdf, const Tuple& tuple) { return pdf(tuple); }
  };

  // to be generalized for tuples in N variables
  template <typename PDF>
  struct LikelihoodEvaluator<PDF, double> {
    static double evaluate(const PDF& pdf, const double val) { return pdf(val); }
  };

  struct NoExtendedLikelihood {};

  template <typename Sample, typename PDF, typename Yield = NoExtendedLikelihood>
  class Likelihood {
  public:
    Likelihood() {}
    Likelihood(const Sample& sample, PDF& pdf, Yield& yield) : pdf_(&pdf), yield_(&yield), sample_(sample) {}
    double operator()() const { return log(); }
    double log() const {
      double l = -(*yield_)();
      for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
        double p = Evaluator::evaluate(*pdf_, *i);
        l += std::log(p);
      }
      sampleSize_ = sample_.size();
      return l;
    }
    double logNFactorial() const { return std::log(TMath::Factorial(sampleSize_)); }
    double absoluteLog() const { return log() - logNFactorial(); }
    PDF& pdf() { return *pdf_; }
    const PDF& pdf() const { return *pdf_; }
    Yield& yield() { return *yield_; }
    const Yield& yield() const { return *yield_; }
    unsigned int sampleSize() const { return sampleSize_; }

  private:
    typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
    PDF* pdf_;
    Yield* yield_;
    Sample sample_;
    mutable unsigned int sampleSize_ = 0u;
  };

  template <typename Sample, typename PDF>
  class Likelihood<Sample, PDF, NoExtendedLikelihood> {
  public:
    Likelihood() {}
    Likelihood(const Sample& sample, PDF& pdf) : pdf_(&pdf), sample_(sample) {}
    double operator()() const { return log(); }
    double log() const {
      double l = 0;
      for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
        l += std::log(Evaluator::evaluate(*pdf_, *i));
      }
      return l;
    }
    PDF& pdf() { return *pdf_; }
    const PDF& pdf() const { return *pdf_; }

  private:
    typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
    PDF* pdf_;
    Sample sample_;
  };

  template <typename Sample, typename PDF, typename Yield>
  struct RootMinuitResultPrinter<Likelihood<Sample, PDF, Yield> > {
    static void print(double amin, unsigned int numberOfFreeParameters, const Likelihood<Sample, PDF, Yield>& f) {
      std::cout << "-2 log(maximum-likelihood) = " << amin << ", free parameters = " << numberOfFreeParameters
                << std::endl;
    }
  };

  template <typename Sample, typename PDF, typename Yield>
  struct RootMinuitFuncEvaluator<Likelihood<Sample, PDF, Yield> > {
    static double evaluate(const Likelihood<Sample, PDF, Yield>& f) { return -2 * f(); }
  };

}  // namespace fit

#endif