Logo ROOT  
Reference Guide
MarkovChain.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Authors: Kevin Belasco 17/06/2009
3// Authors: Kyle Cranmer 17/06/2009
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class RooStats::MarkovChain
13 \ingroup Roostats
14
15Stores the steps in a Markov Chain of points. Allows user to access the
16weight and NLL value (if applicable) with which a point was added to the
17MarkovChain.
18
19*/
20
21#include "Rtypes.h"
22
23#include "TNamed.h"
25#include "RooDataSet.h"
26#include "RooArgSet.h"
27#include "RooRealVar.h"
29#include "RooDataHist.h"
30#include "THnSparse.h"
31
32using namespace std;
33
35
36using namespace RooFit;
37using namespace RooStats;
38
39static const char* WEIGHT_NAME = "weight_MarkovChain_local_";
40static const char* NLL_NAME = "nll_MarkovChain_local_";
41static const char* DATASET_NAME = "dataset_MarkovChain_local_";
42static const char* DEFAULT_NAME = "_markov_chain";
43static const char* DEFAULT_TITLE = "Markov Chain";
44
45MarkovChain::MarkovChain() :
47{
48 fParameters = nullptr;
49 fDataEntry = nullptr;
50 fChain = nullptr;
51 fNLL = nullptr;
52 fWeight = nullptr;
53}
54
57{
58 fParameters = nullptr;
59 fDataEntry = nullptr;
60 fChain = nullptr;
61 fNLL = nullptr;
62 fWeight = nullptr;
63 SetParameters(parameters);
64}
65
66MarkovChain::MarkovChain(const char* name, const char* title,
67 RooArgSet& parameters) : TNamed(name, title)
68{
69 fParameters = nullptr;
70 fDataEntry = nullptr;
71 fChain = nullptr;
72 fNLL = nullptr;
73 fWeight = nullptr;
74 SetParameters(parameters);
75}
76
78{
79 delete fChain;
80 delete fParameters;
81 delete fDataEntry;
82
83 fParameters = new RooArgSet();
84 fParameters->addClone(parameters);
85
86 // kbelasco: consider setting fDataEntry = fChain->get()
87 // to see if that makes it possible to get values of variables without
88 // doing string comparison
89 RooRealVar nll(NLL_NAME, "-log Likelihood", 0);
90 RooRealVar weight(WEIGHT_NAME, "weight", 0);
91
92 fDataEntry = new RooArgSet();
93 fDataEntry->addClone(parameters);
94 fDataEntry->addClone(nll);
95 fDataEntry->addClone(weight);
98
99 fChain = new RooDataSet(DATASET_NAME, "Markov Chain", *fDataEntry,WEIGHT_NAME);
100}
101
102void MarkovChain::Add(RooArgSet& entry, double nllValue, double weight)
103{
104 if (fParameters == nullptr)
105 SetParameters(entry);
107 fNLL->setVal(nllValue);
108 //kbelasco: this is stupid, but some things might require it, so be doubly sure
109 fWeight->setVal(weight);
110 fChain->add(*fDataEntry, weight);
111 //fChain->add(*fDataEntry);
112}
113
115{
116 // Discards the first n accepted points.
117
118 if(fParameters == nullptr) SetParameters(*(RooArgSet*)otherChain.Get());
119 int counter = 0;
120 for( int i=0; i < otherChain.Size(); i++ ) {
121 RooArgSet* entry = (RooArgSet*)otherChain.Get(i);
122 counter += 1;
123 if( counter > burnIn ) {
124 AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
125 }
126 }
127}
128void MarkovChain::Add(MarkovChain& otherChain, double discardEntries)
129{
130 // Discards the first entries. This is different to the definition of
131 // burn-in used in the Bayesian calculator where the first n accepted
132 // terms from the proposal function are discarded.
133
134 if(fParameters == nullptr) SetParameters(*(RooArgSet*)otherChain.Get());
135 double counter = 0.0;
136 for( int i=0; i < otherChain.Size(); i++ ) {
137 RooArgSet* entry = (RooArgSet*)otherChain.Get(i);
138 counter += otherChain.Weight();
139 if( counter > discardEntries ) {
140 AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
141 }
142 }
143}
144
145void MarkovChain::AddFast(RooArgSet& entry, double nllValue, double weight)
146{
148 fNLL->setVal(nllValue);
149 //kbelasco: this is stupid, but some things might require it, so be doubly sure
150 fWeight->setVal(weight);
151 fChain->addFast(*fDataEntry, weight);
152 //fChain->addFast(*fDataEntry);
153}
154
156{
157 RooArgSet args;
158 if (whichVars == nullptr) {
159 //args.add(*fParameters);
160 //args.add(*fNLL);
161 args.add(*fDataEntry);
162 } else {
163 args.add(*whichVars);
164 }
165
167 //data = dynamic_cast<RooDataSet*>(fChain->reduce(args));
168 data = (RooDataSet*)fChain->reduce(args);
169
170 return data;
171}
172
174 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
175 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
176{
178 data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
179 return data;
180}
181
183{
184 RooArgSet args;
185 if (whichVars == nullptr) {
186 args.add(*fParameters);
187 //args.add(*fNLL);
188 //args.add(*fDataEntry);
189 } else {
190 args.add(*whichVars);
191 }
192
194 RooDataHist* hist = data->binnedClone();
195 delete data;
196
197 return hist;
198}
199
201 const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
202 const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
203{
205 RooDataHist* hist;
206 data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
207 hist = data->binnedClone();
208 delete data;
209
210 return hist;
211}
212
214{
215 RooArgList axes;
216 if (whichVars == nullptr)
217 axes.add(*fParameters);
218 else
219 axes.add(*whichVars);
220
221 Int_t dim = axes.getSize();
222 std::vector<double> min(dim);
223 std::vector<double> max(dim);
224 std::vector<Int_t> bins(dim);
225 std::vector<const char *> names(dim);
226 Int_t i = 0;
227 for (auto const *var : static_range_cast<RooRealVar *>(axes)) {
228 names[i] = var->GetName();
229 min[i] = var->getMin();
230 max[i] = var->getMax();
231 bins[i] = var->numBins();
232 ++i;
233 }
234
235 THnSparseF* sparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
236 dim, &bins[0], &min[0], &max[0]);
237
238 // kbelasco: it appears we need to call Sumw2() just to get the
239 // histogram to keep a running total of the weight so that Getsumw doesn't
240 // just return 0
241 sparseHist->Sumw2();
242
243 // Fill histogram
245 const RooArgSet* entry;
246 std::vector<double> x(dim);
247 for ( i = 0; i < size; i++) {
248 entry = fChain->get(i);
249
250 for (Int_t ii = 0; ii < dim; ii++) {
251 //LM: doing this is probably quite slow
252 x[ii] = entry->getRealValue( names[ii]);
253 sparseHist->Fill(x.data(), fChain->weight());
254 }
255 }
256
257 return sparseHist;
258}
259
260double MarkovChain::NLL(Int_t i) const
261{
262 // kbelasco: how to do this?
263 //fChain->get(i);
264 //return fNLL->getVal();
265 return fChain->get(i)->getRealValue(NLL_NAME);
266}
267
268double MarkovChain::NLL() const
269{
270 // kbelasco: how to do this?
271 //fChain->get();
272 //return fNLL->getVal();
273 return fChain->get()->getRealValue(NLL_NAME);
274}
275
277{
278 return fChain->weight();
279}
280
282{
283 fChain->get(i);
284 return fChain->weight();
285}
static const char * DATASET_NAME
Definition: MarkovChain.cxx:41
static const char * NLL_NAME
Definition: MarkovChain.cxx:40
static const char * DEFAULT_TITLE
Definition: MarkovChain.cxx:43
static const char * WEIGHT_NAME
Definition: MarkovChain.cxx:39
static const char * DEFAULT_NAME
Definition: MarkovChain.cxx:42
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition: TGX11.cxx:110
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:220
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
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 RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
RooAbsData * reduce(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg())
Create a reduced copy of this dataset.
Definition: RooAbsData.cxx:450
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:26
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
void add(const RooArgSet &row, double weight=1.0, double weightError=0.0) override
Add one ore more rows of data.
virtual void addFast(const RooArgSet &row, double weight=1.0, double weightError=0.0)
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
double weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:953
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
RooRealVar * fWeight
Definition: MarkovChain.h:119
virtual void AddFast(RooArgSet &entry, double nllValue, double weight=1.0)
add an entry to the chain ONLY IF you have constructed with parameters or called SetParameters
virtual double NLL() const
get the NLL value of the current (last indexed) entry
virtual void AddWithBurnIn(MarkovChain &otherChain, Int_t burnIn=0)
add another markov chain
RooArgSet * fDataEntry
Definition: MarkovChain.h:116
virtual void Add(RooArgSet &entry, double nllValue, double weight=1.0)
safely add an entry to the chain
virtual double NLL(Int_t i) const
get the NLL value of entry at position i
RooArgSet * fParameters
Definition: MarkovChain.h:115
virtual THnSparse * GetAsSparseHist(RooAbsCollection *whichVars=nullptr) const
Get a clone of the markov chain on which this interval is based as a sparse histogram.
RooDataSet * fChain
Definition: MarkovChain.h:117
virtual void SetParameters(RooArgSet &parameters)
set which of your parameters this chain should store
Definition: MarkovChain.cxx:77
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
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
Definition: MarkovChain.h:49
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:149
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:206
Efficient multidimensional histogram.
Definition: THnSparse.h:36
void Sumw2() override
Enable calculation of errors.
Definition: THnSparse.cxx:948
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
Namespace for the RooStats classes.
Definition: Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:63