Logo ROOT   6.10/09
Reference Guide
PdfProposal.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 /** \class RooStats::PdfProposal
15  \ingroup Roostats
16 
17 PdfProposal is a concrete implementation of the ProposalFunction interface.
18 It proposes points across the parameter space in the distribution of the
19 given PDF.
20 
21 To make Propose(xPrime, x) dependent on x, configure with
22 PdfProposal::AddMapping(varToUpdate, valueToUse). For example, suppose we have:
23 
24 ~~~{.cpp}
25 // our parameter
26 RooRealVar p("p", "p", 5, 0, 10);
27 
28 // create mean and sigma for gaussian proposal function
29 RooRealVar meanP("meanP", "meanP", 0, 10);
30 RooRealVar sigma("sigma", "sigma", 1, 0, 5);
31 RooGaussian pGaussian("pGaussian", "pGaussian", p, meanP, sigma);
32 
33 // configure proposal function
34 PdfProposal pdfProposal(pGaussian);
35 pdfProposal.AddMapping(meanP, p); // each call of Propose(xPrime, x), meanP in
36  // the proposal function will be updated to
37  // the value of p in x. this will center the
38  // proposal function about x's p when
39  // proposing for xPrime
40 
41 // To improve performance, PdfProposal has the ability to cache a specified
42 // number of proposals. If you don't call this function, the default cache size
43 // is 1, which can be slow.
44 pdfProposal.SetCacheSize(desiredCacheSize);
45 ~~~
46 
47 PdfProposal currently uses a fixed cache size. Adaptive caching methods are in the works
48 for future versions.
49 */
50 
51 
52 #include "Rtypes.h"
53 
54 #include "RooStats/PdfProposal.h"
55 #include "RooStats/RooStatsUtils.h"
56 #include "RooArgSet.h"
57 #include "RooDataSet.h"
58 #include "RooAbsPdf.h"
59 #include "RooMsgService.h"
60 #include "RooRealVar.h"
61 #include "TIterator.h"
62 
63 #include <map>
64 
66 
67 using namespace RooFit;
68 using namespace RooStats;
69 using namespace std;
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// By default, PdfProposal does NOT own the PDF that serves as the
73 /// proposal density function
74 
75 PdfProposal::PdfProposal() : ProposalFunction()
76 {
77  fPdf = NULL;
78  fOwnsPdf = kFALSE;
79  fCacheSize = 1;
80  fCachePosition = 0;
81  fCache = NULL;
82 }
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 /// By default, PdfProposal does NOT own the PDF that serves as the
86 /// proposal density function
87 
89 {
90  fPdf = &pdf;
91  fOwnsPdf = kFALSE;
92  fCacheSize = 1;
93  fCachePosition = 0;
94  fCache = NULL;
95 }
96 
97 ////////////////////////////////////////////////////////////////////////////////
98 /// determine whether these two RooArgSets represent the same point
99 
101 {
102  if (x1.equals(x2)) {
103  TIterator* it = x1.createIterator();
104  RooRealVar* r;
105  while ((r = (RooRealVar*)it->Next()) != NULL)
106  if (r->getVal() != x2.getRealValue(r->GetName())) {
107  delete it;
108  return kFALSE;
109  }
110  delete it;
111  return kTRUE;
112  }
113  return kFALSE;
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Populate xPrime with a new proposed point
118 
120 {
121  if (fLastX.getSize() == 0) {
122  // fLastX not yet initialized
123  fLastX.addClone(x);
124  // generate initial cache
126  if (fMap.size() > 0) {
127  for (fIt = fMap.begin(); fIt != fMap.end(); fIt++)
128  fIt->first->setVal(fIt->second->getVal(&x));
129  }
130  fCache = fPdf->generate(xPrime, fCacheSize);
131  }
132 
133  Bool_t moved = false;
134  if (fMap.size() > 0) {
135  moved = !Equals(fLastX, x);
136 
137  // if we've moved, set the values of the variables in the PDF to the
138  // corresponding values of the variables in x, according to the
139  // mappings (i.e. let the variables in x set the given values for the
140  // PDF that will generate xPrime)
141  if (moved) {
142  // update the pdf parameters
144 
145  for (fIt = fMap.begin(); fIt != fMap.end(); fIt++)
146  fIt->first->setVal(fIt->second->getVal(&x));
147 
148  // save the new x in fLastX
150  }
151  }
152 
153  // generate new cache if necessary
154  if (moved || fCachePosition >= fCacheSize) {
155  delete fCache;
156  fCache = fPdf->generate(xPrime, fCacheSize);
157  fCachePosition = 0;
158  }
159 
160  const RooArgSet* proposal = fCache->get(fCachePosition);
161  fCachePosition++;
162  RooStats::SetParameters(proposal, &xPrime);
163 }
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 /// Determine whether or not the proposal density is symmetric for
167 /// points x1 and x2 - that is, whether the probabilty of reaching x2
168 /// from x1 is equal to the probability of reaching x1 from x2
169 
171 {
172  // kbelasco: is there a better way to do this?
173  return false;
174 }
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 /// Return the probability of proposing the point x1 given the starting
178 /// point x2
179 
181 {
183  for (fIt = fMap.begin(); fIt != fMap.end(); fIt++)
184  fIt->first->setVal(fIt->second->getVal(&x2));
185  RooArgSet* temp = fPdf->getObservables(x1);
186  RooStats::SetParameters(&x1, temp);
187  delete temp;
188  return fPdf->getVal(&x1); // could just as well use x2
189 }
190 
191 ////////////////////////////////////////////////////////////////////////////////
192 /// specify a mapping between a parameter of the proposal function and
193 /// a parameter of interest. this mapping is used to set the value of
194 /// proposalParam equal to the value of update to determine the
195 /// proposal function.
196 /// proposalParam is a parameter of the proposal function that must
197 /// be set to the value of update (from the current point) in order to
198 /// propose a new point.
199 
201 {
203  if (update.getParameters((RooAbsData*)NULL)->getSize() == 0)
204  fMaster.add(update);
205  fMap.insert(pair<RooRealVar*, RooAbsReal*>(&proposalParam, &update));
206 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
TIterator * createIterator(Bool_t dir=kIterForward) 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
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:58
virtual Bool_t Equals(RooArgSet &x1, RooArgSet &x2)
whether we own the proposal density function
bool Bool_t
Definition: RtypesCore.h:59
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically named contents.
STL namespace.
#define NULL
Definition: RtypesCore.h:88
Iterator abstract base class.
Definition: TIterator.h:30
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
virtual Double_t GetProposalDensity(RooArgSet &x1, RooArgSet &x2)
Return the probability of proposing the point x1 given the starting point x2.
virtual void AddMapping(RooRealVar &proposalParam, RooAbsReal &update)
specify a mapping between a parameter of the proposal function and a parameter of interest...
Int_t fCacheSize
the last point we were at
Definition: PdfProposal.h:108
std::map< RooRealVar *, RooAbsReal * >::iterator fIt
map of values in pdf to update
Definition: PdfProposal.h:106
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
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
Int_t getSize() const
TRandom2 r(17)
RooDataSet * fCache
our position in the cached proposal data set
Definition: PdfProposal.h:110
RooArgSet fLastX
pdf iterator
Definition: PdfProposal.h:107
static void update(gsl_integration_workspace *workspace, double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
Bool_t fOwnsPdf
pointers to master variables needed for updates
Definition: PdfProposal.h:112
const Bool_t kFALSE
Definition: RtypesCore.h:92
Namespace for the RooStats classes.
Definition: Asimov.h:20
virtual void Propose(RooArgSet &xPrime, RooArgSet &x)
Populate xPrime with a new proposed point.
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:560
RooArgSet fMaster
the cached proposal data set
Definition: PdfProposal.h:111
PdfProposal is a concrete implementation of the ProposalFunction interface.
Definition: PdfProposal.h:30
std::map< RooRealVar *, RooAbsReal * > fMap
the proposal density function
Definition: PdfProposal.h:105
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
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
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
Generate a new dataset containing the specified variables with events sampled from our distribution...
Definition: RooAbsPdf.cxx:1702
virtual TObject * Next()=0
PdfProposal()
By default, PdfProposal does NOT own the PDF that serves as the proposal density function.
Definition: PdfProposal.cxx:75
virtual Bool_t IsSymmetric(RooArgSet &x1, RooArgSet &x2)
Determine whether or not the proposal density is symmetric for points x1 and x2 - that is...
const Bool_t kTRUE
Definition: RtypesCore.h:91
Int_t fCachePosition
how many points to generate each time
Definition: PdfProposal.h:109