96MCMCInterval::MCMCInterval(
const char*
name)
189 fDataHist->get(bin1);
190 double n1 = fDataHist->weight();
191 fDataHist->get(bin2);
192 double n2 = fDataHist->weight();
201 double n1 = fSparseHist->GetBinContent(bin1);
202 double n2 = fSparseHist->GetBinContent(bin2);
210 fChain(chain), fParam(param) {}
212 double xi = fChain->Get(i)->getRealValue(fParam->GetName());
213 double xj = fChain->Get(j)->getRealValue(fParam->GetName());
251 x[i] =
fAxes[i]->getVal();
280 <<
"Interval type not set. Returning false." << endl;
321 "number of variables in axes (" <<
size <<
322 ") doesn't match number of parameters ("
339 <<
"parameters have not been set." << endl;
345 "MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: " <<
346 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
347 "in Markov chain." << endl;
379 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist(): " <<
380 "Crucial data member was nullptr." << endl;
381 coutE(
Eval) <<
"Make sure to fully construct/initialize." << endl;
384 if (
fHist !=
nullptr) {
391 "MCMCInterval::CreateHist: creation of histogram failed: " <<
392 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
393 "in Markov chain." << endl;
398 fHist =
new TH1F(
"posterior",
"MCMC Posterior Histogram",
402 fHist =
new TH2F(
"posterior",
"MCMC Posterior Histogram",
407 fHist =
new TH3F(
"posterior",
"MCMC Posterior Histogram",
413 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist() : " <<
414 "TH1* couldn't handle dimension: " <<
fDimension << endl;
451 <<
"Crucial data member was nullptr." << endl;
468 fDimension, bins.data(), min.data(), max.data());
477 "MCMCInterval::CreateSparseHist: creation of histogram failed: " <<
478 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
479 "in Markov chain." << endl;
499 coutE(
Eval) <<
"* Error in MCMCInterval::CreateDataHist(): " <<
500 "Crucial data member was nullptr or empty." << endl;
501 coutE(
Eval) <<
"Make sure to fully construct/initialize." << endl;
507 "MCMCInterval::CreateDataHist: creation of histogram failed: " <<
508 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
509 "in Markov chain." << endl;
527 "Crucial data member (Markov chain) was nullptr." << endl;
535 "MCMCInterval::CreateVector: creation of vector failed: " <<
536 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
537 "in Markov chain." << endl;
545 for (i = 0; i <
size; i++) {
562 if (
fAxes !=
nullptr)
567 if (
dynamic_cast<RooRealVar*
>(obj) !=
nullptr)
570 coutE(
Eval) <<
"* Error in MCMCInterval::SetParameters: " <<
571 obj->GetName() <<
" not a RooRealVar*" << std::endl;
589 "Error: Interval type not set" << endl;
608 if (fLeftSideTF < 0 || fLeftSideTF > 1) {
610 <<
"Fraction must be in the range [0, 1]. "
617 <<
"Error: Can only find a tail-fraction interval for 1-D intervals"
622 if (
fAxes ==
nullptr) {
624 <<
"Crucial data member was nullptr." << endl;
657 double leftTailSum = 0;
658 double rightTailSum = 0;
661 double ll = param->
getMin();
662 double ul = param->
getMax();
691 if (
TMath::Abs(rightTailSum +
w - rightTailCutoff) <
732 coutW(
Eval) <<
"Warning: Integral of Keys PDF came out to " << full
733 <<
" instead of expected value 1. Will continue using this "
734 <<
"factor to normalize further integrals of this PDF." << endl;
743 for (
auto *var : static_range_cast<RooRealVar*>(
fParameters))
744 volume *= (var->getMax() - var->getMin());
746 double topCutoff = full / volume;
747 double bottomCutoff = topCutoff;
754 bool changed =
false;
767 bottomCutoff = topCutoff / 2.0;
788 topCutoff = bottomCutoff * 2.0;
792 coutI(
Eval) <<
"range set: [" << bottomCutoff <<
", " << topCutoff <<
"]"
795 cutoff = (topCutoff + bottomCutoff) / 2.0;
807 bottomCutoff = cutoff;
810 cutoff = (topCutoff + bottomCutoff) / 2.0;
811 coutI(
Eval) <<
"cutoff range: [" << bottomCutoff <<
", "
812 << topCutoff <<
"]" << endl;
847 std::vector<Long_t> bins(numBins);
848 for (
Int_t ibin = 0; ibin < numBins; ibin++)
849 bins[ibin] = (
Long_t)ibin;
857 for (i = numBins - 1; i >= 0; i--) {
875 for ( ; i >= 0; i--) {
884 for ( ; i < numBins; i++) {
891 if (i == numBins - 1)
917 std::vector<Int_t> bins(numBins);
918 for (
Int_t ibin = 0; ibin < numBins; ibin++)
926 for (i = numBins - 1; i >= 0; i--) {
945 for ( ; i >= 0; i--) {
955 for ( ; i < numBins; i++) {
963 if (i == numBins - 1)
986 <<
"not implemented for this type of interval. Returning 0." << endl;
1002 "Error: Interval type not set" << endl;
1018 "Error: Interval type not set" << endl;
1095 <<
"Sorry, will not compute lower limit unless dimension == 1" << endl;
1103 coutE(
Eval) <<
"In MCMCInterval::LowerLimitBySparseHist: "
1104 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1105 <<
"steps in the Markov chain. Returning param.getMin()." << endl;
1113 double lowerLimit = param.
getMax();
1115 for (
Long_t i = 0; i < numBins; i++) {
1118 if (val < lowerLimit)
1139 coutE(
Eval) <<
"In MCMCInterval::LowerLimitByDataHist: "
1140 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1141 <<
"steps in the Markov chain. Returning param.getMin()." << endl;
1148 double lowerLimit = param.
getMax();
1150 for (
Int_t i = 0; i < numBins; i++) {
1154 if (val < lowerLimit)
1172 <<
"Sorry, will not compute upper limit unless dimension == 1" << endl;
1180 coutE(
Eval) <<
"In MCMCInterval::UpperLimitBySparseHist: "
1181 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1182 <<
"steps in the Markov chain. Returning param.getMax()." << endl;
1190 double upperLimit = param.
getMin();
1192 for (
Long_t i = 0; i < numBins; i++) {
1195 if (val > upperLimit)
1216 coutE(
Eval) <<
"In MCMCInterval::UpperLimitByDataHist: "
1217 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1218 <<
"steps in the Markov chain. Returning param.getMax()." << endl;
1225 double upperLimit = param.
getMin();
1227 for (
Int_t i = 0; i < numBins; i++) {
1231 if (val > upperLimit)
1255 coutE(
Eval) <<
"in MCMCInterval::LowerLimitByKeys(): "
1256 <<
"couldn't find lower limit, check that the number of burn in "
1257 <<
"steps < number of total steps in the Markov chain. Returning "
1258 <<
"param.getMin()" << endl;
1265 double lowerLimit = param.
getMax();
1267 for (
Int_t i = 0; i < numBins; i++) {
1271 if (val < lowerLimit)
1295 coutE(
Eval) <<
"in MCMCInterval::UpperLimitByKeys(): "
1296 <<
"couldn't find upper limit, check that the number of burn in "
1297 <<
"steps < number of total steps in the Markov chain. Returning "
1298 <<
"param.getMax()" << endl;
1305 double upperLimit = param.
getMin();
1307 for (
Int_t i = 0; i < numBins; i++) {
1311 if (val > upperLimit)
1334 coutE(
Eval) <<
"in MCMCInterval::KeysMax(): "
1335 <<
"couldn't find Keys max value, check that the number of burn in "
1336 <<
"steps < number of total steps in the Markov chain. Returning 0"
1344 for (
Int_t i = 0; i < numBins; i++) {
1387 coutI(
Eval) <<
"cutoff = " << cutoff <<
", conf = " << confLevel << endl;
1399 <<
"confidence level not set " << endl;
1400 if (
fHist ==
nullptr)
1403 if (
fHist ==
nullptr)
1416 <<
"confidence level not set " << endl;
1433 <<
"confidence level not set " << endl;
1472 if (
fAxes ==
nullptr)
1498 bool tempChangeBinning =
true;
1500 if (!
fAxes[i]->getBinning(
nullptr,
false,
false).isUniform()) {
1501 tempChangeBinning =
false;
1510 tempChangeBinning =
false;
1512 if (tempChangeBinning) {
1524 "Keys PDF & Heaviside Product Data Hist",
fParameters);
1527 if (tempChangeBinning) {
1532 fAxes[i]->setBins(savedBins[i],
nullptr);
1543 coutE(
Eval) <<
"MCMCInterval: size is wrong, parameters don't match" << std::endl;
1547 coutE(
Eval) <<
"MCMCInterval: size is ok, but parameters don't match" << std::endl;
static const double DEFAULT_EPSILON
static const double DEFAULT_DELTA
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
THnSparseT< TArrayF > THnSparseF
TRObject operator()(const T1 &t1) const
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Int_t numBins() const
Return number of bins.
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Int_t numBins(const char *rangeName=nullptr) const override
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, double scaleFactor, bool correctForBinVolume=false, bool showProgress=false) const
Fill a RooDataHist with values sampled from this function at the bin centers.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables std::liste...
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
The RooDataHist is a container class to hold N-dimensional binned data.
double sum(bool correctForBinSize, bool inverseCorr=false) const
Return the sum of the weights of all bins in the histogram.
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
double weight(std::size_t i) const
Return weight of i-th bin.
const RooArgSet * get() const override
Get bin centre of current bin.
RooDataSet is a container class to hold unbinned data.
Generic N-dimensional implementation of a kernel estimation p.d.f.
static double infinity()
Return internal infinity representation.
A RooProduct represents the product of a given set of RooAbsReal objects.
RooRealVar represents a variable that can be changed from the outside.
void setVal(double value) override
Set value of variable to 'value'.
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
ConfInterval is an interface class for a generic interval in the RooStats framework.
Represents the Heaviside function.
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual void CreateDataHist()
virtual void DetermineByDataHist()
virtual void DetermineShortestInterval()
virtual double GetActualConfidenceLevel()
virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
double fKeysConfLevel
the actual conf level determined by keys
virtual void CreateVector(RooRealVar *param)
RooDataHist * fDataHist
the binned Markov Chain data
virtual double LowerLimitByDataHist(RooRealVar ¶m)
determine lower limit using histogram
double fDelta
topCutoff (a) considered == bottomCutoff (b) iff
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
TH1 * fHist
the binned Markov Chain data
virtual double UpperLimitBySparseHist(RooRealVar ¶m)
determine upper limit using histogram
enum IntervalType fIntervalType
virtual double UpperLimit(RooRealVar ¶m)
get the highest value of param that is within the confidence interval
RooRealVar * fCutoffVar
cutoff variable to use for integrating keys pdf
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void DetermineByKeys()
double fTFLower
lower limit of the tail-fraction interval
bool fUseSparseHist
whether to use sparse hist (vs. RooDataHist)
double fVecWeight
sum of weights of all entries in fVector
double fLeftSideTF
left side tail-fraction for interval
bool AcceptableConfLevel(double confLevel)
virtual void DetermineBySparseHist()
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
RooProduct * fProduct
the (keysPdf * heaviside) product
virtual double LowerLimitBySparseHist(RooRealVar ¶m)
determine lower limit using histogram
bool WithinDeltaFraction(double a, double b)
double fKeysCutoff
cutoff keys pdf value to be in interval
virtual double LowerLimitShortest(RooRealVar ¶m)
get the lower limit of param in the shortest confidence interval Note that this works better for some...
virtual double LowerLimitTailFraction(RooRealVar ¶m)
determine lower limit of the lower confidence interval
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
Int_t fDimension
number of variables
void SetConfidenceLevel(double cl) override
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
MCMCInterval(const char *name=nullptr)
default constructor
RooRealVar ** fAxes
array of pointers to RooRealVars representing the axes of the histogram fAxes[0] represents x-axis,...
virtual void DetermineTailFractionInterval()
double fTFConfLevel
the actual conf level of tail-fraction interval
double fHistConfLevel
the actual conf level determined by hist
virtual void CreateSparseHist()
Heaviside * fHeaviside
the Heaviside function
double fFull
Value of intergral of fProduct.
virtual double UpperLimitByDataHist(RooRealVar ¶m)
determine upper limit using histogram
virtual void DetermineInterval()
virtual double CalcConfLevel(double cutoff, double full)
virtual double LowerLimitByHist(RooRealVar ¶m)
determine lower limit using histogram
virtual double LowerLimit(RooRealVar ¶m)
get the lowest value of param that is within the confidence interval
virtual double UpperLimitShortest(RooRealVar ¶m)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
THnSparse * fSparseHist
the binned Markov Chain data
bool CheckParameters(const RooArgSet &point) const override
check if parameters are correct. (dummy implementation to start)
std::vector< Int_t > fVector
vector containing the Markov chain data
virtual double LowerLimitByKeys(RooRealVar ¶m)
determine lower limit in the shortest interval by using keys pdf
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
double fTFUpper
upper limit of the tail-fraction interval
double fHistCutoff
cutoff bin size to be in interval
MarkovChain * fChain
the markov chain
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
virtual double UpperLimitTailFraction(RooRealVar ¶m)
determine upper limit of the lower confidence interval
bool fIsHistStrict
whether the specified confidence level is a
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual double UpperLimitByHist(RooRealVar ¶m)
determine upper limit using histogram
virtual void CreateKeysPdf()
virtual double UpperLimitByKeys(RooRealVar ¶m)
determine upper limit in the shortest interval by using keys pdf
Int_t fNumBurnInSteps
number of steps to discard as burn in, starting from the first
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
RooDataHist * fKeysDataHist
data hist representing product
RooArgSet fParameters
parameters of interest for this interval
bool IsInInterval(const RooArgSet &point) const override
determine whether this point is in the confidence interval
virtual void DetermineByHist()
virtual void SetParameters(const RooArgSet ¶meters)
Set the parameters of interest for this interval and change other internal data members accordingly.
virtual void CreateKeysDataHist()
RooNDKeysPdf * fKeysPdf
the kernel estimation pdf
double fEpsilon
acceptable error for Keys interval determination
virtual double GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual double GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Stores the steps in a Markov Chain of points.
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
virtual double Weight() const
get the weight of the current (last indexed) entry
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
virtual Int_t Size() const
get the number of steps in the chain
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
1-D histogram with a float per channel (see TH1 documentation)}
TH1 is the base class of all histogram classes in ROOT.
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
2-D histogram with a float per channel (see TH1 documentation)}
3-D histogram with a float per channel (see TH1 documentation)}
Long64_t Fill(const Double_t *x, Double_t w=1.)
TAxis * GetAxis(Int_t dim) const
Efficient multidimensional histogram.
Double_t GetBinContent(const Int_t *idx) const
Forwards to THnBase::GetBinContent() overload.
Long64_t GetBin(const Int_t *idx) const override
void Sumw2() override
Enable calculation of errors.
Long64_t GetNbins() const override
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
const char * GetName() const override
Returns name of object.
RooCmdArg SelectVars(const RooArgSet &vars)
RooCmdArg EventRange(Int_t nStart, Int_t nStop)
RooCmdArg NormSet(Args_t &&... argsOrArgSet)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
CompareDataHistBins(RooDataHist *hist)
CompareSparseHistBins(THnSparse *hist)
CompareVectorIndices(MarkovChain *chain, RooRealVar *param)