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 
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 }
RooCmdArg
Definition: RooCmdArg.h:27
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
RooAbsRealLValue::getMax
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
Definition: RooAbsRealLValue.h:89
DEFAULT_TITLE
static const char * DEFAULT_TITLE
Definition: MarkovChain.cxx:43
RooArgSet::getRealValue
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:474
RooArgSet.h
RooStats::MarkovChain::Size
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:61
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooStats::MarkovChain::fDataEntry
RooArgSet * fDataEntry
Definition: MarkovChain.h:128
DATASET_NAME
static const char * DATASET_NAME
Definition: MarkovChain.cxx:41
THnBase::Fill
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:144
TNamed.h
RooStats::MarkovChain::Get
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:63
RooArgList
Definition: RooArgList.h:21
RooStats::MarkovChain
Definition: MarkovChain.h:36
RooStats::MarkovChain::GetAsDataHist
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
Definition: MarkovChain.cxx:182
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
x
Double_t x[n]
Definition: legend1.C:17
RooStats::MarkovChain::MarkovChain
MarkovChain()
Definition: MarkovChain.cxx:45
RooArgSet::add
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
DEFAULT_NAME
static const char * DEFAULT_NAME
Definition: MarkovChain.cxx:42
RooStats::MarkovChain::fParameters
RooArgSet * fParameters
Definition: MarkovChain.h:127
RooDataSet::addFast
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 'data' argset, to the data set.
Definition: RooDataSet.cxx:1266
RooDataSet.h
TIterator
Definition: TIterator.h:30
RooAbsRealLValue::getMin
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Definition: RooAbsRealLValue.h:86
RooStats::MarkovChain::GetAsSparseHist
virtual THnSparse * GetAsSparseHist(RooAbsCollection *whichVars=NULL) const
Get a clone of the markov chain on which this interval is based as a sparse histogram.
Definition: MarkovChain.cxx:213
RooDataSet::get
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Definition: RooDataSet.cxx:1038
THnSparse
Definition: THnSparse.h:36
THnSparseT
Definition: THnSparse.h:206
RooStats::MarkovChain::fChain
RooDataSet * fChain
Definition: MarkovChain.h:129
RooDataSet::add
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 'data' argset, to the data set.
Definition: RooDataSet.cxx:1164
RooDataHist
Definition: RooDataHist.h:39
MarkovChain.h
RooAbsRealLValue::numBins
virtual Int_t numBins(const char *rangeName=0) const
Definition: RooAbsRealLValue.h:54
RooStats::MarkovChain::AddFast
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
Definition: MarkovChain.cxx:145
WEIGHT_NAME
static const char * WEIGHT_NAME
Definition: MarkovChain.cxx:39
RooFit
Definition: RooCFunction1Binding.h:29
RooStats::MarkovChain::fNLL
RooRealVar * fNLL
Definition: MarkovChain.h:130
TNamed
Definition: TNamed.h:29
RooDataHist.h
RooAbsCollection::createIterator
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
Definition: RooAbsCollection.h:118
RooStats::MarkovChain::Add
virtual void Add(RooArgSet &entry, Double_t nllValue, Double_t weight=1.0)
safely add an entry to the chain
Definition: MarkovChain.cxx:102
RooAbsCollection
Definition: RooAbsCollection.h:30
RooAbsCollection::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:437
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
RooStats::MarkovChain::AddWithBurnIn
virtual void AddWithBurnIn(MarkovChain &otherChain, Int_t burnIn=0)
add another markov chain
Definition: MarkovChain.cxx:114
RooStats::MarkovChain::Weight
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
Definition: MarkovChain.cxx:279
RooRealVar.h
RooStats::MarkovChain::fWeight
RooRealVar * fWeight
Definition: MarkovChain.h:131
TIterator::Next
virtual TObject * Next()=0
RooStats::MarkovChain::NLL
virtual Double_t NLL() const
get the NLL value of the current (last indexed) entry
Definition: MarkovChain.cxx:271
TIterator::Reset
virtual void Reset()=0
RooStats::MarkovChain::NLL
virtual Double_t NLL(Int_t i) const
get the NLL value of entry at position i
Definition: MarkovChain.cxx:263
RooStats::MarkovChain::GetAsDataSet
virtual RooDataSet * GetAsDataSet(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
Definition: MarkovChain.cxx:155
RooStatsUtils.h
RooStats::MarkovChain::SetParameters
virtual void SetParameters(RooArgSet &parameters)
set which of your parameters this chain should store
Definition: MarkovChain.cxx:77
THnSparse.h
Double_t
double Double_t
Definition: RtypesCore.h:59
NLL_NAME
static const char * NLL_NAME
Definition: MarkovChain.cxx:40
RooStats
Definition: Asimov.h:19
RooAbsData::reduce
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:382
RooDataSet::binnedClone
RooDataHist * binnedClone(const char *newName=0, const char *newTitle=0) const
Return binned clone of this dataset.
Definition: RooDataSet.cxx:972
THnSparseF
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:220
RooStats::SetParameters
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:65
name
char name[80]
Definition: TGX11.cxx:110
RooDataSet
Definition: RooDataSet.h:33
RooArgSet::addClone
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
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooRealVar
Definition: RooRealVar.h:36
Rtypes.h
THnSparse::Sumw2
void Sumw2()
Enable calculation of errors.
Definition: THnSparse.cxx:948
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooArgSet
Definition: RooArgSet.h:28
RooDataSet::weight
virtual Double_t weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:994
int