root/RooStats/MaxLikelihoodEstimateTestStat.h

0001 // @(#)root/roostats:$Id$
0002 // Author: Kyle Cranmer    June 2010
0003 /*************************************************************************
0004  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
0005  * All rights reserved.                                                  *
0006  *                                                                       *
0007  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0008  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0009  *************************************************************************/
0010
0011 #ifndef ROOSTATS_MaxLikelihoodEstimateTestStat
0012 #define ROOSTATS_MaxLikelihoodEstimateTestStat
0013
0014
0015
0016
0017 #include "Rtypes.h"
0018
0019 #include "RooFitResult.h"
0020 #include "RooStats/TestStatistic.h"
0021 #include "RooAbsPdf.h"
0022 #include "RooRealVar.h"
0023 #include "RooMinimizer.h"
0024 #include "Math/MinimizerOptions.h"
0025 #include "RooStats/RooStatsUtils.h"
0026
0027
0028
0029 namespace RooStats {
0030
0031 /** \class MaxLikelihoodEstimateTestStat
0032     \ingroup Roostats
0033 MaxLikelihoodEstimateTestStat: TestStatistic that returns maximum likelihood
0034 estimate of a specified parameter.
0035 */
0036
0037 class MaxLikelihoodEstimateTestStat : public TestStatistic {
0038
0039 public:
0040    MaxLikelihoodEstimateTestStat()
0041       : fStrategy(::ROOT::Math::MinimizerOptions::DefaultStrategy()),
0042         fPrintLevel(::ROOT::Math::MinimizerOptions::DefaultPrintLevel())
0043    {
0044    }
0045
0046    MaxLikelihoodEstimateTestStat(RooAbsPdf &pdf, RooRealVar &parameter)
0047       : fPdf(&pdf),
0048         fParameter(&parameter),
0049         fStrategy(::ROOT::Math::MinimizerOptions::DefaultStrategy()),
0050         fPrintLevel(::ROOT::Math::MinimizerOptions::DefaultPrintLevel())
0051    {
0052    }
0053
0054   //______________________________
0055   double Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/) override {
0056
0057
0058     RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
0059     RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
0060
0061     /*
0062     // this is more straight forward, but produces a lot of messages
0063     RooFitResult* res = fPdf.fitTo(data, RooFit::CloneData(false),RooFit::Minos(0),RooFit::Hesse(false), RooFit::Save(1),RooFit::PrintLevel(-1),RooFit::PrintEvalErrors(0));
0064     RooRealVar* mle = (RooRealVar*) res->floatParsFinal().find(fParameter.GetName());
0065     double ret = mle->getVal();
0066     delete res;
0067     return ret;
0068     */
0069
0070     std::unique_ptr<RooArgSet> allParams{fPdf->getParameters(data)};
0071     RooStats::RemoveConstantParameters(&*allParams);
0072
0073     // need to call constrain for RooSimultaneous until stripDisconnected problem fixed
0074     std::unique_ptr<RooAbsReal> nll{fPdf->createNLL(data, RooFit::CloneData(false),RooFit::Constrain(*allParams),RooFit::ConditionalObservables(fConditionalObs))};
0075
0076     //RooAbsReal* nll = fPdf->createNLL(data, RooFit::CloneData(false));
0077
0078     // RooAbsReal* profile = nll->createProfile(RooArgSet());
0079     // profile->getVal();
0080     // RooArgSet* vars = profile->getVariables();
0081     // RooMsgService::instance().setGlobalKillBelow(msglevel);
0082     // double ret = vars->getRealValue(fParameter->GetName());
0083     // return ret;
0084
0085
0086      RooMinimizer minim(*nll);
0087      minim.setStrategy(fStrategy);
0088      //LM: RooMinimizer.setPrintLevel has +1 offset - so subtract  here -1
0089      minim.setPrintLevel(fPrintLevel-1);
0090      int status = -1;
0091      //   minim.optimizeConst(true);
0092      for (int tries = 0, maxtries = 4; tries <= maxtries; ++tries) {
0093      //    status = minim.minimize(fMinimizer, ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str());
0094         status = minim.minimize(fMinimizer, "Minimize");
0095         if (status == 0) {
0096            break;
0097         } else {
0098            if (tries > 1) {
0099          std::cout << "    ----> Doing a re-scan first\n";
0100          minim.minimize(fMinimizer,"Scan");
0101        }
0102            if (tries > 2) {
0103               std::cout << "    ----> trying with strategy = 1\n";
0104               minim.setStrategy(1);
0105            }
0106         }
0107      }
0108      //std::cout << "BEST FIT values " << std::endl;
0109      //allParams->Print("V");
0110
0111      RooMsgService::instance().setGlobalKillBelow(msglevel);
0112
0113      if (status != 0) return -1;
0114      return fParameter->getVal();
0115
0116
0117   }
0118
0119   const TString GetVarName() const override {
0120     TString varName = Form("Maximum Likelihood Estimate of %s",fParameter->GetName());
0121     return varName;
0122   }
0123
0124
0125   virtual void PValueIsRightTail(bool isright) {  fUpperLimit = isright; }
0126   bool PValueIsRightTail(void) const override { return fUpperLimit; }
0127
0128    // set the conditional observables which will be used when creating the NLL
0129    // so the pdf's will not be normalized on the conditional observables when computing the NLL
0130    void SetConditionalObservables(const RooArgSet& set) override {fConditionalObs.removeAll(); fConditionalObs.add(set);}
0131
0132
0133    private:
0134       RooAbsPdf *fPdf = nullptr;
0135       RooRealVar *fParameter = nullptr;
0136       RooArgSet fConditionalObs;
0137       bool fUpperLimit = true;
0138       TString fMinimizer;
0139       Int_t fStrategy;
0140       Int_t fPrintLevel;
0141
0142
0143
0144    protected:
0145    ClassDefOverride(MaxLikelihoodEstimateTestStat,2)
0146 };
0147
0148 }
0149
0150
0151 #endif