Logo ROOT   6.08/07
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 ////////////////////////////////////////////////////////////////////////////////
13 
14 //_________________________________________________
15 
16 #ifndef ROOT_Rtypes
17 #include "Rtypes.h"
18 #endif
19 
20 #ifndef ROOT_TNamed
21 #include "TNamed.h"
22 #endif
23 #ifndef ROOSTATS_MarkovChain
24 #include "RooStats/MarkovChain.h"
25 #endif
26 #ifndef ROO_DATA_SET
27 #include "RooDataSet.h"
28 #endif
29 #ifndef ROO_ARG_SET
30 #include "RooArgSet.h"
31 #endif
32 #ifndef ROO_REAL_VAR
33 #include "RooRealVar.h"
34 #endif
35 #ifndef RooStats_RooStatsUtils
36 #include "RooStats/RooStatsUtils.h"
37 #endif
38 #ifndef ROO_DATA_HIST
39 #include "RooDataHist.h"
40 #endif
41 #ifndef ROOT_THnSparse
42 #include "THnSparse.h"
43 #endif
44 
45 using namespace std;
46 
48 
49 using namespace RooFit;
50 using namespace RooStats;
51 
52 static const char* WEIGHT_NAME = "weight_MarkovChain_local_";
53 static const char* NLL_NAME = "nll_MarkovChain_local_";
54 static const char* DATASET_NAME = "dataset_MarkovChain_local_";
55 static const char* DEFAULT_NAME = "_markov_chain";
56 static const char* DEFAULT_TITLE = "Markov Chain";
57 
58 MarkovChain::MarkovChain() :
60 {
61  fParameters = NULL;
62  fDataEntry = NULL;
63  fChain = NULL;
64  fNLL = NULL;
65  fWeight = NULL;
66 }
67 
70 {
71  fParameters = NULL;
72  fDataEntry = NULL;
73  fChain = NULL;
74  fNLL = NULL;
75  fWeight = NULL;
76  SetParameters(parameters);
77 }
78 
79 MarkovChain::MarkovChain(const char* name, const char* title,
80  RooArgSet& parameters) : TNamed(name, title)
81 {
82  fParameters = NULL;
83  fDataEntry = NULL;
84  fChain = NULL;
85  fNLL = NULL;
86  fWeight = NULL;
87  SetParameters(parameters);
88 }
89 
91 {
92  delete fChain;
93  delete fParameters;
94  delete fDataEntry;
95 
96  fParameters = new RooArgSet();
97  fParameters->addClone(parameters);
98 
99  // kbelasco: consider setting fDataEntry = fChain->get()
100  // to see if that makes it possible to get values of variables without
101  // doing string comparison
102  RooRealVar nll(NLL_NAME, "-log Likelihood", 0);
103  RooRealVar weight(WEIGHT_NAME, "weight", 0);
104 
105  fDataEntry = new RooArgSet();
106  fDataEntry->addClone(parameters);
107  fDataEntry->addClone(nll);
108  fDataEntry->addClone(weight);
111 
112  fChain = new RooDataSet(DATASET_NAME, "Markov Chain", *fDataEntry,WEIGHT_NAME);
113 }
114 
115 void MarkovChain::Add(RooArgSet& entry, Double_t nllValue, Double_t weight)
116 {
117  if (fParameters == NULL)
118  SetParameters(entry);
120  fNLL->setVal(nllValue);
121  //kbelasco: this is stupid, but some things might require it, so be doubly sure
122  fWeight->setVal(weight);
123  fChain->add(*fDataEntry, weight);
124  //fChain->add(*fDataEntry);
125 }
126 
128 {
129  // Discards the first n accepted points.
130 
131  if(fParameters == NULL) SetParameters(*(RooArgSet*)otherChain.Get());
132  int counter = 0;
133  for( int i=0; i < otherChain.Size(); i++ ) {
134  RooArgSet* entry = (RooArgSet*)otherChain.Get(i);
135  counter += 1;
136  if( counter > burnIn ) {
137  AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
138  }
139  }
140 }
141 void MarkovChain::Add(MarkovChain& otherChain, Double_t discardEntries)
142 {
143  // Discards the first entries. This is different to the definition of
144  // burn-in used in the Bayesian calculator where the first n accepted
145  // terms from the proposal function are discarded.
146 
147  if(fParameters == NULL) SetParameters(*(RooArgSet*)otherChain.Get());
148  double counter = 0.0;
149  for( int i=0; i < otherChain.Size(); i++ ) {
150  RooArgSet* entry = (RooArgSet*)otherChain.Get(i);
151  counter += otherChain.Weight();
152  if( counter > discardEntries ) {
153  AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
154  }
155  }
156 }
157 
158 void MarkovChain::AddFast(RooArgSet& entry, Double_t nllValue, Double_t weight)
159 {
161  fNLL->setVal(nllValue);
162  //kbelasco: this is stupid, but some things might require it, so be doubly sure
163  fWeight->setVal(weight);
164  fChain->addFast(*fDataEntry, weight);
165  //fChain->addFast(*fDataEntry);
166 }
167 
169 {
170  RooArgSet args;
171  if (whichVars == NULL) {
172  //args.add(*fParameters);
173  //args.add(*fNLL);
174  args.add(*fDataEntry);
175  } else {
176  args.add(*whichVars);
177  }
178 
179  RooDataSet* data;
180  //data = dynamic_cast<RooDataSet*>(fChain->reduce(args));
181  data = (RooDataSet*)fChain->reduce(args);
182 
183  return data;
184 }
185 
187  const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
188  const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
189 {
190  RooDataSet* data;
191  data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
192  return data;
193 }
194 
196 {
197  RooArgSet args;
198  if (whichVars == NULL) {
199  args.add(*fParameters);
200  //args.add(*fNLL);
201  //args.add(*fDataEntry);
202  } else {
203  args.add(*whichVars);
204  }
205 
206  RooDataSet* data = (RooDataSet*)fChain->reduce(args);
207  RooDataHist* hist = data->binnedClone();
208  delete data;
209 
210  return hist;
211 }
212 
214  const RooCmdArg& arg3, const RooCmdArg& arg4, const RooCmdArg& arg5,
215  const RooCmdArg& arg6, const RooCmdArg& arg7, const RooCmdArg& arg8) const
216 {
217  RooDataSet* data;
218  RooDataHist* hist;
219  data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
220  hist = data->binnedClone();
221  delete data;
222 
223  return hist;
224 }
225 
227 {
228  RooArgList axes;
229  if (whichVars == NULL)
230  axes.add(*fParameters);
231  else
232  axes.add(*whichVars);
233 
234  Int_t dim = axes.getSize();
235  std::vector<Double_t> min(dim);
236  std::vector<Double_t> max(dim);
237  std::vector<Int_t> bins(dim);
238  std::vector<const char *> names(dim);
239  TIterator* it = axes.createIterator();
240  for (Int_t i = 0; i < dim; i++) {
241  RooRealVar * var = dynamic_cast<RooRealVar*>(it->Next() );
242  assert(var != 0);
243  names[i] = var->GetName();
244  min[i] = var->getMin();
245  max[i] = var->getMax();
246  bins[i] = var->numBins();
247  }
248 
249  THnSparseF* sparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
250  dim, &bins[0], &min[0], &max[0]);
251 
252  // kbelasco: it appears we need to call Sumw2() just to get the
253  // histogram to keep a running total of the weight so that Getsumw doesn't
254  // just return 0
255  sparseHist->Sumw2();
256 
257  // Fill histogram
258  Int_t size = fChain->numEntries();
259  const RooArgSet* entry;
260  Double_t* x = new Double_t[dim];
261  for (Int_t i = 0; i < size; i++) {
262  entry = fChain->get(i);
263  it->Reset();
264  for (Int_t ii = 0; ii < dim; ii++) {
265  //LM: doing this is probably quite slow
266  x[ii] = entry->getRealValue( names[ii]);
267  sparseHist->Fill(x, fChain->weight());
268  }
269  }
270  delete[] x;
271  delete it;
272 
273  return sparseHist;
274 }
275 
277 {
278  // kbelasco: how to do this?
279  //fChain->get(i);
280  //return fNLL->getVal();
281  return fChain->get(i)->getRealValue(NLL_NAME);
282 }
283 
285 {
286  // kbelasco: how to do this?
287  //fChain->get();
288  //return fNLL->getVal();
289  return fChain->get()->getRealValue(NLL_NAME);
290 }
291 
293 {
294  return fChain->weight();
295 }
296 
298 {
299  fChain->get(i);
300  return fChain->weight();
301 }
virtual Double_t getMin(const char *name=0) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
TIterator * createIterator(Bool_t dir=kIterForward) const
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:56
virtual Double_t getMax(const char *name=0) const
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:86
virtual void Reset()=0
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:153
RooArgSet * fDataEntry
Definition: MarkovChain.h:139
static const char * DATASET_NAME
Definition: MarkovChain.cxx:54
virtual Int_t numBins(const char *rangeName=0) const
RooArgSet * fParameters
Definition: MarkovChain.h:138
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:69
virtual Int_t Size() const
get the number of steps in the chain
Definition: MarkovChain.h:72
int Int_t
Definition: RtypesCore.h:41
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:343
static const char * DEFAULT_NAME
Definition: MarkovChain.cxx:55
STL namespace.
Iterator abstract base class.
Definition: TIterator.h:32
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
RooDataSet 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:981
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:74
Efficient multidimensional histogram.
Definition: THnSparse.h:52
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
RooDataSet * fChain
Definition: MarkovChain.h:140
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.
virtual Double_t weight() const
Return event weight of current event.
static const char * WEIGHT_NAME
Definition: MarkovChain.cxx:52
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:233
virtual void AddWithBurnIn(MarkovChain &otherChain, Int_t burnIn=0)
add another markov chain
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
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:94
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
Int_t getSize() const
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:29
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual void add(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...
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:53
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:279
RooAbsArg * find(const char *name) const
Find object with given name in list.
double Double_t
Definition: RtypesCore.h:55
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:90
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:142
static const char * NLL_NAME
Definition: MarkovChain.cxx:53
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:527
virtual TObject * Next()=0
#define NULL
Definition: Rtypes.h:82
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:219
char name[80]
Definition: TGX11.cxx:109
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:27