Logo ROOT   6.18/05
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 = 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
66MarkovChain::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
102void 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}
128void 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
145void 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}
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
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
THnSparseT< TArrayF > THnSparseF
Definition: THnSparse.h:220
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
Int_t getSize() const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
RooAbsArg * find(const char *name) const
Find object with given name in list.
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
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Int_t numBins(const char *rangeName=0) const
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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 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 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
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:27
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event 'index'.
virtual Double_t weight() const
Return event weight of current event.
Definition: RooDataSet.cxx:977
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 'data' argset, to the data set.
RooDataHist * binnedClone(const char *newName=0, const char *newTitle=0) const
Return binned clone of this dataset.
Definition: RooDataSet.cxx:955
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.
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:233
Stores the steps in a Markov Chain of points.
Definition: MarkovChain.h:30
RooRealVar * fWeight
Definition: MarkovChain.h:119
virtual Double_t 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
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
RooArgSet * fDataEntry
Definition: MarkovChain.h:116
virtual Double_t Weight() const
get the weight of the current (last indexed) entry
virtual RooDataHist * GetAsDataHist(RooArgSet *whichVars=NULL) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
RooArgSet * fParameters
Definition: MarkovChain.h:115
virtual void Add(RooArgSet &entry, Double_t nllValue, Double_t weight=1.0)
safely add an entry to the chain
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 Double_t NLL(Int_t i) const
get the NLL value of entry at position i
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
Definition: MarkovChain.h:51
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 RooDataSet * GetAsDataSet(RooArgSet *whichVars=NULL) 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:144
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:206
Efficient multidimensional histogram.
Definition: THnSparse.h:36
void Sumw2()
Enable calculation of errors.
Definition: THnSparse.cxx:949
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
Template specialisation used in RooAbsArg:
Namespace for the RooStats classes.
Definition: Asimov.h:20
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:58