Logo ROOT   master
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 
15 Stores the steps in a Markov Chain of points. Allows user to access the
16 weight and NLL value (if applicable) with which a point was added to the
17 MarkovChain.
18 
19 */
20 
21 #include "Rtypes.h"
22 
23 #include "TNamed.h"
24 #include "RooStats/MarkovChain.h"
25 #include "RooDataSet.h"
26 #include "RooArgSet.h"
27 #include "RooRealVar.h"
28 #include "RooStats/RooStatsUtils.h"
29 #include "RooDataHist.h"
30 #include "THnSparse.h"
31 
32 using namespace std;
33 
35 
36 using namespace RooFit;
37 using namespace RooStats;
38 
39 static const char* WEIGHT_NAME = "weight_MarkovChain_local_";
40 static const char* NLL_NAME = "nll_MarkovChain_local_";
41 static const char* DATASET_NAME = "dataset_MarkovChain_local_";
42 static const char* DEFAULT_NAME = "_markov_chain";
43 static const char* DEFAULT_TITLE = "Markov Chain";
44 
45 MarkovChain::MarkovChain() :
47 {
48  fParameters = NULL;
49  fDataEntry = NULL;
50  fChain = NULL;
51  fNLL = NULL;
52  fWeight = NULL;
53 }
54 
57 {
58  fParameters = NULL;
59  fDataEntry = NULL;
60  fChain = NULL;
61  fNLL = NULL;
62  fWeight = NULL;
63  SetParameters(parameters);
64 }
65 
66 MarkovChain::MarkovChain(const char* name, const char* title,
67  RooArgSet& parameters) : TNamed(name, title)
68 {
69  fParameters = NULL;
70  fDataEntry = NULL;
71  fChain = NULL;
72  fNLL = NULL;
73  fWeight = NULL;
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 
102 void MarkovChain::Add(RooArgSet& entry, Double_t nllValue, Double_t weight)
103 {
104  if (fParameters == NULL)
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 == NULL) 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 }
128 void MarkovChain::Add(MarkovChain& otherChain, Double_t 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 == NULL) 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 
145 void MarkovChain::AddFast(RooArgSet& entry, Double_t nllValue, Double_t 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 == NULL) {
159  //args.add(*fParameters);
160  //args.add(*fNLL);
161  args.add(*fDataEntry);
162  } else {
163  args.add(*whichVars);
164  }
165 
166  RooDataSet* data;
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 {
177  RooDataSet* data;
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 == NULL) {
186  args.add(*fParameters);
187  //args.add(*fNLL);
188  //args.add(*fDataEntry);
189  } else {
190  args.add(*whichVars);
191  }
192 
193  RooDataSet* data = (RooDataSet*)fChain->reduce(args);
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 {
204  RooDataSet* data;
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 == NULL)
217  axes.add(*fParameters);
218  else
219  axes.add(*whichVars);
220 
221  Int_t dim = axes.getSize();
222  std::vector<Double_t> min(dim);
223  std::vector<Double_t> max(dim);
224  std::vector<Int_t> bins(dim);
225  std::vector<const char *> names(dim);
226  TIterator* it = axes.createIterator();
227  for (Int_t i = 0; i < dim; i++) {
228  RooRealVar * var = dynamic_cast<RooRealVar*>(it->Next() );
229  assert(var != 0);
230  names[i] = var->GetName();
231  min[i] = var->getMin();
232  max[i] = var->getMax();
233  bins[i] = var->numBins();
234  }
235 
236  THnSparseF* sparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
237  dim, &bins[0], &min[0], &max[0]);
238 
239  // kbelasco: it appears we need to call Sumw2() just to get the
240  // histogram to keep a running total of the weight so that Getsumw doesn't
241  // just return 0
242  sparseHist->Sumw2();
243 
244  // Fill histogram
245  Int_t size = fChain->numEntries();
246  const RooArgSet* entry;
247  Double_t* x = new Double_t[dim];
248  for (Int_t i = 0; i < size; i++) {
249  entry = fChain->get(i);
250  it->Reset();
251  for (Int_t ii = 0; ii < dim; ii++) {
252  //LM: doing this is probably quite slow
253  x[ii] = entry->getRealValue( names[ii]);
254  sparseHist->Fill(x, fChain->weight());
255  }
256  }
257  delete[] x;
258  delete it;
259 
260  return sparseHist;
261 }
262 
264 {
265  // kbelasco: how to do this?
266  //fChain->get(i);
267  //return fNLL->getVal();
268  return fChain->get(i)->getRealValue(NLL_NAME);
269 }
270 
272 {
273  // kbelasco: how to do this?
274  //fChain->get();
275  //return fNLL->getVal();
276  return fChain->get()->getRealValue(NLL_NAME);
277 }
278 
280 {
281  return fChain->weight();
282 }
283 
285 {
286  fChain->get(i);
287  return fChain->weight();
288 }
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
static const char * DEFAULT_TITLE
Definition: MarkovChain.cxx:43
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
virtual void Reset()=0
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:144
RooArgSet * fDataEntry
Definition: MarkovChain.h:116
static const char * DATASET_NAME
Definition: MarkovChain.cxx:41
virtual Int_t numBins(const char *rangeName=0) const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooArgSet * fParameters
Definition: MarkovChain.h:115
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:65
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:49
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:381
static const char * DEFAULT_NAME
Definition: MarkovChain.cxx:42
STL namespace.
Iterator abstract base class.
Definition: TIterator.h:30
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
Double_t x[n]
Definition: legend1.C:17
RooDataHist * binnedClone(const char *newName=0, const char *newTitle=0) const
Return binned clone of this dataset.
Definition: RooDataSet.cxx:951
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
Efficient multidimensional histogram.
Definition: THnSparse.h:36
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
RooDataSet * fChain
Definition: MarkovChain.h:117
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void Add(RooArgSet &entry, Double_t nllValue, Double_t weight=1.0)
safely add an entry to the chain
virtual THnSparse * GetAsSparseHist(RooAbsCollection *whichVars=NULL) const
Get a clone of the markov chain on which this interval is based as a sparse histogram.
static const char * WEIGHT_NAME
Definition: MarkovChain.cxx:39
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:220
virtual void AddWithBurnIn(MarkovChain &otherChain, Int_t burnIn=0)
add another markov chain
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void addClone(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling addOwned() for each element in the source...
Definition: RooArgSet.h:96
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:225
Int_t getSize() const
virtual Double_t weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:973
virtual Double_t NLL(Int_t i) const
get the NLL value of entry at position i
virtual void AddFast(RooArgSet &entry, Double_t nllValue, Double_t weight=1.0)
add an entry to the chain ONLY IF you have constructed with parameters or called SetParameters ...
void Sumw2()
Enable calculation of errors.
Definition: THnSparse.cxx:949
virtual Double_t NLL() const
get the NLL value of the current (last indexed) entry
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
virtual void addFast(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0)
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define ClassImp(name)
Definition: Rtypes.h:361
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
virtual void SetParameters(RooArgSet &parameters)
set which of your parameters this chain should store
Definition: MarkovChain.cxx:77
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects...
RooRealVar * fWeight
Definition: MarkovChain.h:119
static const char * NLL_NAME
Definition: MarkovChain.cxx:40
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:472
virtual TObject * Next()=0
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:206
char name[80]
Definition: TGX11.cxx:109
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:28