root/RooStats/ConfidenceBelt.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_ConfidenceBelt
0012 #define RooStats_ConfidenceBelt
0013
0014 #include "RooArgSet.h"
0015 #include "RooAbsData.h"
0016 #include "RooStats/ConfInterval.h"
0017
0018 #include "RooStats/SamplingDistribution.h"
0019
0020 #include "TRef.h"
0021
0022 #include <vector>
0023 #include <map>
0024
0025
0026 namespace RooStats {
0027
0028 ////////////////////////////////////////////////////////////////////////////////
0029
0030   class SamplingSummaryLookup : public TObject {
0031
0032      typedef std::pair<double, double> AcceptanceCriteria; // defined by Confidence level, leftside tail probability
0033      typedef std::map<Int_t, AcceptanceCriteria> LookupTable; // map ( Index, ( CL, leftside tail prob) )
0034
0035   public:
0036     SamplingSummaryLookup() {}
0037
0038     void Add(double cl, double leftside){
0039       // add cl,leftside pair to lookup table
0040       AcceptanceCriteria tmp(cl, leftside);
0041
0042       // should check to see if this is already in the map
0043       if(GetLookupIndex(cl,leftside) >=0 ){
0044          std::cout<< "SamplingSummaryLookup::Add, already in lookup table" << std::endl;
0045       } else
0046    fLookupTable[fLookupTable.size()]= tmp;
0047     }
0048
0049     Int_t GetLookupIndex(double cl, double leftside){
0050       // get index for cl,leftside pair
0051       AcceptanceCriteria tmp(cl, leftside);
0052
0053       double tolerance = 1E-6; // some small number to protect floating point comparison.  What is better way?
0054       LookupTable::iterator it = fLookupTable.begin();
0055       Int_t index = -1;
0056       for(; it!= fLookupTable.end(); ++it) {
0057    index++;
0058    if( std::abs( (*it).second.first - cl ) < tolerance &&
0059        std::abs( (*it).second.second - leftside ) < tolerance )
0060      break; // exit loop, found
0061       }
0062
0063       // check that it was found
0064       if(index == (Int_t)fLookupTable.size())
0065    index = -1;
0066
0067       return index;
0068     }
0069
0070   double GetConfidenceLevel(Int_t index){
0071     if(index<0 || index>(Int_t)fLookupTable.size()) {
0072        std::cout << "SamplingSummaryLookup::GetConfidenceLevel, index not in lookup table" << std::endl;
0073        return -1;
0074     }
0075     return fLookupTable[index].first;
0076   }
0077
0078   double GetLeftSideTailFraction(Int_t index){
0079     if(index<0 || index>(Int_t)fLookupTable.size()) {
0080        std::cout << "SamplingSummaryLookup::GetLeftSideTailFraction, index not in lookup table" << std::endl;
0081        return -1;
0082     }
0083     return fLookupTable[index].second;
0084   }
0085
0086   private:
0087     LookupTable fLookupTable; ///< map ( Index, ( CL, leftside tail prob) )
0088
0089   protected:
0090     ClassDefOverride(SamplingSummaryLookup,1)  // A simple class used by ConfidenceBelt
0091   };
0092
0093
0094   ///////////////////////////
0095   class AcceptanceRegion : public TObject{
0096   public:
0097      AcceptanceRegion() : fLookupIndex(0), fLowerLimit(0), fUpperLimit(0) {}
0098
0099      AcceptanceRegion(Int_t lu, double ll, double ul) : fLookupIndex(lu), fLowerLimit(ll), fUpperLimit(ul) {}
0100      Int_t GetLookupIndex() { return fLookupIndex; }
0101      double GetLowerLimit() { return fLowerLimit; }
0102      double GetUpperLimit() { return fUpperLimit; }
0103
0104   private:
0105     Int_t fLookupIndex; // want a small footprint reference to the RooArgSet for particular parameter point
0106     double fLowerLimit;  // lower limit on test statistic
0107     double fUpperLimit;  // upper limit on test statistic
0108
0109   protected:
0110     ClassDefOverride(AcceptanceRegion,1)  // A simple class for acceptance regions used for ConfidenceBelt
0111
0112   };
0113
0114 ////////////////////////////////////////////////////////////////////////////////
0115
0116   class SamplingSummary : public TObject {
0117   public:
0118      SamplingSummary() : fParameterPointIndex(0) {}
0119      SamplingSummary(AcceptanceRegion& ar) : fParameterPointIndex(0) {
0120       AddAcceptanceRegion(ar);
0121     }
0122     Int_t GetParameterPointIndex(){return fParameterPointIndex;}
0123     SamplingDistribution* GetSamplingDistribution(){
0124       return (SamplingDistribution*) fSamplingDistribution.GetObject(); // dereference TRef
0125     }
0126     AcceptanceRegion& GetAcceptanceRegion(Int_t index=0){return fAcceptanceRegions[index];}
0127
0128     void AddAcceptanceRegion(AcceptanceRegion& ar){
0129       Int_t index =  ar.GetLookupIndex();
0130       if( fAcceptanceRegions.count(index) !=0) {
0131    std::cout << "SamplingSummary::AddAcceptanceRegion, need to implement merging protocol" << std::endl;
0132       } else {
0133    fAcceptanceRegions[index]=ar;
0134       }
0135     }
0136
0137   private:
0138      Int_t fParameterPointIndex; // want a small footprint reference to the RooArgSet for particular parameter point
0139      TRef fSamplingDistribution; // persistent pointer to a SamplingDistribution
0140      std::map<Int_t, AcceptanceRegion> fAcceptanceRegions;
0141
0142   protected:
0143     ClassDefOverride(SamplingSummary,1)  // A summary of acceptance regions for confidence belt
0144
0145   };
0146
0147 ////////////////////////////////////////////////////////////////////////////////
0148
0149  class ConfidenceBelt : public TNamed {
0150
0151   private:
0152     SamplingSummaryLookup fSamplingSummaryLookup;
0153     std::vector<SamplingSummary> fSamplingSummaries; // composite of several AcceptanceRegions
0154     RooAbsData* fParameterPoints = nullptr;  // either a histogram (RooDataHist) or a tree (RooDataSet)
0155
0156
0157   public:
0158     // constructors,destructors
0159     ConfidenceBelt() = default;
0160     ConfidenceBelt(const char* name);
0161     ConfidenceBelt(const char* name, const char* title);
0162     ConfidenceBelt(const char* name, RooAbsData&);
0163     ConfidenceBelt(const char* name, const char* title, RooAbsData&);
0164
0165     /// add after creating a region
0166     void AddAcceptanceRegion(RooArgSet&, AcceptanceRegion region, double cl=-1., double leftside=-1.);
0167
0168     /// add without creating a region, more useful for clients
0169     void AddAcceptanceRegion(RooArgSet& point, Int_t dataSetIndex, double lower, double upper, double cl=-1., double leftside=-1.);
0170
0171     AcceptanceRegion* GetAcceptanceRegion(RooArgSet&, double cl=-1., double leftside=-1.);
0172     double GetAcceptanceRegionMin(RooArgSet&, double cl=-1., double leftside=-1.);
0173     double GetAcceptanceRegionMax(RooArgSet&, double cl=-1., double leftside=-1.);
0174     std::vector<double> ConfidenceLevels() const ;
0175
0176     /// do we want it to return list of parameters
0177     virtual RooArgSet* GetParameters() const;
0178
0179     /// check if parameters are correct. (dummy implementation to start)
0180     bool CheckParameters(RooArgSet&) const ;
0181
0182   protected:
0183     ClassDefOverride(ConfidenceBelt,1)  // A confidence belt for the Neyman Construction
0184
0185   };
0186 }
0187
0188 #endif