root/RooStats/ToyMCSampler.h

0001 // @(#)root/roostats:$Id$
0002 // Author: Sven Kreiss and Kyle Cranmer    June 2010
0003 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
0004 // Additions and modifications by Mario Pelliccioni
0005 /*************************************************************************
0006  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
0007  * All rights reserved.                                                  *
0008  *                                                                       *
0009  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0010  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0011  *************************************************************************/
0012
0013 #ifndef ROOSTATS_ToyMCSampler
0014 #define ROOSTATS_ToyMCSampler
0015
0016 #include "RooStats/TestStatSampler.h"
0017 #include "RooStats/SamplingDistribution.h"
0018 #include "RooStats/TestStatistic.h"
0019 #include "RooStats/ModelConfig.h"
0020
0021 #include "RooWorkspace.h"
0022 #include "RooMsgService.h"
0023 #include "RooAbsPdf.h"
0024 #include "RooRealVar.h"
0025 #include "RooDataSet.h"
0026
0027 #include <vector>
0028 #include <list>
0029 #include <string>
0030 #include <sstream>
0031 #include <memory>
0032
0033 namespace RooStats {
0034
0035   class DetailedOutputAggregator;
0036
0037 class NuisanceParametersSampler {
0038
0039    public:
0040       NuisanceParametersSampler(RooAbsPdf *prior=nullptr, const RooArgSet *parameters=nullptr, Int_t nToys=1000, bool asimov=false) :
0041          fPrior(prior),
0042          fParams(parameters),
0043          fNToys(nToys),
0044          fExpected(asimov),
0045          fIndex(0)
0046       {
0047          if(prior) Refresh();
0048       }
0049       virtual ~NuisanceParametersSampler() = default;
0050
0051       void NextPoint(RooArgSet& nuisPoint, double& weight);
0052
0053    protected:
0054       void Refresh();
0055
0056    private:
0057       RooAbsPdf *fPrior;           // prior for nuisance parameters
0058       const RooArgSet *fParams;    // nuisance parameters
0059       Int_t fNToys;
0060       bool fExpected;
0061
0062       std::unique_ptr<RooAbsData> fPoints;         // generated nuisance parameter points
0063       Int_t fIndex;                // current index in fPoints array
0064 };
0065
0066 class ToyMCSampler: public TestStatSampler {
0067
0068    public:
0069
0070       ToyMCSampler(TestStatistic &ts, Int_t ntoys);
0071       ~ToyMCSampler() override;
0072
0073       static void SetAlwaysUseMultiGen(bool flag);
0074
0075       void SetUseMultiGen(bool flag) { fUseMultiGen = flag ; }
0076
0077       /// main interface
0078       SamplingDistribution* GetSamplingDistribution(RooArgSet& paramPoint) override;
0079       virtual RooDataSet* GetSamplingDistributions(RooArgSet& paramPoint);
0080       virtual RooDataSet* GetSamplingDistributionsSingleWorker(RooArgSet& paramPoint);
0081
0082       virtual SamplingDistribution* AppendSamplingDistribution(
0083          RooArgSet& allParameters,
0084          SamplingDistribution* last,
0085          Int_t additionalMC
0086       );
0087
0088
0089       /// The pdf can be nullptr in which case the density from SetPdf()
0090       /// is used. The snapshot and TestStatistic is also optional.
0091       virtual void AddTestStatistic(TestStatistic* t = nullptr) {
0092          if( t == nullptr ) {
0093             oocoutI(nullptr,InputArguments) << "No test statistic given. Doing nothing." << std::endl;
0094             return;
0095          }
0096
0097          fTestStatistics.push_back( t );
0098       }
0099
0100       /// generates toy data
0101       ///   without weight
0102       virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, RooAbsPdf& pdf) const {
0103          if(fExpectedNuisancePar) oocoutE(nullptr,InputArguments) << "ToyMCSampler: using expected nuisance parameters but ignoring weight. Use GetSamplingDistribution(paramPoint, weight) instead." << std::endl;
0104          double weight;
0105          return GenerateToyData(paramPoint, weight, pdf);
0106       }
0107       virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint) const { return GenerateToyData(paramPoint,*fPdf); }
0108       /// generates toy data
0109       ///   with weight
0110       virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight, RooAbsPdf& pdf) const;
0111       virtual RooAbsData* GenerateToyData(RooArgSet& paramPoint, double& weight) const { return GenerateToyData(paramPoint,weight,*fPdf); }
0112
0113       /// generate global observables
0114       virtual void GenerateGlobalObservables(RooAbsPdf& pdf) const;
0115
0116
0117       /// Main interface to evaluate the test statistic on a dataset
0118       virtual double EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI, int i ) {
0119          return fTestStatistics[i]->Evaluate(data, nullPOI);
0120       }
0121       double EvaluateTestStatistic(RooAbsData& data, RooArgSet& nullPOI) override { return EvaluateTestStatistic( data,nullPOI, 0 ); }
0122       virtual RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi);
0123
0124
0125       virtual TestStatistic* GetTestStatistic(unsigned int i) const {
0126          if( fTestStatistics.size() <= i ) return nullptr;
0127          return fTestStatistics[i];
0128       }
0129       TestStatistic* GetTestStatistic(void) const override { return GetTestStatistic(0); }
0130
0131       double ConfidenceLevel() const override { return 1. - fSize; }
0132       void Initialize(
0133          RooAbsArg& /*testStatistic*/,
0134          RooArgSet& /*paramsOfInterest*/,
0135          RooArgSet& /*nuisanceParameters*/
0136       ) override {}
0137
0138       virtual Int_t GetNToys(void) { return fNToys; }
0139       virtual void SetNToys(const Int_t ntoy) { fNToys = ntoy; }
0140       /// Forces the generation of exactly `n` events even for extended PDFs. Set to 0 to
0141       /// use the Poisson-distributed events from the extended PDF.
0142       virtual void SetNEventsPerToy(const Int_t nevents) {
0143          fNEvents = nevents;
0144       }
0145
0146       // Getter for the number of events per toy
0147       int nEventsPerToy() const {
0148          return fNEvents;
0149       }
0150
0151
0152       /// Set the Pdf, add to the workspace if not already there
0153       void SetParametersForTestStat(const RooArgSet& nullpoi) override {
0154          auto params = std::make_unique<RooArgSet>();
0155          nullpoi.snapshot(*params);
0156          fParametersForTestStat = std::move(params);
0157       }
0158
0159       void SetPdf(RooAbsPdf& pdf) override { fPdf = &pdf; ClearCache(); }
0160
0161       /// How to randomize the prior. Set to nullptr to deactivate randomization.
0162       void SetPriorNuisance(RooAbsPdf* pdf) override {
0163          fPriorNuisance = pdf;
0164          if (fNuisanceParametersSampler) {
0165             delete fNuisanceParametersSampler;
0166             fNuisanceParametersSampler = nullptr;
0167          }
0168       }
0169       /// specify the nuisance parameters (eg. the rest of the parameters)
0170       void SetNuisanceParameters(const RooArgSet& np) override { fNuisancePars = &np; }
0171       /// specify the observables in the dataset (needed to evaluate the test statistic)
0172       void SetObservables(const RooArgSet& o) override { fObservables = &o; }
0173       /// specify the conditional observables
0174       void SetGlobalObservables(const RooArgSet& o) override { fGlobalObservables = &o; }
0175
0176
0177       /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
0178       void SetTestSize(double size) override { fSize = size; }
0179       /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
0180       void SetConfidenceLevel(double cl) override { fSize = 1. - cl; }
0181
0182       /// Set the TestStatistic (want the argument to be a function of the data & parameter points
0183       virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i) {
0184          if( fTestStatistics.size() < i ) {
0185             oocoutE(nullptr,InputArguments) << "Cannot set test statistic for this index." << std::endl;
0186             return;
0187          }
0188          if (fTestStatistics.size() == i) {
0189             fTestStatistics.push_back(testStatistic);
0190          } else {
0191             fTestStatistics[i] = testStatistic;
0192          }
0193       }
0194       void SetTestStatistic(TestStatistic *t) override { return SetTestStatistic(t,0); }
0195
0196       virtual void SetExpectedNuisancePar(bool i = true) { fExpectedNuisancePar = i; }
0197       virtual void SetAsimovNuisancePar(bool i = true) { fExpectedNuisancePar = i; }
0198
0199       /// Checks for sufficient information to do a GetSamplingDistribution(...).
0200       bool CheckConfig(void);
0201
0202       /// control to use bin data generation (=> see RooFit::AllBinned() option)
0203       void SetGenerateBinned(bool binned = true) { fGenerateBinned = binned; }
0204       /// name of the tag for individual components to be generated binned (=> see RooFit::GenBinned() option)
0205       void SetGenerateBinnedTag( const char* binnedTag = "" ) { fGenerateBinnedTag = binnedTag; }
0206       /// set auto binned generation (=> see RooFit::AutoBinned() option)
0207       void SetGenerateAutoBinned( bool autoBinned = true ) { fGenerateAutoBinned = autoBinned; }
0208
0209       /// Set the name of the sampling distribution used for plotting
0210       void SetSamplingDistName(const char* name) override { if(name) fSamplingDistName = name; }
0211       std::string GetSamplingDistName(void) { return fSamplingDistName; }
0212
0213       /// This option forces a maximum number of total toys.
0214       void SetMaxToys(double t) { fMaxToys = t; }
0215
0216       void SetToysLeftTail(double toys, double threshold) {
0217          fToysInTails = toys;
0218          fAdaptiveLowLimit = threshold;
0219          fAdaptiveHighLimit = RooNumber::infinity();
0220       }
0221       void SetToysRightTail(double toys, double threshold) {
0222          fToysInTails = toys;
0223          fAdaptiveHighLimit = threshold;
0224          fAdaptiveLowLimit = -RooNumber::infinity();
0225       }
0226       void SetToysBothTails(double toys, double low_threshold, double high_threshold) {
0227          fToysInTails = toys;
0228          fAdaptiveHighLimit = high_threshold;
0229          fAdaptiveLowLimit = low_threshold;
0230       }
0231
0232       void SetProtoData(const RooDataSet* d) { fProtoData = d; }
0233
0234    protected:
0235
0236       const RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg);
0237
0238       /// helper for GenerateToyData
0239       std::unique_ptr<RooAbsData> Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=nullptr, int forceEvents=0) const;
0240
0241       /// helper method for clearing  the cache
0242       virtual void ClearCache();
0243
0244
0245       /// densities, snapshots, and test statistics to reweight to
0246       RooAbsPdf *fPdf = nullptr; ///< model (can be alt or null)
0247       std::unique_ptr<const RooArgSet> fParametersForTestStat;
0248       std::vector<TestStatistic*> fTestStatistics;
0249
0250       std::string fSamplingDistName; ///< name of the model
0251       RooAbsPdf *fPriorNuisance = nullptr; ///< prior pdf for nuisance parameters
0252       const RooArgSet *fNuisancePars = nullptr;
0253       const RooArgSet *fObservables = nullptr;
0254       const RooArgSet *fGlobalObservables = nullptr;
0255       Int_t fNToys;   ///< number of toys to generate
0256       Int_t fNEvents = 0; ///< number of events per toy (may be ignored depending on settings)
0257       double fSize = 0.05;
0258       bool fExpectedNuisancePar =
0259          false; ///< whether to use expectation values for nuisance parameters (ie Asimov data set)
0260       bool fGenerateBinned = false;
0261       TString fGenerateBinnedTag = "";
0262       bool fGenerateAutoBinned = true;
0263
0264       /// minimum no of toys in tails for adaptive sampling
0265       /// (taking weights into account, therefore double)
0266       /// Default: 0.0 which means no adaptive sampling
0267       double fToysInTails = 0.0;
0268       /// maximum no of toys
0269       /// (taking weights into account, therefore double)
0270       double fMaxToys;
0271       /// tails
0272       double fAdaptiveLowLimit;
0273       double fAdaptiveHighLimit;
0274
0275       const RooDataSet *fProtoData = nullptr; ///< in dev
0276
0277       mutable NuisanceParametersSampler *fNuisanceParametersSampler = nullptr; ///<!
0278
0279       // objects below cache information and are mutable and non-persistent
0280       mutable std::unique_ptr<RooArgSet> _allVars; ///<!
0281       mutable std::vector<RooAbsPdf*> _pdfList;    ///<! We don't own those objects
0282       mutable std::vector<std::unique_ptr<RooArgSet>> _obsList;         ///<!
0283       mutable std::vector<std::unique_ptr<RooAbsPdf::GenSpec>> _gsList; ///<!
0284       mutable std::unique_ptr<RooAbsPdf::GenSpec> _gs1; ///<! GenSpec #1
0285       mutable std::unique_ptr<RooAbsPdf::GenSpec> _gs2; ///<! GenSpec #2
0286       mutable std::unique_ptr<RooAbsPdf::GenSpec> _gs3; ///<! GenSpec #3
0287       mutable std::unique_ptr<RooAbsPdf::GenSpec> _gs4; ///<! GenSpec #4
0288
0289       static bool fgAlwaysUseMultiGen ;  ///< Use PrepareMultiGen always
0290       bool fUseMultiGen = false;         ///< Use PrepareMultiGen?
0291
0292    ClassDefOverride(ToyMCSampler, 0) // A simple implementation of the TestStatSampler interface
0293 };
0294 }
0295
0296
0297 #endif