/usr/include/root/RooStats/MCMCInterval.h is in libroot-roofit-dev 5.34.19+dfsg-1.2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 | // @(#)root/roostats:$Id$
// Authors: Kevin Belasco 17/06/2009
// Authors: Kyle Cranmer 17/06/2009
/*************************************************************************
* 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_MCMCInterval
#define RooStats_MCMCInterval
#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif
#ifndef ROOSTATS_ConfInterval
#include "RooStats/ConfInterval.h"
#endif
#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif
#ifndef ROO_ARG_LIST
#include "RooArgList.h"
#endif
#ifndef ROOSTATS_MarkovChain
#include "RooStats/MarkovChain.h"
#endif
class RooNDKeysPdf;
class RooProduct;
namespace RooStats {
class Heaviside;
class MCMCInterval : public ConfInterval {
public:
// default constructor
explicit MCMCInterval(const char* name = 0);
// constructor from parameter of interest and Markov chain object
MCMCInterval(const char* name, const RooArgSet& parameters,
MarkovChain& chain);
enum {DEFAULT_NUM_BINS = 50};
enum IntervalType {kShortest, kTailFraction};
virtual ~MCMCInterval();
// determine whether this point is in the confidence interval
virtual Bool_t IsInInterval(const RooArgSet& point) const;
// set the desired confidence level (see GetActualConfidenceLevel())
// Note: calling this function triggers the algorithm that determines
// the interval, so call this after initializing all other aspects
// of this IntervalCalculator
// Also, calling this function again with a different confidence level
// retriggers the calculation of the interval
virtual void SetConfidenceLevel(Double_t cl);
// get the desired confidence level (see GetActualConfidenceLevel())
virtual Double_t ConfidenceLevel() const {return fConfidenceLevel;}
// return a set containing the parameters of this interval
// the caller owns the returned RooArgSet*
virtual RooArgSet* GetParameters() const;
// get the cutoff bin height for being considered in the
// confidence interval
virtual Double_t GetHistCutoff();
// get the cutoff RooNDKeysPdf value for being considered in the
// confidence interval
virtual Double_t GetKeysPdfCutoff();
//virtual Double_t GetKeysPdfCutoff() { return fKeysCutoff; }
// get the actual value of the confidence level for this interval.
virtual Double_t GetActualConfidenceLevel();
// whether the specified confidence level is a floor for the actual
// confidence level (strict), or a ceiling (not strict)
virtual void SetHistStrict(Bool_t isHistStrict)
{ fIsHistStrict = isHistStrict; }
// check if parameters are correct. (dummy implementation to start)
Bool_t CheckParameters(const RooArgSet& point) const;
// Set the parameters of interest for this interval
// and change other internal data members accordingly
virtual void SetParameters(const RooArgSet& parameters);
// Set the MarkovChain that this interval is based on
virtual void SetChain(MarkovChain& chain) { fChain = &chain; }
// Set which parameters go on which axis. The first list element
// goes on the x axis, second (if it exists) on y, third (if it
// exists) on z, etc
virtual void SetAxes(RooArgList& axes);
// return a list of RooRealVars representing the axes
// you own the returned RooArgList
virtual RooArgList* GetAxes()
{
RooArgList* axes = new RooArgList();
for (Int_t i = 0; i < fDimension; i++)
axes->addClone(*fAxes[i]);
return axes;
}
// get the lowest value of param that is within the confidence interval
virtual Double_t LowerLimit(RooRealVar& param);
// determine lower limit of the lower confidence interval
virtual Double_t LowerLimitTailFraction(RooRealVar& param);
// get the lower limit of param in the shortest confidence interval
// Note that this works better for some distributions (ones with exactly
// one maximum) than others, and sometimes has little value.
virtual Double_t LowerLimitShortest(RooRealVar& param);
// determine lower limit in the shortest interval by using keys pdf
virtual Double_t LowerLimitByKeys(RooRealVar& param);
// determine lower limit using histogram
virtual Double_t LowerLimitByHist(RooRealVar& param);
// determine lower limit using histogram
virtual Double_t LowerLimitBySparseHist(RooRealVar& param);
// determine lower limit using histogram
virtual Double_t LowerLimitByDataHist(RooRealVar& param);
// get the highest value of param that is within the confidence interval
virtual Double_t UpperLimit(RooRealVar& param);
// determine upper limit of the lower confidence interval
virtual Double_t UpperLimitTailFraction(RooRealVar& param);
// get the upper limit of param in the confidence interval
// Note that this works better for some distributions (ones with exactly
// one maximum) than others, and sometimes has little value.
virtual Double_t UpperLimitShortest(RooRealVar& param);
// determine upper limit in the shortest interval by using keys pdf
virtual Double_t UpperLimitByKeys(RooRealVar& param);
// determine upper limit using histogram
virtual Double_t UpperLimitByHist(RooRealVar& param);
// determine upper limit using histogram
virtual Double_t UpperLimitBySparseHist(RooRealVar& param);
// determine upper limit using histogram
virtual Double_t UpperLimitByDataHist(RooRealVar& param);
// Determine the approximate maximum value of the Keys PDF
Double_t GetKeysMax();
// set the number of steps in the chain to discard as burn-in,
// starting from the first
virtual void SetNumBurnInSteps(Int_t numBurnInSteps)
{ fNumBurnInSteps = numBurnInSteps; }
// set whether to use kernel estimation to determine the interval
virtual void SetUseKeys(Bool_t useKeys) { fUseKeys = useKeys; }
// set whether to use a sparse histogram. you MUST also call
// SetUseKeys(kFALSE) to use a histogram.
virtual void SetUseSparseHist(Bool_t useSparseHist)
{ fUseSparseHist = useSparseHist; }
// get whether we used kernel estimation to determine the interval
virtual Bool_t GetUseKeys() { return fUseKeys; }
// get the number of steps in the chain to disard as burn-in,
// get the number of steps in the chain to disard as burn-in,
// starting from the first
virtual Int_t GetNumBurnInSteps() { return fNumBurnInSteps; }
// set the number of bins to use (same for all axes, for now)
//virtual void SetNumBins(Int_t numBins);
// Get a clone of the histogram of the posterior
virtual TH1* GetPosteriorHist();
// Get a clone of the keys pdf of the posterior
virtual RooNDKeysPdf* GetPosteriorKeysPdf();
// Get a clone of the (keyspdf * heaviside) product of the posterior
virtual RooProduct* GetPosteriorKeysProduct();
// Get the number of parameters of interest in this interval
virtual Int_t GetDimension() const { return fDimension; }
// Get the markov chain on which this interval is based
// You do not own the returned MarkovChain*
virtual const MarkovChain* GetChain() { return fChain; }
// Get a clone of the markov chain on which this interval is based
// as a RooDataSet. You own the returned RooDataSet*
virtual RooDataSet* GetChainAsDataSet(RooArgSet* whichVars = NULL)
{ return fChain->GetAsDataSet(whichVars); }
// Get the markov chain on which this interval is based
// as a RooDataSet. You do not own the returned RooDataSet*
virtual const RooDataSet* GetChainAsConstDataSet()
{ return fChain->GetAsConstDataSet(); }
// Get a clone of the markov chain on which this interval is based
// as a RooDataHist. You own the returned RooDataHist*
virtual RooDataHist* GetChainAsDataHist(RooArgSet* whichVars = NULL)
{ return fChain->GetAsDataHist(whichVars); }
// Get a clone of the markov chain on which this interval is based
// as a THnSparse. You own the returned THnSparse*
virtual THnSparse* GetChainAsSparseHist(RooArgSet* whichVars = NULL)
{ return fChain->GetAsSparseHist(whichVars); }
// Get a clone of the NLL variable from the markov chain
virtual RooRealVar* GetNLLVar() const
{ return fChain->GetNLLVar(); }
// Get a clone of the weight variable from the markov chain
virtual RooRealVar* GetWeightVar() const
{ return fChain->GetWeightVar(); }
// set the acceptable level or error for Keys interval determination
virtual void SetEpsilon(Double_t epsilon)
{
if (epsilon < 0)
coutE(InputArguments) << "MCMCInterval::SetEpsilon will not allow "
<< "negative epsilon value" << std::endl;
else
fEpsilon = epsilon;
}
// Set the type of interval to find. This will only have an effect for
// 1-D intervals. If is more than 1 parameter of interest, then a
// "shortest" interval will always be used, since it generalizes directly
// to N dimensions
virtual void SetIntervalType(enum IntervalType intervalType)
{ fIntervalType = intervalType; }
// Return the type of this interval
virtual enum IntervalType GetIntervalType() { return fIntervalType; }
// set the left-side tail fraction for a tail-fraction interval
virtual void SetLeftSideTailFraction(Double_t a) { fLeftSideTF = a; }
// kbelasco: The inner-workings of the class really should not be exposed
// like this in a comment, but it seems to be the only way to give
// the user any control over this process, if he desires it
//
// Set the fraction delta such that
// topCutoff (a) is considered == bottomCutoff (b) iff
// (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2))
// when determining the confidence interval by Keys
virtual void SetDelta(Double_t delta)
{
if (delta < 0.)
coutE(InputArguments) << "MCMCInterval::SetDelta will not allow "
<< "negative delta value" << std::endl;
else
fDelta = delta;
}
private:
inline Bool_t AcceptableConfLevel(Double_t confLevel);
inline Bool_t WithinDeltaFraction(Double_t a, Double_t b);
protected:
// data members
RooArgSet fParameters; // parameters of interest for this interval
MarkovChain* fChain; // the markov chain
Double_t fConfidenceLevel; // Requested confidence level (eg. 0.95 for 95% CL)
RooDataHist* fDataHist; // the binned Markov Chain data
THnSparse* fSparseHist; // the binned Markov Chain data
Double_t fHistConfLevel; // the actual conf level determined by hist
Double_t fHistCutoff; // cutoff bin size to be in interval
RooNDKeysPdf* fKeysPdf; // the kernel estimation pdf
RooProduct* fProduct; // the (keysPdf * heaviside) product
Heaviside* fHeaviside; // the Heaviside function
RooDataHist* fKeysDataHist; // data hist representing product
RooRealVar* fCutoffVar; // cutoff variable to use for integrating keys pdf
Double_t fKeysConfLevel; // the actual conf level determined by keys
Double_t fKeysCutoff; // cutoff keys pdf value to be in interval
Double_t fFull; // Value of intergral of fProduct
Double_t fLeftSideTF; // left side tail-fraction for interval
Double_t fTFConfLevel; // the actual conf level of tail-fraction interval
std::vector<Int_t> fVector; // vector containing the Markov chain data
Double_t fVecWeight; // sum of weights of all entries in fVector
Double_t fTFLower; // lower limit of the tail-fraction interval
Double_t fTFUpper; // upper limit of the tail-fraction interval
TH1* fHist; // the binned Markov Chain data
Bool_t fUseKeys; // whether to use kernel estimation
Bool_t fUseSparseHist; // whether to use sparse hist (vs. RooDataHist)
Bool_t fIsHistStrict; // whether the specified confidence level is a
// floor for the actual confidence level (strict),
// or a ceiling (not strict) for determination by
// histogram
Int_t fDimension; // number of variables
Int_t fNumBurnInSteps; // number of steps to discard as burn in, starting
// from the first
// LM (not used) Double_t fIntervalSum; // sum of heights of bins in the interval
RooRealVar** fAxes; // array of pointers to RooRealVars representing
// the axes of the histogram
// fAxes[0] represents x-axis, [1] y, [2] z, etc
Double_t fEpsilon; // acceptable error for Keys interval determination
Double_t fDelta; // topCutoff (a) considered == bottomCutoff (b) iff
// (TMath::Abs(a - b) < TMath::Abs(fDelta * (a + b)/2));
// Theoretically, the Abs is not needed here, but
// floating-point arithmetic does not always work
// perfectly, and the Abs doesn't hurt
enum IntervalType fIntervalType;
// functions
virtual void DetermineInterval();
virtual void DetermineShortestInterval();
virtual void DetermineTailFractionInterval();
virtual void DetermineByHist();
virtual void DetermineBySparseHist();
virtual void DetermineByDataHist();
virtual void DetermineByKeys();
virtual void CreateHist();
virtual void CreateSparseHist();
virtual void CreateDataHist();
virtual void CreateKeysPdf();
virtual void CreateKeysDataHist();
virtual void CreateVector(RooRealVar* param);
inline virtual Double_t CalcConfLevel(Double_t cutoff, Double_t full);
ClassDef(MCMCInterval,1) // Concrete implementation of a ConfInterval based on MCMC calculation
};
}
#endif
|