root/RooStats/BayesianCalculator.h

0001 // @(#)root/roostats:$Id$
0002 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
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_BayesianCalculator
0012 #define ROOSTATS_BayesianCalculator
0013
0014 #include "TNamed.h"
0015
0016 #include "Math/IFunctionfwd.h"
0017
0018 #include "RooArgSet.h"
0019
0020 #include "RooStats/IntervalCalculator.h"
0021
0022 #include "RooStats/SimpleInterval.h"
0023
0024 class RooAbsData;
0025 class RooAbsPdf;
0026 class RooPlot;
0027 class RooAbsReal;
0028 class TF1;
0029 class TH1;
0030
0031
0032 namespace RooStats {
0033
0034    class ModelConfig;
0035    class SimpleInterval;
0036
0037    class BayesianCalculator : public IntervalCalculator, public TNamed {
0038
0039    public:
0040
0041       /// constructor
0042       BayesianCalculator( );
0043
0044       BayesianCalculator( RooAbsData& data,
0045                           RooAbsPdf& pdf,
0046                           const RooArgSet& POI,
0047                           RooAbsPdf& priorPdf,
0048                           const RooArgSet* nuisanceParameters = nullptr );
0049
0050       BayesianCalculator( RooAbsData& data,
0051                           ModelConfig& model );
0052
0053       /// destructor
0054       ~BayesianCalculator() override;
0055
0056       /// get the plot with option to get it normalized
0057       RooPlot* GetPosteriorPlot(bool norm = false, double precision = 0.01) const;
0058
0059       /// return posterior pdf (object is managed by the user)
0060       RooAbsPdf* GetPosteriorPdf() const;
0061       /// return posterior function (object is managed by the BayesianCalculator class)
0062       RooAbsReal* GetPosteriorFunction() const;
0063
0064       /// return the approximate posterior as histogram (TH1 object). Note the object is managed by the BayesianCalculator class
0065       TH1 * GetPosteriorHistogram() const;
0066
0067       /// compute the interval. By Default a central interval is computed
0068       /// By using SetLeftTileFraction can control if central/ upper/lower interval
0069       /// For shortest interval use SetShortestInterval(true)
0070       SimpleInterval* GetInterval() const override ;
0071
0072       void SetData( RooAbsData & data ) override {
0073          fData = &data;
0074          ClearAll();
0075       }
0076
0077
0078       /// set the model via the ModelConfig
0079       void SetModel( const ModelConfig& model ) override;
0080
0081       /// specify the parameters of interest in the interval
0082       virtual void SetParameters(const RooArgSet& set) { fPOI.removeAll(); fPOI.add(set); }
0083
0084       /// specify the nuisance parameters (eg. the rest of the parameters)
0085       virtual void SetNuisanceParameters(const RooArgSet& set) {fNuisanceParameters.removeAll(); fNuisanceParameters.add(set);}
0086
0087       /// Set only the Prior Pdf
0088       virtual void SetPriorPdf(RooAbsPdf& pdf) { fPriorPdf = &pdf; }
0089
0090       /// set the conditional observables which will be used when creating the NLL
0091       /// so the pdf's will not be normalized on the conditional observables when computing the NLL
0092       virtual void SetConditionalObservables(const RooArgSet& set) {fConditionalObs.removeAll(); fConditionalObs.add(set);}
0093
0094       /// set the global observables which will be used when creating the NLL
0095       /// so the constraint pdf's will be normalized correctly on the global observables when computing the NLL
0096       virtual void SetGlobalObservables(const RooArgSet& set) {fGlobalObs.removeAll(); fGlobalObs.add(set);}
0097
0098       /// set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval)
0099       void SetTestSize( double size ) override {
0100          fSize = size;
0101          fValidInterval = false;
0102       }
0103       /// set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval)
0104       void SetConfidenceLevel( double cl ) override { SetTestSize(1.-cl); }
0105       /// Get the size of the test (eg. rate of Type I error)
0106       double Size() const override { return fSize; }
0107       /// Get the Confidence level for the test
0108       double ConfidenceLevel() const override { return 1.-fSize; }
0109
0110       /// set the fraction of probability content on the left tail
0111       /// Central limits use 0.5 (default case)
0112       /// for upper limits it is 0 and 1 for lower limit
0113       /// For shortest intervals a negative value (i.e. -1) must be given
0114       void SetLeftSideTailFraction(double leftSideFraction )  {fLeftSideFraction = leftSideFraction;}
0115
0116       /// set the Bayesian calculator to compute the shortest interval (default is central interval)
0117       /// to switch off SetLeftSideTailFraction to the right value
0118       void SetShortestInterval() { fLeftSideFraction = -1; }
0119
0120       /// set the precision of the Root Finder
0121       void SetBrfPrecision( double precision ) { fBrfPrecision = precision; }
0122
0123       /// use directly the approximate posterior function obtained by binning it in nbins
0124       /// by default the cdf is used by integrating the posterior
0125       /// if a value of nbin <= 0 the cdf function will be used
0126       void SetScanOfPosterior(int nbin = 100) { fNScanBins = nbin; }
0127
0128       /// set the number of iterations when running a MC integration algorithm
0129       /// If not set use default algorithmic values
0130       /// In case of ToyMC sampling of the nuisance the value is 100
0131       /// In case of using the GSL MCintegrations types the default value is
0132       /// defined in ROOT::Math::IntegratorMultiDimOptions::DefaultNCalls()
0133       virtual void SetNumIters(Int_t numIters)  { fNumIterations = numIters; }
0134
0135       /// set the integration type (possible type are) :
0136       void SetIntegrationType(const char * type);
0137
0138       /// return the mode (most probable value of the posterior function)
0139       double GetMode() const;
0140
0141       // force the nuisance pdf when using the toy mc sampling
0142       void ForceNuisancePdf(RooAbsPdf & pdf) { fNuisancePdf = &pdf; }
0143
0144    protected:
0145
0146       void ClearAll() const;
0147
0148       void ApproximatePosterior() const;
0149
0150       void ComputeIntervalFromApproxPosterior(double c1, double c2) const;
0151
0152       void ComputeIntervalFromCdf(double c1, double c2) const;
0153
0154       void ComputeIntervalUsingRooFit(double c1, double c2) const;
0155
0156       void ComputeShortestInterval() const;
0157
0158    private:
0159
0160       // plan to replace the above: return a SimpleInterval integrating
0161       // over all other parameters except the one specified as argument
0162       // virtual SimpleInterval* GetInterval( RooRealVar* parameter  ) const { return 0; }
0163
0164       RooAbsData* fData;                         ///< data set
0165       RooAbsPdf* fPdf;                           ///< model pdf  (could contain the nuisance pdf as constraint term)
0166       RooArgSet fPOI;                            ///< POI
0167       RooAbsPdf* fPriorPdf;                      ///< prior pdf (typically for the POI)
0168       RooAbsPdf* fNuisancePdf;                   ///< nuisance pdf (needed when using nuisance sampling technique)
0169       RooArgSet fNuisanceParameters;             ///< nuisance parameters
0170       RooArgSet fConditionalObs    ;             ///< conditional observables
0171       RooArgSet fGlobalObs;                      ///< global observables
0172
0173       mutable RooAbsPdf* fProductPdf;            ///< internal pointer to model * prior
0174       mutable std::unique_ptr<RooAbsReal> fLogLike; ///< internal pointer to log likelihood function
0175       mutable RooAbsReal* fLikelihood;           ///< internal pointer to likelihood function
0176       mutable RooAbsReal* fIntegratedLikelihood; ///< integrated likelihood function, i.e - unnormalized posterior function
0177       mutable RooAbsPdf* fPosteriorPdf;          ///< normalized (on the poi) posterior pdf
0178       mutable ROOT::Math::IGenFunction * fPosteriorFunction;   ///< function representing the posterior
0179       mutable TF1 * fApproxPosterior;    ///< TF1 representing the scanned posterior function
0180       mutable double  fLower;          ///< computer lower interval bound
0181       mutable double  fUpper;          ///< upper interval bound
0182       mutable double  fNLLMin;         ///< minimum value of Nll
0183       double fSize;                      ///< size used for getting the interval
0184       double fLeftSideFraction;          ///< fraction of probability content on left side of interval
0185       double fBrfPrecision;              ///< root finder precision
0186       mutable int fNScanBins;            ///< number of bins to scan, if = -1 no scan is done (default)
0187       int fNumIterations;                ///< number of iterations (when using ToyMC)
0188       mutable bool    fValidInterval;
0189
0190       TString fIntegrationType;
0191
0192    protected:
0193
0194       ClassDefOverride(BayesianCalculator,3)  // BayesianCalculator class
0195
0196    };
0197 }
0198
0199 #endif