Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
32
33using namespace RooFit;
34using namespace RooStats;
35
36static const char* WEIGHT_NAME = "weight_MarkovChain_local_";
37static const char* NLL_NAME = "nll_MarkovChain_local_";
38static const char* DATASET_NAME = "dataset_MarkovChain_local_";
39static const char* DEFAULT_NAME = "_markov_chain";
40static const char* DEFAULT_TITLE = "Markov Chain";
41
43
48
49MarkovChain::MarkovChain(const char *name, const char *title, RooArgSet &parameters) : TNamed(name, title)
50{
51 SetParameters(parameters);
52}
53
55{
56 delete fChain;
57 delete fParameters;
58 delete fDataEntry;
59
60 fParameters = new RooArgSet();
61 fParameters->addClone(parameters);
62
63 // kbelasco: consider setting fDataEntry = fChain->get()
64 // to see if that makes it possible to get values of variables without
65 // doing string comparison
66 RooRealVar nll(NLL_NAME, "-log Likelihood", 0);
67
68 fDataEntry = new RooArgSet();
69 fDataEntry->addClone(parameters);
70 fDataEntry->addClone(nll);
71 fNLL = static_cast<RooRealVar*>(fDataEntry->find(NLL_NAME));
72
74}
75
76void MarkovChain::Add(RooArgSet& entry, double nllValue, double weight)
77{
78 if (fParameters == nullptr)
82 //kbelasco: this is stupid, but some things might require it, so be doubly sure
83 fChain->add(*fDataEntry, weight);
84 //fChain->add(*fDataEntry);
85}
86
88{
89 // Discards the first n accepted points.
90
91 if(fParameters == nullptr) SetParameters(*const_cast<RooArgSet*>(otherChain.Get()));
92 int counter = 0;
93 for( int i=0; i < otherChain.Size(); i++ ) {
94 RooArgSet* entry = const_cast<RooArgSet*>(otherChain.Get(i));
95 counter += 1;
96 if( counter > burnIn ) {
97 AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
98 }
99 }
100}
102{
103 // Discards the first entries. This is different to the definition of
104 // burn-in used in the Bayesian calculator where the first n accepted
105 // terms from the proposal function are discarded.
106
107 if(fParameters == nullptr) SetParameters(*const_cast<RooArgSet*>(otherChain.Get()));
108 double counter = 0.0;
109 for( int i=0; i < otherChain.Size(); i++ ) {
110 RooArgSet* entry = const_cast<RooArgSet*>(otherChain.Get(i));
111 counter += otherChain.Weight();
112 if( counter > discardEntries ) {
113 AddFast( *entry, otherChain.NLL(), otherChain.Weight() );
114 }
115 }
116}
117
118void MarkovChain::AddFast(RooArgSet& entry, double nllValue, double weight)
119{
122 //kbelasco: this is stupid, but some things might require it, so be doubly sure
123 fChain->addFast(*fDataEntry, weight);
124 //fChain->addFast(*fDataEntry);
125}
126
128{
129 RooArgSet args;
130 if (whichVars == nullptr) {
131 //args.add(*fParameters);
132 //args.add(*fNLL);
133 args.add(*fDataEntry);
134 } else {
135 args.add(*whichVars);
136 }
137
138 return RooFit::makeOwningPtr<RooDataSet>(std::unique_ptr<RooAbsData>{fChain->reduce(args)});
139}
140
142{
143 RooArgSet args;
144 if (whichVars == nullptr) {
145 args.add(*fParameters);
146 //args.add(*fNLL);
147 //args.add(*fDataEntry);
148 } else {
149 args.add(*whichVars);
150 }
151
152 std::unique_ptr<RooAbsData> data{fChain->reduce(args)};
153 return RooFit::makeOwningPtr(std::unique_ptr<RooDataHist>{static_cast<RooDataSet&>(*data).binnedClone()});
154}
155
157{
158 RooArgList axes;
159 if (whichVars == nullptr) {
160 axes.add(*fParameters);
161 } else {
162 axes.add(*whichVars);
163 }
164
165 Int_t dim = axes.size();
166 std::vector<double> min(dim);
167 std::vector<double> max(dim);
168 std::vector<Int_t> bins(dim);
169 std::vector<const char *> names(dim);
170 Int_t i = 0;
171 for (auto const *var : static_range_cast<RooRealVar *>(axes)) {
172 names[i] = var->GetName();
173 min[i] = var->getMin();
174 max[i] = var->getMax();
175 bins[i] = var->numBins();
176 ++i;
177 }
178
179 THnSparseF* sparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
180 dim, &bins[0], &min[0], &max[0]);
181
182 // kbelasco: it appears we need to call Sumw2() just to get the
183 // histogram to keep a running total of the weight so that Getsumw doesn't
184 // just return 0
185 sparseHist->Sumw2();
186
187 // Fill histogram
189 const RooArgSet* entry;
190 std::vector<double> x(dim);
191 for ( i = 0; i < size; i++) {
192 entry = fChain->get(i);
193
194 for (Int_t ii = 0; ii < dim; ii++) {
195 //LM: doing this is probably quite slow
196 x[ii] = entry->getRealValue( names[ii]);
197 sparseHist->Fill(x.data(), fChain->weight());
198 }
199 }
200
201 return sparseHist;
202}
203
204double MarkovChain::NLL(Int_t i) const
205{
206 // kbelasco: how to do this?
207 //fChain->get(i);
208 //return fNLL->getVal();
209 return fChain->get(i)->getRealValue(NLL_NAME);
210}
211
212double MarkovChain::NLL() const
213{
214 // kbelasco: how to do this?
215 //fChain->get();
216 //return fNLL->getVal();
217 return fChain->get()->getRealValue(NLL_NAME);
218}
219
221{
222 return fChain->weight();
223}
224
226{
227 fChain->get(i);
228 return fChain->weight();
229}
static const char * DATASET_NAME
static const char * NLL_NAME
static const char * DEFAULT_TITLE
static const char * WEIGHT_NAME
static const char * DEFAULT_NAME
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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:222
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
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.
RooFit::OwningPtr< RooAbsData > reduce(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a reduced copy of this dataset.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
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:24
Container class to hold unbinned data.
Definition RooDataSet.h:34
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
RooFit::OwningPtr< RooDataHist > binnedClone(const char *newName=nullptr, const char *newTitle=nullptr) const
Return binned clone of this dataset.
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.
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
double weight() const override
Return event weight of current event.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
Stores the steps in a Markov Chain of points.
Definition MarkovChain.h:26
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
virtual RooFit::OwningPtr< RooDataHist > GetAsDataHist(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataHist whose entries contain the values of whichVars.
RooArgSet * fDataEntry
Definition MarkovChain.h:95
virtual void Add(RooArgSet &entry, double nllValue, double weight=1.0)
safely add an entry to the chain
virtual RooFit::OwningPtr< RooDataSet > GetAsDataSet(RooArgSet *whichVars=nullptr) const
get this MarkovChain as a RooDataSet whose entries contain the values of whichVars.
RooArgSet * fParameters
Definition MarkovChain.h:94
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:96
virtual void SetParameters(RooArgSet &parameters)
set which of your parameters this chain should store
virtual double Weight() const
get the weight of the current (last indexed) entry
Efficient multidimensional histogram.
Definition THnSparse.h:37
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
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 CodegenImpl.h:64
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
Namespace for the RooStats classes.
Definition CodegenImpl.h:58
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)