ROOT logo
// @(#)root/roostats:$Id: ConfidenceBelt.h 31276 2009-11-18 15:06:42Z moneta $
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef RooStats_ConfidenceBelt
#define RooStats_ConfidenceBelt

#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif
#ifndef ROO_TREE_DATA
#include "RooAbsData.h"
#endif
#ifndef RooStats_ConfInterval
#include "RooStats/ConfInterval.h"
#endif

#include "RooStats/SamplingDistribution.h"

#include "TRef.h"

#include <vector>
#include <map>

using namespace std;

namespace RooStats {

  ///////////////////////////
  class SamplingSummaryLookup : public TObject {

    typedef pair<Double_t, Double_t> AcceptanceCriteria; // defined by Confidence level, leftside tail probability
    typedef map<Int_t, AcceptanceCriteria> LookupTable; // map ( Index, ( CL, leftside tail prob) )

  public:
    SamplingSummaryLookup() {}
    virtual ~SamplingSummaryLookup() {}

    void Add(Double_t cl, Double_t leftside){
      // add cl,leftside pair to lookup table
      AcceptanceCriteria tmp(cl, leftside);

      // should check to see if this is already in the map
      if(GetLookupIndex(cl,leftside) >=0 ){
	cout<< "SamplingSummaryLookup::Add, already in lookup table" << endl;
      } else
	fLookupTable[fLookupTable.size()]= tmp;
    }

    Int_t GetLookupIndex(Double_t cl, Double_t leftside){
      // get index for cl,leftside pair
      AcceptanceCriteria tmp(cl, leftside);

      Double_t tolerance = 1E-6; // some small number to protect floating point comparison.  What is better way?
      LookupTable::iterator it = fLookupTable.begin();
      Int_t index = -1;
      for(; it!= fLookupTable.end(); ++it) {
	index++;
	if( TMath::Abs( (*it).second.first - cl ) < tolerance &&
	    TMath::Abs( (*it).second.second - leftside ) < tolerance )
	  break; // exit loop, found 
      }

      // check that it was found
      if(index == (Int_t)fLookupTable.size())
	index = -1;

      return index;
    }

  Double_t GetConfidenceLevel(Int_t index){
    if(index<0 || index>(Int_t)fLookupTable.size()) {
      cout << "SamplingSummaryLookup::GetConfidenceLevel, index not in lookup table" << endl;
      return -1;
    }
    return fLookupTable[index].first;
  }

  Double_t GetLeftSideTailFraction(Int_t index){
    if(index<0 || index>(Int_t)fLookupTable.size()) {
      cout << "SamplingSummaryLookup::GetLeftSideTailFraction, index not in lookup table" << endl;
      return -1;
    }
    return fLookupTable[index].second;
  }

  private:
    LookupTable fLookupTable; // map ( Index, ( CL, leftside tail prob) )

  protected:
    ClassDef(SamplingSummaryLookup,1)  // A simple class used by ConfidenceBelt
  };


  ///////////////////////////
  class AcceptanceRegion : public TObject{
  public:
    AcceptanceRegion() {}
    virtual ~AcceptanceRegion() {}

    AcceptanceRegion(Int_t lu, Double_t ll, Double_t ul){
      fLookupIndex = lu;
      fLowerLimit = ll;
      fUpperLimit = ul;
    }
    Int_t GetLookupIndex(){return fLookupIndex;}
    Double_t GetLowerLimit(){return fLowerLimit;}
    Double_t GetUpperLimit(){return fUpperLimit;}

  private:
    Int_t fLookupIndex; // want a small footprint reference to the RooArgSet for particular parameter point
    Double_t fLowerLimit;  // lower limit on test statistic
    Double_t fUpperLimit;  // upper limit on test statistic

  protected:
    ClassDef(AcceptanceRegion,1)  // A simple class for acceptance regions used for ConfidenceBelt

  };


  ///////////////////////////
  class SamplingSummary : public TObject {
  public:
    SamplingSummary() {}
    virtual ~SamplingSummary() {}
    SamplingSummary(AcceptanceRegion& ar){
      AddAcceptanceRegion(ar);
    }
    Int_t GetParameterPointIndex(){return fParameterPointIndex;}
    SamplingDistribution* GetSamplingDistribution(){
      return (SamplingDistribution*) fSamplingDistribution.GetObject(); // dereference TRef
    }
    AcceptanceRegion& GetAcceptanceRegion(Int_t index=0){return fAcceptanceRegions[index];}

    void AddAcceptanceRegion(AcceptanceRegion& ar){
      Int_t index =  ar.GetLookupIndex();
      if( fAcceptanceRegions.count(index) !=0) {
	std::cout << "SamplingSummary::AddAcceptanceRegion, need to implement merging protocol" << std::endl;
      } else {
	fAcceptanceRegions[index]=ar;
      }
    }
    
  private:
    Int_t fParameterPointIndex; // want a small footprint reference to the RooArgSet for particular parameter point
    TRef fSamplingDistribution; // persistent pointer to a SamplingDistribution
    map<Int_t, AcceptanceRegion> fAcceptanceRegions;

  protected:
    ClassDef(SamplingSummary,1)  // A summary of acceptance regions for confidence belt

  };

  /////////////////////////////////////////
 class ConfidenceBelt : public TNamed {

  private:
    SamplingSummaryLookup fSamplingSummaryLookup;
    vector<SamplingSummary> fSamplingSummaries; // composite of several AcceptanceRegions
    RooAbsData* fParameterPoints;  // either a histogram (RooDataHist) or a tree (RooDataSet)


  public:
    // constructors,destructors
    ConfidenceBelt();
    ConfidenceBelt(const char* name);
    ConfidenceBelt(const char* name, const char* title);
    ConfidenceBelt(const char* name, RooAbsData&);
    ConfidenceBelt(const char* name, const char* title, RooAbsData&);
    virtual ~ConfidenceBelt();
        
    // add after creating a region
    void AddAcceptanceRegion(RooArgSet&, AcceptanceRegion region, Double_t cl=-1., Double_t leftside=-1.);

    // add without creating a region, more useful for clients
    void AddAcceptanceRegion(RooArgSet& point, Int_t dataSetIndex, Double_t lower, Double_t upper, Double_t cl=-1., Double_t leftside=-1.);

    AcceptanceRegion* GetAcceptanceRegion(RooArgSet&, Double_t cl=-1., Double_t leftside=-1.);
    Double_t GetAcceptanceRegionMin(RooArgSet&, Double_t cl=-1., Double_t leftside=-1.);
    Double_t GetAcceptanceRegionMax(RooArgSet&, Double_t cl=-1., Double_t leftside=-1.);
    vector<Double_t> ConfidenceLevels() const ;
 
    // Method to return lower limit on a given parameter 
    //  Double_t LowerLimit(RooRealVar& param) ; // could provide, but misleading?
    //      Double_t UpperLimit(RooRealVar& param) ; // could provide, but misleading?
    
    // do we want it to return list of parameters
    virtual RooArgSet* GetParameters() const;

    // check if parameters are correct. (dummy implementation to start)
    Bool_t CheckParameters(RooArgSet&) const ;
    
  protected:
    ClassDef(ConfidenceBelt,1)  // A confidence belt for the Neyman Construction
      
  };
}

#endif
 ConfidenceBelt.h:1
 ConfidenceBelt.h:2
 ConfidenceBelt.h:3
 ConfidenceBelt.h:4
 ConfidenceBelt.h:5
 ConfidenceBelt.h:6
 ConfidenceBelt.h:7
 ConfidenceBelt.h:8
 ConfidenceBelt.h:9
 ConfidenceBelt.h:10
 ConfidenceBelt.h:11
 ConfidenceBelt.h:12
 ConfidenceBelt.h:13
 ConfidenceBelt.h:14
 ConfidenceBelt.h:15
 ConfidenceBelt.h:16
 ConfidenceBelt.h:17
 ConfidenceBelt.h:18
 ConfidenceBelt.h:19
 ConfidenceBelt.h:20
 ConfidenceBelt.h:21
 ConfidenceBelt.h:22
 ConfidenceBelt.h:23
 ConfidenceBelt.h:24
 ConfidenceBelt.h:25
 ConfidenceBelt.h:26
 ConfidenceBelt.h:27
 ConfidenceBelt.h:28
 ConfidenceBelt.h:29
 ConfidenceBelt.h:30
 ConfidenceBelt.h:31
 ConfidenceBelt.h:32
 ConfidenceBelt.h:33
 ConfidenceBelt.h:34
 ConfidenceBelt.h:35
 ConfidenceBelt.h:36
 ConfidenceBelt.h:37
 ConfidenceBelt.h:38
 ConfidenceBelt.h:39
 ConfidenceBelt.h:40
 ConfidenceBelt.h:41
 ConfidenceBelt.h:42
 ConfidenceBelt.h:43
 ConfidenceBelt.h:44
 ConfidenceBelt.h:45
 ConfidenceBelt.h:46
 ConfidenceBelt.h:47
 ConfidenceBelt.h:48
 ConfidenceBelt.h:49
 ConfidenceBelt.h:50
 ConfidenceBelt.h:51
 ConfidenceBelt.h:52
 ConfidenceBelt.h:53
 ConfidenceBelt.h:54
 ConfidenceBelt.h:55
 ConfidenceBelt.h:56
 ConfidenceBelt.h:57
 ConfidenceBelt.h:58
 ConfidenceBelt.h:59
 ConfidenceBelt.h:60
 ConfidenceBelt.h:61
 ConfidenceBelt.h:62
 ConfidenceBelt.h:63
 ConfidenceBelt.h:64
 ConfidenceBelt.h:65
 ConfidenceBelt.h:66
 ConfidenceBelt.h:67
 ConfidenceBelt.h:68
 ConfidenceBelt.h:69
 ConfidenceBelt.h:70
 ConfidenceBelt.h:71
 ConfidenceBelt.h:72
 ConfidenceBelt.h:73
 ConfidenceBelt.h:74
 ConfidenceBelt.h:75
 ConfidenceBelt.h:76
 ConfidenceBelt.h:77
 ConfidenceBelt.h:78
 ConfidenceBelt.h:79
 ConfidenceBelt.h:80
 ConfidenceBelt.h:81
 ConfidenceBelt.h:82
 ConfidenceBelt.h:83
 ConfidenceBelt.h:84
 ConfidenceBelt.h:85
 ConfidenceBelt.h:86
 ConfidenceBelt.h:87
 ConfidenceBelt.h:88
 ConfidenceBelt.h:89
 ConfidenceBelt.h:90
 ConfidenceBelt.h:91
 ConfidenceBelt.h:92
 ConfidenceBelt.h:93
 ConfidenceBelt.h:94
 ConfidenceBelt.h:95
 ConfidenceBelt.h:96
 ConfidenceBelt.h:97
 ConfidenceBelt.h:98
 ConfidenceBelt.h:99
 ConfidenceBelt.h:100
 ConfidenceBelt.h:101
 ConfidenceBelt.h:102
 ConfidenceBelt.h:103
 ConfidenceBelt.h:104
 ConfidenceBelt.h:105
 ConfidenceBelt.h:106
 ConfidenceBelt.h:107
 ConfidenceBelt.h:108
 ConfidenceBelt.h:109
 ConfidenceBelt.h:110
 ConfidenceBelt.h:111
 ConfidenceBelt.h:112
 ConfidenceBelt.h:113
 ConfidenceBelt.h:114
 ConfidenceBelt.h:115
 ConfidenceBelt.h:116
 ConfidenceBelt.h:117
 ConfidenceBelt.h:118
 ConfidenceBelt.h:119
 ConfidenceBelt.h:120
 ConfidenceBelt.h:121
 ConfidenceBelt.h:122
 ConfidenceBelt.h:123
 ConfidenceBelt.h:124
 ConfidenceBelt.h:125
 ConfidenceBelt.h:126
 ConfidenceBelt.h:127
 ConfidenceBelt.h:128
 ConfidenceBelt.h:129
 ConfidenceBelt.h:130
 ConfidenceBelt.h:131
 ConfidenceBelt.h:132
 ConfidenceBelt.h:133
 ConfidenceBelt.h:134
 ConfidenceBelt.h:135
 ConfidenceBelt.h:136
 ConfidenceBelt.h:137
 ConfidenceBelt.h:138
 ConfidenceBelt.h:139
 ConfidenceBelt.h:140
 ConfidenceBelt.h:141
 ConfidenceBelt.h:142
 ConfidenceBelt.h:143
 ConfidenceBelt.h:144
 ConfidenceBelt.h:145
 ConfidenceBelt.h:146
 ConfidenceBelt.h:147
 ConfidenceBelt.h:148
 ConfidenceBelt.h:149
 ConfidenceBelt.h:150
 ConfidenceBelt.h:151
 ConfidenceBelt.h:152
 ConfidenceBelt.h:153
 ConfidenceBelt.h:154
 ConfidenceBelt.h:155
 ConfidenceBelt.h:156
 ConfidenceBelt.h:157
 ConfidenceBelt.h:158
 ConfidenceBelt.h:159
 ConfidenceBelt.h:160
 ConfidenceBelt.h:161
 ConfidenceBelt.h:162
 ConfidenceBelt.h:163
 ConfidenceBelt.h:164
 ConfidenceBelt.h:165
 ConfidenceBelt.h:166
 ConfidenceBelt.h:167
 ConfidenceBelt.h:168
 ConfidenceBelt.h:169
 ConfidenceBelt.h:170
 ConfidenceBelt.h:171
 ConfidenceBelt.h:172
 ConfidenceBelt.h:173
 ConfidenceBelt.h:174
 ConfidenceBelt.h:175
 ConfidenceBelt.h:176
 ConfidenceBelt.h:177
 ConfidenceBelt.h:178
 ConfidenceBelt.h:179
 ConfidenceBelt.h:180
 ConfidenceBelt.h:181
 ConfidenceBelt.h:182
 ConfidenceBelt.h:183
 ConfidenceBelt.h:184
 ConfidenceBelt.h:185
 ConfidenceBelt.h:186
 ConfidenceBelt.h:187
 ConfidenceBelt.h:188
 ConfidenceBelt.h:189
 ConfidenceBelt.h:190
 ConfidenceBelt.h:191
 ConfidenceBelt.h:192
 ConfidenceBelt.h:193
 ConfidenceBelt.h:194
 ConfidenceBelt.h:195
 ConfidenceBelt.h:196
 ConfidenceBelt.h:197
 ConfidenceBelt.h:198
 ConfidenceBelt.h:199
 ConfidenceBelt.h:200
 ConfidenceBelt.h:201
 ConfidenceBelt.h:202
 ConfidenceBelt.h:203
 ConfidenceBelt.h:204
 ConfidenceBelt.h:205