BayesianAnalysisToolkit  0.9.3
Classes | Public Member Functions | Protected Attributes | Private Types | Private Member Functions | Private Attributes | List of all members
BCEngineMCMC Class Reference

An engine class for Markov Chain Monte Carlo. More...

#include <BCEngineMCMC.h>

Inheritance diagram for BCEngineMCMC:
Inheritance graph
[legend]
Collaboration diagram for BCEngineMCMC:
Collaboration graph
[legend]

Classes

struct  MCMCThreadLocalStorage
 

Public Types

Enumerators
enum  Precision { kLow, kMedium, kHigh, kVeryHigh }
 

Public Member Functions

void WriteMarkovChain (bool flag)
 
Constructors and destructors
 BCEngineMCMC ()
 
 BCEngineMCMC (const BCEngineMCMC &enginemcmc)
 
virtual ~BCEngineMCMC ()
 
Assignment operators
BCEngineMCMC & operator= (const BCEngineMCMC &engineMCMC)
 
Getters
unsigned MCMCGetNChains () const
 
unsigned MCMCGetNLag () const
 
const std::vector< unsigned > & MCMCGetNIterations () const
 
unsigned MCMCGetCurrentIteration () const
 
unsigned MCMCGetCurrentChain () const
 
unsigned MCMCGetNIterationsConvergenceGlobal () const
 
bool MCMCGetFlagConvergenceGlobal () const
 
unsigned MCMCGetNIterationsMax () const
 
unsigned MCMCGetNIterationsRun () const
 
unsigned MCMCGetNIterationsPreRunMin () const
 
unsigned MCMCGetNIterationsUpdate () const
 
unsigned MCMCGetNIterationsUpdateMax () const
 
std::vector< int > MCMCGetNTrialsTrue () const
 
int MCMCGetNTrials () const
 
const std::vector< double > & MCMCGetprobMean () const
 
const std::vector< double > & MCMCGetVariance () const
 
const std::vector< double > & MCMCGetTrialFunctionScaleFactor () const
 
std::vector< double > MCMCGetTrialFunctionScaleFactor (unsigned ichain) const
 
double MCMCGetTrialFunctionScaleFactor (unsigned ichain, unsigned ipar)
 
const std::vector< double > & MCMCGetx () const
 
std::vector< double > MCMCGetx (unsigned ichain)
 
double MCMCGetx (unsigned ichain, unsigned ipar) const
 
const std::vector< double > & MCMCGetLogProbx () const
 
double MCMCGetLogProbx (unsigned ichain)
 
int MCMCGetPhase () const
 
const std::vector< double > & MCMCGetMaximumPoints () const
 
std::vector< double > MCMCGetMaximumPoint (unsigned i) const
 
const std::vector< double > & MCMCGetMaximumLogProb () const
 
int MCMCGetFlagInitialPosition () const
 
double MCMCGetRValueCriterion () const
 
double MCMCGetRValueParametersCriterion () const
 
double MCMCGetRValue () const
 
double MCMCGetRValueParameters (unsigned i)
 
bool MCMCGetRValueStrict () const
 
bool MCMCGetFlagRun () const
 
TTree * MCMCGetMarkovChainTree (unsigned i)
 
BCH1D * MCMCGetH1Marginalized (unsigned i)
 
BCH2D * MCMCGetH2Marginalized (unsigned i, unsigned j)
 
BCParameter * GetParameter (int index) const
 
BCParameter * GetParameter (const char *name) const
 
unsigned int GetNParameters () const
 
unsigned int GetNFixedParameters ()
 
unsigned int GetNFreeParameters ()
 
const std::vector< double > & GetBestFitParametersMarginalized () const
 
Setters
void MCMCSetTrialFunctionScaleFactor (std::vector< double > scale)
 
void MCMCSetNChains (unsigned n)
 
void MCMCSetNLag (unsigned n)
 
void MCMCSetNIterationsMax (unsigned n)
 
void MCMCSetNIterationsRun (unsigned n)
 
void MCMCSetNIterationsPreRunMin (unsigned n)
 
void MCMCSetNIterationsUpdate (unsigned n)
 
void MCMCSetNIterationsUpdateMax (unsigned n)
 
void MCMCSetMinimumEfficiency (double efficiency)
 
void MCMCSetMaximumEfficiency (double efficiency)
 
void MCMCSetRandomSeed (unsigned seed)
 
void MCMCSetWriteChainToFile (bool flag)
 
void MCMCSetWritePreRunToFile (bool flag)
 
void MCMCSetInitialPositions (const std::vector< double > &x0s)
 
void MCMCSetInitialPositions (std::vector< std::vector< double > > x0s)
 
void MCMCSetFlagInitialPosition (int flag)
 
void MCMCSetFlagOrderParameters (bool flag)
 
void MCMCSetFlagFillHistograms (bool flag)
 
void MCMCSetFlagPreRun (bool flag)
 
void MCMCSetRValueCriterion (double r)
 
void MCMCSetRValueParametersCriterion (double r)
 
void MCMCSetRValueStrict (bool strict=true)
 
void MCMCSetMarkovChainTrees (const std::vector< TTree * > &trees)
 
void MCMCInitializeMarkovChainTrees ()
 
int SetMarginalized (unsigned index, TH1D *h)
 
int SetMarginalized (unsigned index1, unsigned index2, TH2D *h)
 
void MCMCSetValuesDefault ()
 
void MCMCSetValuesQuick ()
 
void MCMCSetValuesDetail ()
 
void MCMCSetPrecision (BCEngineMCMC::Precision precision)
 
void SetNbins (unsigned int nbins)
 
Miscellaneous methods
void Copy (const BCEngineMCMC &enginemcmc)
 
virtual int AddParameter (const char *name, double min, double max, const char *latexname="")
 
virtual int AddParameter (BCParameter *parameter)
 
virtual void MCMCTrialFunction (unsigned ichain, std::vector< double > &x)
 
virtual double MCMCTrialFunctionSingle (unsigned ichain, unsigned ipar)
 
bool MCMCGetProposalPointMetropolis (unsigned chain, std::vector< double > &x)
 
bool MCMCGetProposalPointMetropolis (unsigned chain, unsigned parameter, std::vector< double > &x)
 
bool MCMCGetNewPointMetropolis (unsigned chain=0)
 
bool MCMCGetNewPointMetropolis (unsigned chain, unsigned parameter)
 
void MCMCInChainCheckMaximum ()
 
void MCMCInChainUpdateStatistics ()
 
void MCMCInChainFillHistograms ()
 
void MCMCInChainTestConvergenceAllChains ()
 
void MCMCInChainWriteChains ()
 
virtual double LogEval (const std::vector< double > &parameters)
 
int MCMCMetropolis ()
 
int MCMCMetropolisPreRun ()
 
void MCMCResetRunStatistics ()
 
void MCMCInitializeMarkovChains ()
 
int MCMCInitialize ()
 
virtual void ResetResults ()
 
void ClearParameters (bool hard=false)
 
virtual void MCMCIterationInterface ()
 
virtual void MCMCUserIterationInterface ()
 
virtual void MCMCCurrentPointInterface (std::vector< double > &, int, bool)
 

Protected Attributes

BCParameterSet fParameters
 
unsigned fMCMCNChains
 
unsigned fMCMCNLag
 
std::vector< unsigned > fMCMCNIterations
 
int fMCMCCurrentIteration
 
int fMCMCCurrentChain
 
unsigned fMCMCNIterationsUpdate
 
unsigned fMCMCNIterationsUpdateMax
 
int fMCMCNIterationsConvergenceGlobal
 
bool fMCMCFlagConvergenceGlobal
 
unsigned fMCMCNIterationsMax
 
unsigned fMCMCNIterationsRun
 
unsigned fMCMCNIterationsPreRunMin
 
std::vector< int > fMCMCNTrialsTrue
 
unsigned fMCMCNTrials
 
bool fMCMCFlagWriteChainToFile
 
bool fMCMCFlagWritePreRunToFile
 
std::vector< double > fMCMCTrialFunctionScaleFactor
 
std::vector< double > fMCMCTrialFunctionScaleFactorStart
 
bool fMCMCFlagPreRun
 
bool fMCMCFlagRun
 
std::vector< double > fMCMCInitialPosition
 
double fMCMCEfficiencyMin
 
double fMCMCEfficiencyMax
 
int fMCMCFlagInitialPosition
 
bool fMCMCFlagOrderParameters
 
int fMCMCPhase
 
std::vector< double > fMCMCx
 
std::vector< double > fMCMCxMax
 
std::vector< double > fMCMCxMean
 
std::vector< double > fMCMCxVar
 
std::vector< double > fMCMCprob
 
std::vector< double > fMCMCprobMax
 
std::vector< double > fMCMCprobMean
 
std::vector< double > fMCMCprobVar
 
bool fMCMCRValueUseStrict
 
double fMCMCRValueCriterion
 
double fMCMCRValueParametersCriterion
 
double fMCMCRValue
 
std::vector< double > fMCMCRValueParameters
 
TRandom3 * fRandom
 
std::vector< TH1D * > fMCMCH1Marginalized
 
std::vector< TH2D * > fMCMCH2Marginalized
 
std::vector< TTree * > fMCMCTrees
 
std::vector< double > fMCMCBestFitParameters
 
double fMCMCLogMaximum
 
std::vector< double > fMarginalModes
 

Private Types

typedef bool(BCEngineMCMC::* MCMCPointerToGetProposalPoint )(int chain, std::vector< double > xnew, std::vector< double > xold) const
 

Private Member Functions

void SyncThreadStorage ()
 

Private Attributes

MCMCPointerToGetProposalPoint fMCMCPointerToGetProposalPoint
 
std::vector
< MCMCThreadLocalStorage > 		fMCMCThreadLocalStorage
 
bool fFlagMarginalized
 

Detailed Description

An engine class for Markov Chain Monte Carlo.

Author
Daniel Kollar
Kevin Kröninger
Version
1.0
Date
08.2008 This class represents an engine class for Markov Chain Monte Carlo (MCMC). One or more chains can be defined simultaneously.

Definition at line 41 of file BCEngineMCMC.h.

Member Typedef Documentation

typedef bool(BCEngineMCMC::* BCEngineMCMC::MCMCPointerToGetProposalPoint)(int chain, std::vector< double > xnew, std::vector< double > xold) const
private

Defines a type of a pointer to a member function.

Definition at line 667 of file BCEngineMCMC.h.

Member Enumeration Documentation

enum BCEngineMCMC::Precision

An enumerator for the status of a test.

Enumerator
kLow 
kMedium 
kHigh 
kVeryHigh 

Definition at line 50 of file BCEngineMCMC.h.

50 { kLow, kMedium, kHigh, kVeryHigh };

Constructor & Destructor Documentation

BCEngineMCMC::BCEngineMCMC ( )

Default constructor.

Definition at line 28 of file BCEngineMCMC.cxx.

29 {
30  // set default parameters for the mcmc
31  MCMCSetValuesDefault();
32 
33  // initialize random number generator
34  fRandom = new TRandom3();
35  MCMCSetRandomSeed(0);
36 }
BCEngineMCMC::BCEngineMCMC ( const BCEngineMCMC &  enginemcmc)

Default copy constructor.

Definition at line 174 of file BCEngineMCMC.cxx.

175 {
176  Copy(other);
177 }
BCEngineMCMC::~BCEngineMCMC ( )
virtual

Default destructor.

Definition at line 157 of file BCEngineMCMC.cxx.

158 {
159  // delete random number generator
160  delete fRandom;
161 
162  // delete 1-d marginalized distributions
163  for (unsigned i = 0; i < fMCMCH1Marginalized.size(); ++i)
164  delete fMCMCH1Marginalized[i];
165  fMCMCH1Marginalized.clear();
166 
167  // delete 2-d marginalized distributions
168  for (unsigned i = 0; i < fMCMCH2Marginalized.size(); ++i)
169  delete fMCMCH2Marginalized[i];
170  fMCMCH2Marginalized.clear();
171 }

Member Function Documentation

int BCEngineMCMC::AddParameter ( const char *  name,
double  min,
double  max,
const char *  latexname = "" 
)
virtual

Adds a parameter.

Parameters
minminimum value of the parameter
maxmaximum value of the parameter
latexnameOptional latexname used for plotting
Returns
number of parameters after adding

Definition at line 1489 of file BCEngineMCMC.cxx.

1490 {
1491  // todo memory leak:
1492  // ==8243== 64 (44 direct, 20 indirect) bytes in 1 blocks are definitely lost in loss record 20,913 of 29,122
1493  // ==8243== at 0x402C9B4: operator new(unsigned int) (in /usr/lib/valgrind/vgpreload_memcheck-x86-linux.so)
1494  // ==8243== by 0x409DED1: BCEngineMCMC::AddParameter(char const*, double, double) (BCEngineMCMC.cxx:1480)
1495  return AddParameter(new BCParameter(name, min, max, latexname));
1496 }
int BCEngineMCMC::AddParameter ( BCParameter *  parameter)
virtual

Adds a parameter to the model.

Parameters
parameterA model parameter
See Also
AddParameter(const char * name, double lowerlimit, double upperlimit);

Reimplemented in BCModel.

Definition at line 1499 of file BCEngineMCMC.cxx.

1500 {
1501  return fParameters.Add(par);
1502 }
void BCEngineMCMC::ClearParameters ( bool  hard = false)
inline

Empty the sequence of parameters.

Parameters
hardDelete the parameters.
Warning
Could lead to problems if multiple modes use identical parameters.

Definition at line 632 of file BCEngineMCMC.h.

633  { fParameters.Clear(hard); }
void BCEngineMCMC::Copy ( const BCEngineMCMC &  enginemcmc)

Copy object

Parameters
enginemcmcObject to copy from

Definition at line 180 of file BCEngineMCMC.cxx.

181 {
182  fMCMCPointerToGetProposalPoint = other.fMCMCPointerToGetProposalPoint;
183  fMCMCNChains = other.fMCMCNChains;
184  fMCMCNLag = other.fMCMCNLag;
185  fMCMCNIterations = other.fMCMCNIterations;
186  fMCMCCurrentIteration = other.fMCMCCurrentIteration;
187  fMCMCCurrentChain = other.fMCMCCurrentChain;
188  fMCMCNIterationsUpdate = other.fMCMCNIterationsUpdate;
189  fMCMCNIterationsUpdateMax = other.fMCMCNIterationsUpdateMax;
190  fMCMCNIterationsConvergenceGlobal = other.fMCMCNIterationsConvergenceGlobal;
191  fMCMCFlagConvergenceGlobal = other.fMCMCFlagConvergenceGlobal;
192  fMCMCNIterationsMax = other.fMCMCNIterationsMax;
193  fMCMCNIterationsRun = other.fMCMCNIterationsRun;
194  fMCMCNIterationsPreRunMin = other.fMCMCNIterationsPreRunMin;
195  fMCMCNTrialsTrue = other.fMCMCNTrialsTrue;
196  fMCMCNTrials = other.fMCMCNTrials;
197  fMCMCFlagWriteChainToFile = other.fMCMCFlagWriteChainToFile;
198  fMCMCFlagWritePreRunToFile = other.fMCMCFlagWritePreRunToFile;
199  fMCMCTrialFunctionScaleFactor = other.fMCMCTrialFunctionScaleFactor;
200  fMCMCTrialFunctionScaleFactorStart = other.fMCMCTrialFunctionScaleFactorStart;
201  fMCMCFlagPreRun = other.fMCMCFlagPreRun;
202  fMCMCFlagRun = other.fMCMCFlagRun;
203  fMCMCInitialPosition = other.fMCMCInitialPosition;
204  fMCMCEfficiencyMin = other.fMCMCEfficiencyMin;
205  fMCMCEfficiencyMax = other.fMCMCEfficiencyMax;
206  fMCMCFlagInitialPosition = other.fMCMCFlagInitialPosition;
207  fMCMCFlagOrderParameters = other.fMCMCFlagOrderParameters;
208  fMCMCPhase = other.fMCMCPhase;
209  fMCMCx = other.fMCMCx;
210  fMCMCxMax = other.fMCMCxMax;
211  fMCMCxMean = other.fMCMCxMean;
212  fMCMCxVar = other.fMCMCxVar;
213  fMCMCprob = other.fMCMCprob;
214  fMCMCprobMax = other.fMCMCprobMax;
215  fMCMCprobMean = other.fMCMCprobMean;
216  fMCMCprobVar = other.fMCMCprobVar;
217  fMCMCRValueUseStrict = other.fMCMCRValueUseStrict;
218  fMCMCRValueCriterion = other.fMCMCRValueCriterion ;
219  fMCMCRValueParametersCriterion = other.fMCMCRValueParametersCriterion;
220  fMCMCRValue = other.fMCMCRValue;
221  fMCMCRValueParameters = other.fMCMCRValueParameters;
222  if (other.fRandom)
223  {
224  fRandom = new TRandom3(*other.fRandom);
225  }
226  else
227  {
228  fRandom = NULL;
229  }
230 
231  fMCMCThreadLocalStorage = other.fMCMCThreadLocalStorage;
232 
233  for (unsigned i = 0; i < other.fMCMCH1Marginalized.size(); ++i) {
234  if (other.fMCMCH1Marginalized.at(i))
235  fMCMCH1Marginalized.push_back(new TH1D(*(other.fMCMCH1Marginalized.at(i))));
236  else
237  fMCMCH1Marginalized.push_back(0);
238  }
239 
240  for (unsigned i = 0; i < other.fMCMCH2Marginalized.size(); ++i) {
241  if (other.fMCMCH2Marginalized.at(i))
242  fMCMCH2Marginalized.push_back(new TH2D(*(other.fMCMCH2Marginalized.at(i))));
243  else
244  fMCMCH2Marginalized.push_back(0);
245  }
246 
247  for (unsigned i = 0; i < other.fMCMCTrees.size(); ++i) {
248  fMCMCTrees.push_back(0);
249  }
250 
251  fMarginalModes = other.fMarginalModes;
252  fMCMCBestFitParameters = other.fMCMCBestFitParameters;
253  fMCMCLogMaximum = other.fMCMCLogMaximum;
254 }
const std::vector< double > & BCEngineMCMC::GetBestFitParametersMarginalized ( ) const

Returns the set of values of the parameters at the modes of the marginalized posterior pdfs.

Returns
best fit parameters

Definition at line 352 of file BCEngineMCMC.cxx.

353 {
354  if(fMarginalModes.empty())
355  BCLog::OutError("BCIntegrate::GetBestFitParameterMarginalized : MCMC not yet run, returning center of the range.");
356 
357  return fMarginalModes;
358 }
unsigned int BCEngineMCMC::GetNFixedParameters ( )
Returns
The number of fixed parameters.

Definition at line 270 of file BCEngineMCMC.cxx.

271 {
272  int n = 0;
273  for (unsigned int i = 0; i < fParameters.Size(); ++i) {
274  if (fParameters[i]->Fixed())
275  ++n;
276  }
277 
278  return n;
279 }
unsigned int BCEngineMCMC::GetNFreeParameters ( )
Returns
The number of free parameters.

Definition at line 264 of file BCEngineMCMC.cxx.

265 {
266  return (GetNParameters() - GetNFixedParameters());
267 }
unsigned int BCEngineMCMC::GetNParameters ( ) const
inline
Returns
The number of parameters of the model.

Definition at line 297 of file BCEngineMCMC.h.

298  { return fParameters.Size(); }
BCParameter* BCEngineMCMC::GetParameter ( int  index) const
inline
Parameters
indexThe index of the parameter in the parameter set.
Returns
The parameter.

Definition at line 286 of file BCEngineMCMC.h.

287  { return fParameters.Get(index); }
BCParameter* BCEngineMCMC::GetParameter ( const char *  name) const
inline
Parameters
nameThe name of the parameter in the parameter set.
Returns
The parameter.

Definition at line 292 of file BCEngineMCMC.h.

293  { return fParameters.Get(name); }
double BCEngineMCMC::LogEval ( const std::vector< double > &  parameters)
virtual

Needs to be overloaded in the derived class.

Returns
natural logarithm of the function to map with MCMC

Reimplemented in BCIntegrate, and BCModel.

Definition at line 880 of file BCEngineMCMC.cxx.

881 {
882  // test function for now
883  // this will be overloaded by the user
884  return 0.0;
885 }
virtual void BCEngineMCMC::MCMCCurrentPointInterface ( std::vector< double > &  ,
int  ,
bool   
)
inlinevirtual

Interface allowing to execute arbitrary code for each new point of the MCMC. This method needs to be overloaded in the derived class

Parameters
pointpoint that was generated and checked
ichainindex of the chain
acceptedflag whether or not the point was accepted for the chain

Definition at line 657 of file BCEngineMCMC.h.

658  {}
unsigned BCEngineMCMC::MCMCGetCurrentChain ( ) const
inline
Returns
current chain index

Definition at line 102 of file BCEngineMCMC.h.

103  { return fMCMCCurrentChain; }
unsigned BCEngineMCMC::MCMCGetCurrentIteration ( ) const
inline
Returns
current iterations

Definition at line 97 of file BCEngineMCMC.h.

98  { return fMCMCCurrentIteration; }
bool BCEngineMCMC::MCMCGetFlagConvergenceGlobal ( ) const
inline
Returns
flag if converged or not

Definition at line 113 of file BCEngineMCMC.h.

114  { return fMCMCFlagConvergenceGlobal; }
int BCEngineMCMC::MCMCGetFlagInitialPosition ( ) const
inline
Returns
flag which defined initial position

Definition at line 227 of file BCEngineMCMC.h.

228  { return fMCMCFlagInitialPosition; }
bool BCEngineMCMC::MCMCGetFlagRun ( ) const
inline
Returns
the flag if MCMC has been performed or not

Definition at line 257 of file BCEngineMCMC.h.

258  { return fMCMCFlagRun; }
BCH1D * BCEngineMCMC::MCMCGetH1Marginalized ( unsigned  i)

Retrieve a histogram of the 1D marginalized distribution of a single parameter.

Parameters
iindex of the parameter
Returns
pointer to the histogram

Definition at line 289 of file BCEngineMCMC.cxx.

290 {
291  if ( !fParameters.ValidIndex(index)) {
292  BCLog::OutError(Form("BCEngineMCMC::MCMCGetH1Marginalized. Index %u out of range.", index));
293  return 0;
294  }
295 
296  // use when BCIntegrate acts MCMC marginalization and only some marginals have been computed
297  if (!fMCMCH1Marginalized[index]) {
298  BCLog::OutWarning(Form("BCEngineMCMC::MCMCGetH1Marginalized: marginal distribution not computed/stored for par. %d", index));
299  return 0;
300  }
301 
302  // set histogram
303  BCH1D * hprob = new BCH1D();
304  hprob->SetHistogram(fMCMCH1Marginalized[index]);
305 
306  if (fMarginalModes.empty())
307  fMarginalModes.assign(fParameters.Size(), 0.0);
308  fMarginalModes[index] = hprob->GetMode();
309 
310  return hprob;
311 }
BCH2D * BCEngineMCMC::MCMCGetH2Marginalized ( unsigned  i,
unsigned  j 
)

Retrieve a histogram of the 2D marginalized distribution for two parameters.

Parameters
iindex of the first parameter
jindex of the second parameter
Returns
pointer to the histogram

Definition at line 314 of file BCEngineMCMC.cxx.

315 {
316  if ( !fParameters.ValidIndex(i)) {
317  BCLog::OutError(Form("BCEngineMCMC::MCMCGetH2Marginalized. Index %u out of range.", i));
318  return 0;
319  }
320  if ( !fParameters.ValidIndex(j)) {
321  BCLog::OutError(Form("BCEngineMCMC::MCMCGetH2Marginalized. Index %u out of range.", j));
322  return 0;
323  }
324  if (i == j) {
325  BCLog::OutError(Form("BCEngineMCMC::MCMCGetH2Marginalized. Called with identical indices %u.", i));
326  return 0;
327  }
328 
329  // swap indices
330  if (i > j) {
331  unsigned indexTemp = i;
332  i = j;
333  j = indexTemp;
334  }
335 
336  // memory layout for n parameters and indices i, j:
337  // first (n-1) elements for first parameter vs. all other parameters
338  // then (n-2) elements for second parameter and all others etc
339  // so the first combination for which i is the lower index is at n*i - i(i+1)/2
340  // and the offset is given by (j-i-1)
341  TH2D * h = fMCMCH2Marginalized.at(GetNParameters() * i - (i * i + 3 * i) / 2 + j - 1);
342  if ( !h)
343  return 0;
344 
345  BCH2D * hprob = new BCH2D();
346  hprob->SetHistogram(h);
347 
348  return hprob;
349 }
const std::vector<double>& BCEngineMCMC::MCMCGetLogProbx ( ) const
inline
Returns
log of the probability of the current points of each Markov chain

Definition at line 197 of file BCEngineMCMC.h.

198  { return fMCMCprob; }
double BCEngineMCMC::MCMCGetLogProbx ( unsigned  ichain)
Returns
log of the probability of the current points of the Markov chain.
Parameters
ichainchain index

Definition at line 568 of file BCEngineMCMC.cxx.

569 {
570  // check if ichain is in range
571  if (ichain >= fMCMCNChains)
572  return -1;
573 
574  // return log of the probability at the current point in the ith chain
575  return fMCMCprob.at(ichain);
576 }
TTree* BCEngineMCMC::MCMCGetMarkovChainTree ( unsigned  i)
inline

Retrieve the tree containing the Markov chain.

Parameters
iindex of the Markov chain
Returns
pointer to the tree

Definition at line 264 of file BCEngineMCMC.h.

265  { if (i < fMCMCTrees.size())
266  return fMCMCTrees.at(i);
267  else
268  return 0; }
const std::vector<double>& BCEngineMCMC::MCMCGetMaximumLogProb ( ) const
inline
Returns
maximum (log) probability of each Markov chain

Definition at line 222 of file BCEngineMCMC.h.

223  { return fMCMCprobMax; }
std::vector< double > BCEngineMCMC::MCMCGetMaximumPoint ( unsigned  i) const
Returns
maximum point of Markov chain
Parameters
iThe index of the Markov chain

Definition at line 361 of file BCEngineMCMC.cxx.

362 {
363  // create a new vector with the length of fMCMCNParameters
364  std::vector<double> x;
365 
366  // check if i is in range
367  if (i >= fMCMCNChains)
368  return x;
369 
370  // copy the point in the ith chain into the temporary vector
371  for (unsigned j = 0; j < fParameters.Size(); ++j)
372  x.push_back(fMCMCxMax.at(i * fParameters.Size() + j));
373 
374  return x;
375 }
const std::vector<double>& BCEngineMCMC::MCMCGetMaximumPoints ( ) const
inline
Returns
maximum points of each Markov chain

Definition at line 212 of file BCEngineMCMC.h.

213  { return fMCMCxMax; }
unsigned BCEngineMCMC::MCMCGetNChains ( ) const
inline
Returns
number of Markov chains

Definition at line 82 of file BCEngineMCMC.h.

83  { return fMCMCNChains; }
bool BCEngineMCMC::MCMCGetNewPointMetropolis ( unsigned  chain = 0)

Generates a new point using the Metropolis algorithm.

Parameters
chainchain index

Definition at line 691 of file BCEngineMCMC.cxx.

692 {
693  // calculate index
694  unsigned index = chain * fParameters.Size();
695 
696  fMCMCCurrentChain = chain;
697 
698  // increase counter
699  fMCMCNIterations[chain]++;
700 
701  // get proposal point
702  if (!MCMCGetProposalPointMetropolis(chain, fMCMCThreadLocalStorage[chain].xLocal))
703  {
704  // execute user code for every point
705  MCMCCurrentPointInterface(fMCMCThreadLocalStorage[chain].xLocal, chain, false);
706 
707  return false;
708  }
709 
710  // calculate probabilities of the old and new points
711  double p0 = fMCMCprob[chain];
712  double p1 = LogEval(fMCMCThreadLocalStorage[chain].xLocal);
713 
714  // flag for accept
715  bool accept = false;
716 
717  // if the new point is more probable, keep it ...
718  if (p1 >= p0)
719  accept = true;
720 
721  // ... or else throw dice.
722  else
723  {
724  double r = log(fMCMCThreadLocalStorage[chain].rng->Rndm());
725 
726  if(r < p1 - p0)
727  accept = true;
728  }
729 
730  // fill the new point
731  if(accept)
732  {
733  // increase counter
734  for (unsigned i = 0; i < fParameters.Size(); ++i)
735  fMCMCNTrialsTrue[chain * fParameters.Size() + i]++;
736 
737  // copy the point
738  for(unsigned i = 0; i < fParameters.Size(); ++i)
739  {
740  // save the point
741  fMCMCx[index + i] = fMCMCThreadLocalStorage[chain].xLocal[i];
742 
743  // save the probability of the point
744  fMCMCprob[chain] = p1;
745  }
746  }
747 
748  // execute user code for every point
749  MCMCCurrentPointInterface(fMCMCThreadLocalStorage[chain].xLocal, chain, accept);
750 
751  return accept;
752 }
bool BCEngineMCMC::MCMCGetNewPointMetropolis ( unsigned  chain,
unsigned  parameter 
)

Definition at line 629 of file BCEngineMCMC.cxx.

630 {
631  // calculate index
632  unsigned index = chain * fParameters.Size();
633 
634  fMCMCCurrentChain = chain;
635 
636  // increase counter
637  fMCMCNIterations[chain]++;
638 
639  // get proposal point
640  if (!MCMCGetProposalPointMetropolis(chain, parameter, fMCMCThreadLocalStorage[chain].xLocal))
641  {
642  // execute user code for every point
643  MCMCCurrentPointInterface(fMCMCThreadLocalStorage[chain].xLocal, chain, false);
644 
645  return false;
646  }
647 
648  // calculate probabilities of the old and new points
649  double p0 = fMCMCprob[chain];
650  double p1 = LogEval(fMCMCThreadLocalStorage[chain].xLocal);
651 
652  // flag for accept
653  bool accept = false;
654 
655  // if the new point is more probable, keep it ...
656  if (p1 >= p0)
657  accept = true;
658  // ... or else throw dice.
659  else
660  {
661  double r = log(fMCMCThreadLocalStorage[chain].rng->Rndm());
662 
663  if(r < p1 - p0)
664  accept = true;
665  }
666 
667  // fill the new point
668  if(accept)
669  {
670  // increase counter
671  fMCMCNTrialsTrue[chain * fParameters.Size() + parameter]++;
672 
673  // copy the point
674  for(unsigned i = 0; i < fParameters.Size(); ++i)
675  {
676  // save the point
677  fMCMCx[index + i] = fMCMCThreadLocalStorage[chain].xLocal[i];
678 
679  // save the probability of the point
680  fMCMCprob[chain] = p1;
681  }
682  }
683 
684  // execute user code for every point
685  MCMCCurrentPointInterface(fMCMCThreadLocalStorage[chain].xLocal, chain, accept);
686 
687  return accept;
688 }
const std::vector<unsigned>& BCEngineMCMC::MCMCGetNIterations ( ) const
inline
Returns
number of iterations

Definition at line 92 of file BCEngineMCMC.h.

93  { return fMCMCNIterations; }
unsigned BCEngineMCMC::MCMCGetNIterationsConvergenceGlobal ( ) const
inline
Returns
number of iterations needed for all chains to converge simultaneously

Definition at line 108 of file BCEngineMCMC.h.

109  { return fMCMCNIterationsConvergenceGlobal; }
unsigned BCEngineMCMC::MCMCGetNIterationsMax ( ) const
inline
Returns
maximum number of iterations for a Markov chain

Definition at line 118 of file BCEngineMCMC.h.

119  { return fMCMCNIterationsMax; }
unsigned BCEngineMCMC::MCMCGetNIterationsPreRunMin ( ) const
inline
Returns
minimum number of pre-run iterations for a Markov chain

Definition at line 128 of file BCEngineMCMC.h.

129  { return fMCMCNIterationsPreRunMin; }
unsigned BCEngineMCMC::MCMCGetNIterationsRun ( ) const
inline
Returns
number of iterations for a Markov chain

Definition at line 123 of file BCEngineMCMC.h.

124  { return fMCMCNIterationsRun; }
unsigned BCEngineMCMC::MCMCGetNIterationsUpdate ( ) const
inline
Returns
number of iterations after statistics update.

Definition at line 133 of file BCEngineMCMC.h.

134  { return fMCMCNIterationsUpdate; }
unsigned BCEngineMCMC::MCMCGetNIterationsUpdateMax ( ) const
inline
Returns
maximum number of iterations after statistics update.

Definition at line 138 of file BCEngineMCMC.h.

139  { return fMCMCNIterationsUpdateMax; }
unsigned BCEngineMCMC::MCMCGetNLag ( ) const
inline
Returns
lag of the Markov chains

Definition at line 87 of file BCEngineMCMC.h.

88  { return fMCMCNLag; }
int BCEngineMCMC::MCMCGetNTrials ( ) const
inline
Returns
number of trials

Definition at line 148 of file BCEngineMCMC.h.

149  { return fMCMCNTrials; }
std::vector<int> BCEngineMCMC::MCMCGetNTrialsTrue ( ) const
inline
Returns
number of accepted trials for each chain

Definition at line 143 of file BCEngineMCMC.h.

144  { return fMCMCNTrialsTrue; }
int BCEngineMCMC::MCMCGetPhase ( ) const
inline
Returns
pointer to the phase of a run.

Definition at line 207 of file BCEngineMCMC.h.

208  { return fMCMCPhase; }
const std::vector<double>& BCEngineMCMC::MCMCGetprobMean ( ) const
inline
Returns
mean value of the probability for each chain up to the current iteration

Definition at line 154 of file BCEngineMCMC.h.

155  { return fMCMCprobMean; }
bool BCEngineMCMC::MCMCGetProposalPointMetropolis ( unsigned  chain,
std::vector< double > &  x 
)

Returns a trial point for the Metropolis algorithm.

Parameters
chainchain index
xproposal point
Returns
flag indicating whether the new point lies within the allowed range

Definition at line 579 of file BCEngineMCMC.cxx.

580 {
581  // get unscaled random point. this point might not be in the correct volume.
582  MCMCTrialFunction(chain, x);
583 
584  // get a proposal point from the trial function and scale it
585  for (unsigned i = 0; i < fParameters.Size(); ++i) {
586  // check if parameter is fixed
587  if (fParameters[i]->Fixed()) {
588  x[i] = 0;
589  }
590  x[i] = fMCMCx[chain * fParameters.Size() + i] + x[i] * fParameters[i]->GetRangeWidth();
591  }
592 
593  // check if the point is in the correct volume.
594  for (unsigned i = 0; i < fParameters.Size(); ++i)
595  if (!fParameters[i]->IsValid(x[i]))
596  return false;
597 
598  return true;
599 }
bool BCEngineMCMC::MCMCGetProposalPointMetropolis ( unsigned  chain,
unsigned  parameter,
std::vector< double > &  x 
)

Returns a trial point for the Metropolis algorithm.

Parameters
chainchain index
xproposal point
Returns
flag indicating whether the new point lies within the allowed range

Definition at line 602 of file BCEngineMCMC.cxx.

603 {
604  // copy the old point into the new
605  for (unsigned i = 0; i < fParameters.Size(); ++i)
606  x[i] = fMCMCx[ichain * fParameters.Size() + i];
607 
608  // check if parameter is fixed
609  if (fParameters[ipar]->Fixed()) {
610  x[ipar] = fParameters[ipar]->GetFixedValue();
611  return true; // assume that value is inside allowed region
612  }
613 
614  // get unscaled random point in the dimension of the chosen
615  // parameter. this point might not be in the correct volume.
616  double proposal = MCMCTrialFunctionSingle(ichain, ipar);
617 
618  // modify the parameter under study
619  x[ipar] += proposal * fParameters[ipar]->GetRangeWidth();
620 
621  // check if the point is in the correct volume.
622  if (fParameters[ipar]->IsValid(x[ipar]))
623  return true;
624 
625  return false;
626 }
double BCEngineMCMC::MCMCGetRValue ( ) const
inline
Returns
R-value

Definition at line 242 of file BCEngineMCMC.h.

243  { return fMCMCRValue; }
double BCEngineMCMC::MCMCGetRValueCriterion ( ) const
inline
Returns
R-value criterion

Definition at line 232 of file BCEngineMCMC.h.

233  { return fMCMCRValueCriterion; }
double BCEngineMCMC::MCMCGetRValueParameters ( unsigned  i)
inline
Returns
R-value for a parameter
Parameters
iparameter index

Definition at line 248 of file BCEngineMCMC.h.

249  { return fMCMCRValueParameters.at(i); }
double BCEngineMCMC::MCMCGetRValueParametersCriterion ( ) const
inline
Returns
R-value criterion for parameters

Definition at line 237 of file BCEngineMCMC.h.

238  { return fMCMCRValueParametersCriterion; }
bool BCEngineMCMC::MCMCGetRValueStrict ( ) const
inline

Use strict or relaxed rule for Gelman/Rubin R-value

Definition at line 252 of file BCEngineMCMC.h.

253  { return fMCMCRValueUseStrict; }
const std::vector<double>& BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( ) const
inline
Returns
scale factor for all parameters and chains

Definition at line 165 of file BCEngineMCMC.h.

166  { return fMCMCTrialFunctionScaleFactor; }
std::vector< double > BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( unsigned  ichain) const
Returns
scale factor for all parameters and achain.
Parameters
ichainchain index

Definition at line 504 of file BCEngineMCMC.cxx.

505 {
506  // create a new vector with the length of fParameters
507  std::vector<double> x;
508 
509  // check if ichain is in range
510  if (ichain >= fMCMCNChains)
511  return x;
512 
513  // copy the scale factors into the temporary vector
514  for (unsigned j = 0; j < fParameters.Size(); ++j)
515  x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fParameters.Size() + j));
516 
517  return x;
518 }
double BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( unsigned  ichain,
unsigned  ipar 
)
Returns
scale factor for a parameter and a chain.
Parameters
ichainchain index
iparparameter index

Definition at line 521 of file BCEngineMCMC.cxx.

522 {
523  // check if ichain is in range
524  if (ichain >= fMCMCNChains)
525  return 0;
526 
527  // check if ipar is in range
528  if (ipar >= fParameters.Size())
529  return 0;
530 
531  // return component of ipar point in the ichain chain
532  return fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNChains + ipar);
533 }
const std::vector<double>& BCEngineMCMC::MCMCGetVariance ( ) const
inline
Returns
mean value of the probability for each chain up to the current iteration

Definition at line 160 of file BCEngineMCMC.h.

161  { return fMCMCprobVar; }
const std::vector<double>& BCEngineMCMC::MCMCGetx ( ) const
inline
Returns
current point of each Markov chain

Definition at line 181 of file BCEngineMCMC.h.

182  { return fMCMCx; }
std::vector< double > BCEngineMCMC::MCMCGetx ( unsigned  ichain)
Parameters
ichainindex of the Markov chain
Returns
current point of the Markov chain

Definition at line 536 of file BCEngineMCMC.cxx.

537 {
538  // create a new vector with the length of fParameters.Size()
539  std::vector<double> x;
540 
541  // check if ichain is in range
542  if (ichain >= fMCMCNChains)
543  return x;
544 
545  // copy the point in the ichain chain into the temporary vector
546  for (unsigned j = 0; j < fParameters.Size(); ++j)
547  x.push_back(fMCMCx.at(ichain * fParameters.Size() + j));
548 
549  return x;
550 }
double BCEngineMCMC::MCMCGetx ( unsigned  ichain,
unsigned  ipar 
) const
Parameters
ichainchain index
iparparameter index
Returns
parameter of the Markov chain

Definition at line 553 of file BCEngineMCMC.cxx.

554 {
555  // check if ichain is in range
556  if (ichain >= fMCMCNChains)
557  return 0;
558 
559  // check if ipar is in range
560  if (ipar >= fParameters.Size())
561  return 0;
562 
563  // return component of jth point in the ith chain
564  return fMCMCx.at(ichain * fParameters.Size() + ipar);
565 }
void BCEngineMCMC::MCMCInChainCheckMaximum ( )

Updates statistics: find new maximum

Definition at line 755 of file BCEngineMCMC.cxx.

756 {
757  // loop over all chains
758  for (unsigned i = 0; i < fMCMCNChains; ++i)
759  {
760  // check if new maximum is found or chain is at the beginning
761  if (fMCMCprob[i] > fMCMCprobMax[i] || fMCMCNIterations[i] == 1)
762  {
763  // copy maximum value
764  fMCMCprobMax[i] = fMCMCprob[i];
765 
766  // copy mode of chain
767  for (unsigned j = 0; j < fParameters.Size(); ++j)
768  fMCMCxMax[i * fParameters.Size() + j] = fMCMCx[i * fParameters.Size() + j];
769  }
770  }
771 }
void BCEngineMCMC::MCMCInChainFillHistograms ( )

Updates statistics: fill marginalized distributions

Definition at line 804 of file BCEngineMCMC.cxx.

805 {
806  // loop over chains
807  for (unsigned i = 0; i < fMCMCNChains; ++i)
808  {
809  // fill each 1-dimensional histogram (if supposed to be filled)
810  for (unsigned j = 0; j < fParameters.Size(); ++j)
811  if (TH1 * h = fMCMCH1Marginalized[j])
812  h->Fill(fMCMCx[i * fParameters.Size() + j]);
813 
814  // fill each 2-dimensional histogram (if supposed to be filled)
815  unsigned counter = 0;
816 
817  for (unsigned j = 0; j < fParameters.Size(); ++j)
818  for (unsigned k = j+1; k < fParameters.Size(); ++k)
819  {
820  if (TH2D * h = fMCMCH2Marginalized[counter])
821  h->Fill(fMCMCx[i*fParameters.Size()+j],fMCMCx[i* fParameters.Size()+k]);
822  counter ++;
823  }
824  }
825 }
void BCEngineMCMC::MCMCInChainTestConvergenceAllChains ( )

Updates statistics: check convergence

Definition at line 828 of file BCEngineMCMC.cxx.

829 {
830  if (fMCMCNChains > 1 && fMCMCNTrials > 1)
831  {
832  // define flag for convergence
833  bool flag_convergence = true;
834 
835  // extract means and variances
836  std::vector<double> means(fMCMCNChains);
837  std::vector<double> variances(fMCMCNChains);
838 
839  // loop over parameters
840  for (unsigned iparameters = 0; iparameters < fParameters.Size(); ++iparameters){
841  if (fParameters[iparameters]->Fixed())
842  continue;
843 
844  // loop over chains
845  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains) {
846  // get parameter index
847  unsigned index = ichains * fParameters.Size() + iparameters;
848  means[ichains] = fMCMCxMean[index];
849  variances[ichains] = fMCMCxVar[index];
850  }
851  fMCMCRValueParameters[iparameters] = BCMath::Rvalue(means, variances, fMCMCNTrials, fMCMCRValueUseStrict);
852 
853  // set flag to false if convergence criterion is not fulfilled for the parameter
854  if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
855  flag_convergence = false;
856 
857  // else: leave convergence flag true for that parameter
858  }
859 
860  fMCMCRValue = BCMath::Rvalue(fMCMCprobMean, fMCMCprobVar, fMCMCNTrials, fMCMCRValueUseStrict);
861 
862  // set flag to false if convergence criterion is not fulfilled for the log-likelihood
863  if (!((fMCMCRValue - 1.0) < fMCMCRValueCriterion))
864  flag_convergence = false;
865 
866  // remember number of iterations needed to converge
867  if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
868  fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / GetNFreeParameters();
869  }
870 }
void BCEngineMCMC::MCMCInChainUpdateStatistics ( )

Updates statistics:

Definition at line 774 of file BCEngineMCMC.cxx.

775 {
776  // length of vectors
777  unsigned nentries = fParameters.Size() * fMCMCNChains;
778 
779  // loop over all parameters of all chains
780  for (unsigned i = 0; i < nentries; ++i) {
781  // calculate mean value of each parameter in the chain for this part
782  fMCMCxMean[i] += (fMCMCx[i] - fMCMCxMean[i]) / double(fMCMCNTrials);
783 
784  // calculate variance of each chain for this part
785  if (fMCMCNTrials > 1)
786  fMCMCxVar[i] = (1.0 - 1./double(fMCMCNTrials)) * fMCMCxVar[i]
787  + (fMCMCx[i] - fMCMCxMean[i]) * (fMCMCx[i] - fMCMCxMean[i]) / double(fMCMCNTrials - 1);
788  }
789 
790  // loop over chains
791  for (unsigned i = 0; i < fMCMCNChains; ++i) {
792  // calculate mean value of each chain for this part
793  fMCMCprobMean[i] += (fMCMCprob[i] - fMCMCprobMean[i]) / double(fMCMCNTrials);
794 
795  // calculate variance of each chain for this part
796  if (fMCMCNTrials > 1)
797  fMCMCprobVar[i] = (1.0 - 1/double(fMCMCNTrials)) * fMCMCprobVar[i]
798  + (fMCMCprob[i] - fMCMCprobMean[i]) * (fMCMCprob[i] - fMCMCprobMean[i]) / double(fMCMCNTrials - 1);
799  }
800 
801 }
void BCEngineMCMC::MCMCInChainWriteChains ( )

Updates statistics: write chains to file

Definition at line 872 of file BCEngineMCMC.cxx.

873 {
874  // loop over all chains
875  for (unsigned i = 0; i < fMCMCNChains; ++i)
876  fMCMCTrees[i]->Fill();
877 }
int BCEngineMCMC::MCMCInitialize ( )

Initializes the engine.

Returns
An error code

Definition at line 1563 of file BCEngineMCMC.cxx.

1564 {
1565  // resource allocation must be done only by one thread
1566  // reset values
1567  ResetResults();
1568 
1569  // free memory for vectors
1570  fMCMCNIterations.assign(fMCMCNChains, 0);
1571  fMCMCprobMean.assign(fMCMCNChains, 0);
1572  fMCMCprobVar.assign(fMCMCNChains, 0);
1573  fMCMCprob.assign(fMCMCNChains, -1.0);
1574  fMCMCprobMax.assign(fMCMCNChains, -1.0);
1575 
1576  fMCMCNTrialsTrue.assign(fMCMCNChains * fParameters.Size(), 0);
1577  fMCMCNTrials = 0;
1578  fMCMCxMax.assign(fMCMCNChains * fParameters.Size(), 0.);
1579  fMCMCxMean.assign(fMCMCNChains * fParameters.Size(), 0);
1580  fMCMCxVar.assign(fMCMCNChains * fParameters.Size(), 0);
1581 
1582  fMCMCRValueParameters.assign(fParameters.Size(), 0.);
1583 
1584  SyncThreadStorage();
1585 
1586  if (fMCMCTrialFunctionScaleFactorStart.size() == 0)
1587  fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fParameters.Size(), 1.0);
1588  else
1589  for (unsigned i = 0; i < fMCMCNChains; ++i)
1590  for (unsigned j = 0; j < fParameters.Size(); ++j)
1591  fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j));
1592 
1593  // set initial position
1594  if (fMCMCFlagInitialPosition == 2) // user defined points
1595  {
1596  // define flag
1597  bool flag = true;
1598 
1599  // check the length of the array of initial positions
1600  if (fMCMCInitialPosition.size() != (fMCMCNChains * fParameters.Size()))
1601  {
1602  BCLog::OutError("BCEngine::MCMCInitialize : Length of vector containing initial positions does not have required length.");
1603  flag = false;
1604  }
1605 
1606  // check the boundaries
1607  if (flag)
1608  {
1609  for (unsigned j = 0; j < fMCMCNChains; ++j)
1610  for (unsigned i = 0; i < fParameters.Size(); ++i)
1611  if (fParameters[i]->IsValid(fMCMCInitialPosition[j * fParameters.Size() + i]))
1612  {
1613  BCLog::OutError("BCEngine::MCMCInitialize : Initial position out of boundaries.");
1614  flag = false;
1615  }
1616  }
1617 
1618  // check flag
1619  if (!flag)
1620  fMCMCFlagInitialPosition = 1;
1621  }
1622 
1623  if (fMCMCFlagInitialPosition == 0) // center of the interval
1624  for (unsigned j = 0; j < fMCMCNChains; ++j)
1625  for (unsigned i = 0; i < fParameters.Size(); ++i) {
1626  if (fParameters[i]->Fixed())
1627  fMCMCx.push_back(fParameters[i]->GetFixedValue());
1628  else
1629  fMCMCx.push_back(fParameters[i]->GetLowerLimit() + .5 * fParameters[i]->GetRangeWidth());
1630  }
1631 
1632  else if (fMCMCFlagInitialPosition == 2) // user defined
1633  {
1634  for (unsigned j = 0; j < fMCMCNChains; ++j)
1635  for (unsigned i = 0; i < fParameters.Size(); ++i) {
1636  if (fParameters[i]->Fixed()) {
1637  fMCMCx.push_back(fParameters[i]->GetFixedValue());
1638  BCLog::OutWarning("BCEngineMCMC::MCMCInitialize. Inconsisten start value. Changed parameter value to fixed value.");
1639  }
1640  else
1641  fMCMCx.push_back(fMCMCInitialPosition.at(j * fParameters.Size() + i));
1642  }
1643  }
1644 
1645  else
1646  {
1647  for (unsigned j = 0; j < fMCMCNChains; ++j) // random number (default)
1648  for (unsigned i = 0; i < fParameters.Size(); ++i) {
1649  if (fParameters[i]->Fixed())
1650  fMCMCx.push_back(fParameters[i]->GetFixedValue());
1651  else
1652  fMCMCx.push_back(fParameters[i]->GetLowerLimit() + fMCMCThreadLocalStorage[j].rng->Rndm() * fParameters[i]->GetRangeWidth());
1653  }
1654  }
1655 
1656  // copy the point of the first chain
1657  std::copy(fMCMCx.begin(), fMCMCx.begin() + fParameters.Size(), fMCMCThreadLocalStorage.at(0).xLocal.begin());
1658 
1659  // define 1-dimensional histograms for marginalization
1660  bool fillAny=false;
1661  for(unsigned i = 0; i < fParameters.Size(); ++i)
1662  {
1663  const BCParameter * p = fParameters[i];
1664  TH1D * h1 = NULL;
1665  if (p->FillHistograms() && ! p->Fixed()) {
1666  h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i),
1667  TString::Format(";%s;", p->GetLatexName().c_str()),
1668  p->GetNbins(), p->GetLowerLimit(), p->GetUpperLimit());
1669  h1->SetStats(kFALSE);
1670 
1671  fillAny = true;
1672  }
1673  fMCMCH1Marginalized.push_back(h1);
1674  }
1675  // if filling no histograms, set H1 vector to zero size, implies no 2D histograms either
1676  if (!fillAny) {
1677  fMCMCH1Marginalized.clear();
1678  }
1679  else {
1680  // define 2-dimensional histograms for marginalization
1681  for(unsigned i = 0; i < fParameters.Size(); ++i) {
1682  BCParameter * p1 = fParameters[i];
1683  for (unsigned j = i + 1; j < fParameters.Size(); ++j) {
1684  TH2D * h2 = 0;
1685  BCParameter * p2 = fParameters[j];
1686  if (p2->FillHistograms() && p1->FillHistograms() && ! p2->Fixed() && ! p1->Fixed()) {
1687  h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, j), "",
1688  p1->GetNbins(), p1->GetLowerLimit(), p1->GetUpperLimit(),
1689  p2->GetNbins(), p2->GetLowerLimit(), p2->GetUpperLimit());
1690 
1691  // decorate histogram
1692  h2->SetXTitle(Form("%s", p1->GetLatexName().data()));
1693  h2->SetYTitle(Form("%s", p2->GetLatexName().data()));
1694  h2->SetStats(kFALSE);
1695  }
1696  fMCMCH2Marginalized.push_back(h2);
1697  }
1698  }
1699  }
1700  return 1;
1701 }
void BCEngineMCMC::MCMCInitializeMarkovChains ( )

Initializes Markov chains.

Definition at line 1505 of file BCEngineMCMC.cxx.

1506 {
1507  // evaluate function at the starting point
1508  std::vector<double> x0;
1509 
1510  for (unsigned j = 0; j < fMCMCNChains; ++j)
1511  {
1512  x0.clear();
1513  for (unsigned i = 0; i < fParameters.Size(); ++i)
1514  x0.push_back(fMCMCx[j * fParameters.Size() + i]);
1515  fMCMCprob[j] = LogEval(x0);
1516  }
1517 
1518  x0.clear();
1519 }
void BCEngineMCMC::MCMCInitializeMarkovChainTrees ( )

Initialize trees containing the Markov chains.

Definition at line 461 of file BCEngineMCMC.cxx.

462 {
463  // clear vector
464  fMCMCTrees.clear();
465 
466  // create new trees
467  for (unsigned i = 0; i < fMCMCNChains; ++i) {
468  fMCMCTrees.push_back(new TTree(TString::Format("MarkovChainTree_%i", i), "MarkovChainTree"));
469  fMCMCTrees[i]->Branch("Iteration", &fMCMCNIterations[i], "iteration/i");
470  // todo check example and parallel_TEST how to automatically determine #parameters
471 // fMCMCTrees[i]->Branch("NParameters", fParameters.Size(), "parameters/I");
472  fMCMCTrees[i]->Branch("LogProbability", &fMCMCprob[i], "log(probability)/D");
473  fMCMCTrees[i]->Branch("Phase", &fMCMCPhase, "phase/I");
474 
475  for (unsigned j = 0; j < fParameters.Size(); ++j)
476  fMCMCTrees[i]->Branch(TString::Format("Parameter%i", j),
477  &fMCMCx[i * fParameters.Size() + j],
478  TString::Format("parameter %i/D", j));
479  }
480 }
void BCEngineMCMC::MCMCIterationInterface ( )
virtual

Interface allowing to execute arbitrary code for each iteration of the MCMC. The frequency of calling this method is influenced by the setup of the Lag and whether or not the MCMC is run with ordered parameters. This method needs to be overloaded in the derived class.

Reimplemented in BCFitter, and BCRooInterface.

Definition at line 1704 of file BCEngineMCMC.cxx.

1705 {
1706  // do user defined stuff
1707  MCMCUserIterationInterface();
1708 }
int BCEngineMCMC::MCMCMetropolis ( )

Runs Metropolis algorithm.

Definition at line 1293 of file BCEngineMCMC.cxx.

1294 {
1295  // check the number of free parameters
1296  if (GetNFreeParameters() <= 0) {
1297  BCLog::OutWarning("BCEngineMCMC::MCMCMetropolis. Number of free parameters <= 0. Do not run Metropolis.");
1298  return 0;
1299  }
1300 
1301  // check if prerun should be performed
1302  if (fMCMCFlagPreRun)
1303  MCMCMetropolisPreRun();
1304  else
1305  BCLog::OutWarning("BCEngineMCMC::MCMCMetropolis. Not running prerun. This can cause trouble if the data have changed.");
1306 
1307  // print to screen
1308  BCLog::OutSummary( "Run Metropolis MCMC...");
1309 
1310  // reset run statistics
1311  MCMCResetRunStatistics();
1312 
1313  // set phase and cycle number
1314  fMCMCPhase = 2;
1315 
1316  // perform run
1317  BCLog::OutSummary(Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));
1318 
1319  int nwrite = fMCMCNIterationsRun/10;
1320  if(nwrite < 100)
1321  nwrite=100;
1322  else if(nwrite < 500)
1323  nwrite=1000;
1324  else if(nwrite < 10000)
1325  nwrite=1000;
1326  else
1327  nwrite=10000;
1328 
1329  // start the run
1330  for (fMCMCCurrentIteration = 1; fMCMCCurrentIteration <= int(fMCMCNIterationsRun); ++fMCMCCurrentIteration)
1331  {
1332  if ( (fMCMCCurrentIteration)%nwrite == 0 )
1333  BCLog::OutDetail(Form(" --> iteration number %i (%.2f%%)", fMCMCCurrentIteration, (double)(fMCMCCurrentIteration)/(double)fMCMCNIterationsRun*100.));
1334 
1335  // if the flag is set then run over the parameters one after the other.
1336  if (fMCMCFlagOrderParameters)
1337  {
1338  // loop over parameters
1339  for (unsigned iparameters = 0; iparameters < fParameters.Size(); ++iparameters)
1340  {
1341  if (fParameters[iparameters]->Fixed())
1342  continue;
1343 
1344  // loop over chains
1345  {
1346  unsigned chunk = 1; (void) chunk;
1347  unsigned ichains; (void) ichains;
1348 #pragma omp parallel for shared(chunk) private(ichains) schedule(static, chunk)
1349  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains)
1350  {
1351  MCMCGetNewPointMetropolis(ichains, iparameters);
1352  }
1353  }
1354  // reset current chain
1355  fMCMCCurrentChain = -1;
1356 
1357  // update search for maximum
1358  MCMCInChainCheckMaximum();
1359 
1360  } // end loop over all parameters
1361 
1362  // check if the current iteration is consistent with the lag
1363  if ( fMCMCCurrentIteration % fMCMCNLag == 0)
1364  {
1365  // do anything interface
1366  MCMCIterationInterface();
1367 
1368  // fill histograms
1369  if ( ! fMCMCH1Marginalized.empty() or ! fMCMCH2Marginalized.empty())
1370  MCMCInChainFillHistograms();
1371 
1372  // write chain to file
1373  if (fMCMCFlagWriteChainToFile)
1374  MCMCInChainWriteChains();
1375  }
1376  }
1377  // if the flag is not set then run over the parameters at the same time.
1378  else
1379  {
1380  // loop over chains
1381  {
1382  unsigned chunk = 1; (void) chunk;
1383  unsigned ichains; (void) ichains;
1384 #pragma omp parallel for shared(chunk) private(ichains) schedule(static, chunk)
1385  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains)
1386  {
1387  // get new point
1388  MCMCGetNewPointMetropolis(ichains);
1389  }
1390  }
1391  // reset current chain
1392  fMCMCCurrentChain = -1;
1393 
1394  // update search for maximum
1395  MCMCInChainCheckMaximum();
1396 
1397  // check if the current iteration is consistent with the lag
1398  if (fMCMCCurrentIteration % fMCMCNLag == 0)
1399  {
1400  // do anything interface
1401  MCMCIterationInterface();
1402 
1403  // fill histograms
1404  if ( ! fMCMCH1Marginalized.empty() or ! fMCMCH2Marginalized.empty())
1405  MCMCInChainFillHistograms();
1406 
1407  // write chain to file
1408  if (fMCMCFlagWriteChainToFile)
1409  MCMCInChainWriteChains();
1410  }
1411  }
1412 
1413  } // end run
1414 
1415  // print convergence status
1416  BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));
1417 
1418  // print modes
1419 
1420  // find global maximum
1421  double probmax = fMCMCprobMax.at(0);
1422  unsigned probmaxindex = 0;
1423 
1424  // loop over all chains and find the maximum point
1425  for (unsigned i = 1; i < fMCMCNChains; ++i) {
1426  if (fMCMCprobMax.at(i) > probmax)
1427  {
1428  probmax = fMCMCprobMax.at(i);
1429  probmaxindex = i;
1430  }
1431  }
1432 
1433  // save if improved the log posterior
1434  if (fMCMCBestFitParameters.empty() || probmax > fMCMCLogMaximum) {
1435  fMCMCLogMaximum = probmax;
1436  fMCMCBestFitParameters.assign(fParameters.Size(), 0.0);
1437  for (unsigned i = 0; i < fParameters.Size(); ++i)
1438  fMCMCBestFitParameters[i] = fMCMCxMax[probmaxindex * fParameters.Size() + i];
1439  }
1440 
1441  BCLog::OutDetail(" --> Global mode from MCMC:");
1442  BCLog::OutDebug(Form(" --> Posterior value: %g", probmax));
1443  int ndigits = (int) log10(fParameters.Size());
1444  for (unsigned i = 0; i < fParameters.Size(); ++i)
1445  BCLog::OutDetail(TString::Format(" --> parameter %*i: %.4g", ndigits+1, i, fMCMCBestFitParameters[i]));
1446 
1447  // reset counter
1448  fMCMCCurrentIteration = -1;
1449 
1450  // reset current chain
1451  fMCMCCurrentChain = -1;
1452 
1453  // set flags
1454  fMCMCFlagRun = true;
1455 
1456  return 1;
1457 }
int BCEngineMCMC::MCMCMetropolisPreRun ( )

Runs a pre run for the Metropolis algorithm.

Definition at line 888 of file BCEngineMCMC.cxx.

889 {
890  // print on screen
891  BCLog::OutSummary("Pre-run Metropolis MCMC...");
892 
893  // initialize Markov chain
894  MCMCInitialize();
895  MCMCInitializeMarkovChains();
896 
897  // helper variable containing number of digits in the number of parameters
898  int ndigits = (int)log10(fParameters.Size()) +1;
899  if(ndigits<4)
900  ndigits=4;
901 
902  // reset run statistics
903  MCMCResetRunStatistics();
904  fMCMCNIterationsConvergenceGlobal = -1;
905 
906  // perform run
907  BCLog::OutSummary(Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax));
908 
909  // don't write to file during pre run
910  bool tempflag_writetofile = fMCMCFlagWriteChainToFile;
911  fMCMCFlagWriteChainToFile = false;
912 
913  // initialize counter variables and flags
914  fMCMCCurrentIteration = 1; // counts the number of iterations
915  unsigned counterupdate = 1; // after how many iterations is an update needed?
916  bool convergence = false; // convergence reached?
917  bool flagefficiency = false; // efficiency reached?
918 
919  // array of efficiencies
920  std::vector<double> efficiency;
921  efficiency.assign(fParameters.Size() * fMCMCNChains, 0.0);
922 
923  // how often to check convergence and efficiencies?
924  // it's either every fMCMCNParameters*nMCMCNIterationsUpdate (for 5 parameters the default would be 5000)
925  // or it's fMCMCNIterationsUpdateMax (10000 by default)
926  // whichever of the two is smaller
927  unsigned updateLimit = ( fMCMCNIterationsUpdateMax<fMCMCNIterationsUpdate*(fParameters.Size()) && fMCMCNIterationsUpdateMax>0 ) ?
928  fMCMCNIterationsUpdateMax : fMCMCNIterationsUpdate*(fParameters.Size());
929 
930  // loop over chains
931  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains) {
932  // loop over parameters
933  for (unsigned iparameter = 0; iparameter < fParameters.Size(); ++iparameter){
934  // global index of the parameter (throughout all the chains)
935  unsigned index = ichains * fParameters.Size() + iparameter;
936  // reset counters
937  fMCMCxMean[index] = fMCMCx[index];
938  }
939  fMCMCprobMean[ichains] = fMCMCprob[ichains];
940  }
941 
942  // set phase and cycle number
943  fMCMCPhase = 1;
944 
945  // run chain ...
946  // (a) for at least a minimum number of iterations,
947  // (b) until a maximum number of iterations is reached,
948  // (c) or until convergence is reached and the efficiency is in the
949  // specified region
950  while (fMCMCCurrentIteration < int(fMCMCNIterationsPreRunMin) ||
951  (fMCMCCurrentIteration < int(fMCMCNIterationsMax) && !(convergence && flagefficiency)))
952  {
953  //-------------------------------------------
954  // reset flags and counters
955  //-------------------------------------------
956 
957  // set convergence to false by default
958  convergence = false;
959 
960  // set number of iterations needed to converge to negative
961  fMCMCNIterationsConvergenceGlobal = -1;
962 
963  //-------------------------------------------
964  // get new point in n-dim space
965  //-------------------------------------------
966 
967  ++fMCMCNTrials;
968 
969  // if the flag is set then run over the parameters one after the other.
970  if (fMCMCFlagOrderParameters)
971  {
972  // loop over parameters
973  {
974  for (unsigned iparameters = 0; iparameters < fParameters.Size(); ++iparameters)
975  {
976  if (fParameters[iparameters]->Fixed())
977  continue;
978  // loop over chains
979 
980  unsigned chunk = 1; (void) chunk;
981  unsigned ichains; (void) ichains;
982 #pragma omp parallel for shared(chunk) private(ichains) schedule(static, chunk)
983  for (ichains = 0; ichains < fMCMCNChains; ++ichains){
984  MCMCGetNewPointMetropolis(ichains, iparameters);
985  }
986  // search for global maximum
987  MCMCInChainCheckMaximum();
988  }
989  }
990  }
991 
992  // if the flag is not set then run over the parameters at the same time.
993  else
994  {
995  // loop over chains
996  {
997  unsigned chunk = 1; (void) chunk;
998  unsigned ichains; (void) ichains;
999 #pragma omp parallel for shared(chunk) private(ichains) schedule(static, chunk)
1000  for (ichains = 0; ichains < fMCMCNChains; ++ichains){
1001  MCMCGetNewPointMetropolis(ichains);
1002  }
1003  }
1004  // search for global maximum
1005  MCMCInChainCheckMaximum();
1006  }
1007 
1008  //-------------------------------------------
1009  // print out message to log
1010  //-------------------------------------------
1011 
1012  // progress printout
1013  if ( fMCMCCurrentIteration > 0 && fMCMCCurrentIteration % fMCMCNIterationsUpdate == 0 )
1014  BCLog::OutDetail(Form(" --> Iteration %i", fMCMCNIterations[0] / GetNFreeParameters()));
1015 
1016  //-------------------------------------------
1017  // update statistics
1018  //-------------------------------------------
1019 
1020  if (counterupdate > 1)
1021  MCMCInChainUpdateStatistics();
1022 
1023  //-------------------------------------------
1024  // update scale factors and check convergence
1025  //-------------------------------------------
1026 
1027  // debugKK
1028  // check if this line makes sense
1029  if ( counterupdate % updateLimit == 0 && counterupdate > 0 && fMCMCCurrentIteration >= int(fMCMCNIterationsPreRunMin))
1030  {
1031  // -----------------------------
1032  // reset flags and counters
1033  // -----------------------------
1034 
1035  bool rvalues_ok = true;
1036 
1037  static bool has_converged = false;
1038 
1039  // reset the number of iterations needed for convergence to
1040  // negative
1041  fMCMCNIterationsConvergenceGlobal = -1;
1042 
1043  // -----------------------------
1044  // check convergence status
1045  // -----------------------------
1046 
1047  // test convergence
1048  MCMCInChainTestConvergenceAllChains();
1049 
1050  // set convergence flag
1051  if (fMCMCNIterationsConvergenceGlobal > 0)
1052  convergence = true;
1053 
1054  // print convergence status:
1055  if (convergence && fMCMCNChains > 1)
1056  BCLog::OutDetail(Form(" * Convergence status: Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
1057  else if (!convergence && fMCMCNChains > 1)
1058  {
1059  BCLog::OutDetail(Form(" * Convergence status: Set of %i Markov chains did not converge after %i iterations.", fMCMCNChains, fMCMCCurrentIteration));
1060 
1061  BCLog::OutDetail(" - R-Values:");
1062  for (unsigned iparameter = 0; iparameter < fParameters.Size(); ++iparameter)
1063  {
1064  if (fParameters[iparameter]->Fixed())
1065  continue;
1066  if(fabs(fMCMCRValueParameters[iparameter]-1.) < fMCMCRValueParametersCriterion)
1067  BCLog::OutDetail(TString::Format(" parameter %*i : %.06f",ndigits, iparameter, fMCMCRValueParameters.at(iparameter)));
1068  else
1069  {
1070  if ( fMCMCRValueParameters.at(iparameter) != std::numeric_limits<double>::max() )
1071  BCLog::OutDetail(TString::Format(" parameter %*i : %.06f <--",ndigits, iparameter, fMCMCRValueParameters.at(iparameter)));
1072  else
1073  BCLog::OutDetail(TString::Format(" parameter %*i : MAX_DOUBLE <--",ndigits, iparameter));
1074  rvalues_ok = false;
1075  }
1076  }
1077  if(fabs(fMCMCRValue-1.) < fMCMCRValueCriterion)
1078  BCLog::OutDetail(Form(" log-likelihood : %.06f", fMCMCRValue));
1079  else
1080  {
1081  if ( fMCMCRValue != std::numeric_limits<double>::max() )
1082  BCLog::OutDetail(Form(" log-likelihood : %.06f <--", fMCMCRValue));
1083  else
1084  BCLog::OutDetail(" log-likelihood : MAX_DOUBLE <--");
1085  rvalues_ok = false;
1086  }
1087  }
1088 
1089  // set convergence flag
1090  if(!has_converged)
1091  if(rvalues_ok)
1092  has_converged = true;
1093 
1094  // -----------------------------
1095  // check efficiency status
1096  // -----------------------------
1097 
1098  // set flag
1099  flagefficiency = true;
1100 
1101  bool flagprintefficiency = true;
1102 
1103  // loop over chains
1104  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains)
1105  {
1106  // loop over parameters
1107  for (unsigned iparameter = 0; iparameter < fParameters.Size(); ++iparameter)
1108  {
1109  if (fParameters[iparameter]->Fixed())
1110  continue;
1111 
1112  // global index of the parameter (throughout all the chains)
1113  unsigned index = ichains * fParameters.Size() + iparameter;
1114 
1115  // calculate efficiency
1116  efficiency[index] = double(fMCMCNTrialsTrue[index]) / double(fMCMCNTrials);
1117 
1118  // adjust scale factors if efficiency is too low
1119  if (efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01)
1120  {
1121  if (flagprintefficiency)
1122  {
1123  BCLog::OutDetail(Form(" * Efficiency status: Efficiencies not within pre-defined range."));
1124  BCLog::OutDetail(Form(" - Efficiencies:"));
1125  flagprintefficiency = false;
1126  }
1127 
1128  double fscale=2.;
1129  if(has_converged && fMCMCEfficiencyMin/efficiency[index] > 2.)
1130  fscale = 4.;
1131  fMCMCTrialFunctionScaleFactor[index] /= fscale;
1132 
1133  BCLog::OutDetail(Form(" Efficiency of parameter %i dropped below %.2f%% (eps = %.2f%%) in chain %i. Set scale to %.4g",
1134  iparameter, 100. * fMCMCEfficiencyMin, 100. * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index]));
1135  }
1136 
1137  // adjust scale factors if efficiency is too high
1138  else if (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.0)
1139  {
1140  if (flagprintefficiency)
1141  {
1142  BCLog::OutDetail(Form(" * Efficiency status: Efficiencies not within pre-defined ranges."));
1143  BCLog::OutDetail(Form(" - Efficiencies:"));
1144  flagprintefficiency = false;
1145  }
1146 
1147  fMCMCTrialFunctionScaleFactor[index] *= 2.;
1148 
1149  BCLog::OutDetail(Form(" Efficiency of parameter %i above %.2f%% (eps = %.2f%%) in chain %i. Set scale to %.4g",
1150  iparameter, 100.0 * fMCMCEfficiencyMax, 100.0 * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index]));
1151  }
1152 
1153  // check flag
1154  if ((efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01)
1155  || (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.))
1156  flagefficiency = false;
1157  } // end of running over all parameters
1158  } // end of running over all chains
1159 
1160  // print to screen
1161  if (flagefficiency)
1162  BCLog::OutDetail(Form(" * Efficiency status: Efficiencies within pre-defined ranges."));
1163 
1164  // -----------------------------
1165  // reset counters
1166  // -----------------------------
1167 
1168  counterupdate = 0;
1169 
1170  // loop over chains
1171  fMCMCNTrials = 0;
1172  for (unsigned ichains = 0; ichains < fMCMCNChains; ++ichains) {
1173  // loop over parameters
1174  for (unsigned iparameter = 0; iparameter < fParameters.Size(); ++iparameter){
1175  // global index of the parameter (throughout all the chains)
1176  unsigned index = ichains * fParameters.Size() + iparameter;
1177  // reset counters
1178  fMCMCNTrialsTrue[index] = 0;
1179  fMCMCxMean[index] = fMCMCx[index];
1180  fMCMCxVar[index] = 0;
1181  }
1182  fMCMCprobMean[ichains] = fMCMCprob[ichains];
1183  fMCMCprobVar[ichains] = 0;
1184  }
1185  } // end if update scale factors and check convergence
1186 
1187  //-------------------------------------------
1188  // write chain to file
1189  //-------------------------------------------
1190 
1191  // write chain to file
1192  if (fMCMCFlagWritePreRunToFile)
1193  MCMCInChainWriteChains();
1194 
1195  //-------------------------------------------
1196  // increase counters
1197  //-------------------------------------------
1198  fMCMCCurrentIteration++;
1199  counterupdate++;
1200 
1201  } // end of running
1202 
1203  // did we check convergence at least once ?
1204  if (fMCMCCurrentIteration < int(updateLimit))
1205  {
1206  BCLog::OutWarning(" Convergence never checked !");
1207  BCLog::OutWarning(" Increase maximum number of iterations in the pre-run /MCMCSetNIterationsMax()/");
1208  BCLog::OutWarning(" or decrease maximum number of iterations for update /MCMCSetNIterationsUpdateMax()/");
1209  }
1210 
1211  // ---------------
1212  // after chain run
1213  // ---------------
1214 
1215  // define convergence status
1216  if (fMCMCNIterationsConvergenceGlobal > 0)
1217  fMCMCFlagConvergenceGlobal = true;
1218  else
1219  fMCMCFlagConvergenceGlobal = false;
1220 
1221  // print convergence status
1222  if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && !flagefficiency)
1223  BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations but could not adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
1224 
1225  else if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && flagefficiency)
1226  BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations and all scales are adjusted.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
1227 
1228  else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && flagefficiency)
1229  BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));
1230 
1231  else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && !flagefficiency)
1232  BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations and could not adjust scales.", fMCMCNChains, fMCMCNIterationsMax));
1233 
1234  else if(fMCMCNChains == 1)
1235  BCLog::OutSummary(" --> No convergence criterion for a single chain defined.");
1236 
1237  else
1238  BCLog::OutSummary(" --> Only one Markov chain. No global convergence criterion defined.");
1239 
1240  BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCCurrentIteration));
1241 
1242 
1243  // print efficiencies
1244  std::vector<double> efficiencies;
1245 
1246  for (unsigned i = 0; i < fParameters.Size(); ++i)
1247  efficiencies.push_back(0.);
1248 
1249  BCLog::OutDetail(" --> Average efficiencies:");
1250  for (unsigned i = 0; i < fParameters.Size(); ++i)
1251  {
1252  if (fParameters[i]->Fixed())
1253  continue;
1254 
1255  for (unsigned j = 0; j < fMCMCNChains; ++j)
1256  efficiencies[i] += efficiency[j * fParameters.Size() + i] / double(fMCMCNChains);
1257 
1258  BCLog::OutDetail(TString::Format(" --> parameter %*d : %.02f%%",ndigits, i, 100. * efficiencies[i]));
1259  }
1260 
1261 
1262  // print scale factors
1263  std::vector<double> scalefactors;
1264 
1265  for (unsigned i = 0; i < fParameters.Size(); ++i)
1266  scalefactors.push_back(0.0);
1267 
1268  BCLog::OutDetail(" --> Average scale factors:");
1269  for (unsigned i = 0; i < fParameters.Size(); ++i)
1270  {
1271  if (fParameters[i]->Fixed())
1272  continue;
1273  for (unsigned j = 0; j < fMCMCNChains; ++j)
1274  scalefactors[i] += fMCMCTrialFunctionScaleFactor[j * fParameters.Size() + i] / double(fMCMCNChains);
1275 
1276  BCLog::OutDetail(TString::Format(" --> parameter %*i : %.02f%%",ndigits, i, 100. * scalefactors[i]));
1277  }
1278 
1279  // reset flag
1280  fMCMCFlagWriteChainToFile = tempflag_writetofile;
1281 
1282  // reset current iteration
1283  fMCMCCurrentIteration = -1;
1284 
1285  // reset current chain
1286  fMCMCCurrentChain = -1;
1287 
1288  // no error
1289  return 1;
1290 }
void BCEngineMCMC::MCMCResetRunStatistics ( )

Resets the run statistics.

Definition at line 1460 of file BCEngineMCMC.cxx.

1461 {
1462  fMCMCNTrials = 0;
1463 
1464  for (unsigned j = 0; j < fMCMCNChains; ++j)
1465  {
1466  fMCMCNIterations[j] = 0;
1467  fMCMCprobMean[j] = 0;
1468  fMCMCprobVar[j] = 0;
1469 
1470  for (unsigned k = 0; k < fParameters.Size(); ++k)
1471  {
1472  fMCMCNTrialsTrue[j * fParameters.Size() + k] = 0;
1473  }
1474  }
1475 
1476  // reset marginalized distributions
1477  for (unsigned i = 0; i < fMCMCH1Marginalized.size(); ++i)
1478  if (fMCMCH1Marginalized[i])
1479  fMCMCH1Marginalized[i]->Reset();
1480 
1481  for (unsigned i = 0; i < fMCMCH2Marginalized.size(); ++i)
1482  if (fMCMCH2Marginalized[i])
1483  fMCMCH2Marginalized[i]->Reset();
1484 
1485  fMCMCRValue = 100;
1486 }
void BCEngineMCMC::MCMCSetFlagFillHistograms ( bool  flag)

Sets the flag for all parameters to either fill histograms or not.

Definition at line 417 of file BCEngineMCMC.cxx.

418 {
419  for (unsigned i = 0; i < fParameters.Size(); ++i)
420  fParameters[i]->FillHistograms(flag);
421 }
void BCEngineMCMC::MCMCSetFlagInitialPosition ( int  flag)
inline

Sets flag which defines initial position.

Definition at line 399 of file BCEngineMCMC.h.

400  { fMCMCFlagInitialPosition = flag; }
void BCEngineMCMC::MCMCSetFlagOrderParameters ( bool  flag)
inline

Sets the flag which controls the sequence parameters during the running of the MCMC.

Definition at line 405 of file BCEngineMCMC.h.

406  { fMCMCFlagOrderParameters = flag; }
void BCEngineMCMC::MCMCSetFlagPreRun ( bool  flag)
inline

Sets the flag if a prerun should be performed or not.

Definition at line 412 of file BCEngineMCMC.h.

413  { fMCMCFlagPreRun = flag; }
void BCEngineMCMC::MCMCSetInitialPositions ( const std::vector< double > &  x0s)

Sets the initial positions for all chains.

Parameters
x0sinitial positions for all chains.

Definition at line 387 of file BCEngineMCMC.cxx.

388 {
389  // clear the existing initial position vector
390  fMCMCInitialPosition.clear();
391 
392  // copy the initial positions
393  unsigned n = x0s.size();
394 
395  for (unsigned i = 0; i < n; ++i)
396  fMCMCInitialPosition.push_back(x0s.at(i));
397 
398  // use these initial positions for the Markov chain
399  MCMCSetFlagInitialPosition(2);
400 }
void BCEngineMCMC::MCMCSetInitialPositions ( std::vector< std::vector< double > >  x0s)

Sets the initial positions for all chains.

Parameters
x0sinitial positions for all chains.

Definition at line 403 of file BCEngineMCMC.cxx.

404 {
405  // create new vector
406  std::vector<double> y0s;
407 
408  // loop over vector elements
409  for (unsigned i = 0; i < x0s.size(); ++i)
410  for (unsigned j = 0; j < x0s.at(i).size(); ++j)
411  y0s.push_back((x0s.at(i)).at(j));
412 
413  MCMCSetInitialPositions(y0s);
414 }
void BCEngineMCMC::MCMCSetMarkovChainTrees ( const std::vector< TTree * > &  trees)

Sets the tree containing the Markov chains.

Definition at line 424 of file BCEngineMCMC.cxx.

425 {
426  // clear vector
427  fMCMCTrees.clear();
428 
429  // copy tree
430  for (unsigned i = 0; i < trees.size(); ++i)
431  fMCMCTrees.push_back(trees[i]);
432 }
void BCEngineMCMC::MCMCSetMaximumEfficiency ( double  efficiency)
inline

Sets the maximum efficiency required for a chain.

Definition at line 370 of file BCEngineMCMC.h.

371  { fMCMCEfficiencyMax = efficiency; }
void BCEngineMCMC::MCMCSetMinimumEfficiency ( double  efficiency)
inline

Sets the minimum efficiency required for a chain.

Definition at line 365 of file BCEngineMCMC.h.

366  { fMCMCEfficiencyMin = efficiency; }
void BCEngineMCMC::MCMCSetNChains ( unsigned  n)

Sets the number of Markov chains which are run in parallel.

Definition at line 378 of file BCEngineMCMC.cxx.

379 {
380  fMCMCNChains = n;
381 
382  // re-initialize
383  MCMCInitialize();
384 }
void BCEngineMCMC::MCMCSetNIterationsMax ( unsigned  n)
inline

Sets the maximum number of iterations in the pre-run.

Definition at line 335 of file BCEngineMCMC.h.

336  { fMCMCNIterationsMax = n; }
void BCEngineMCMC::MCMCSetNIterationsPreRunMin ( unsigned  n)
inline

Sets the minimum number of iterations in the pre-run

Definition at line 345 of file BCEngineMCMC.h.

346  { fMCMCNIterationsPreRunMin = n; }
void BCEngineMCMC::MCMCSetNIterationsRun ( unsigned  n)
inline

Sets the number of iterations.

Definition at line 340 of file BCEngineMCMC.h.

341  { fMCMCNIterationsRun = n; }
void BCEngineMCMC::MCMCSetNIterationsUpdate ( unsigned  n)
inline

Sets the number of iterations in the pre-run after which an update on the statistics (convergence, efficiency, etc.) is done.

Parameters
nThe number of iterations.

Definition at line 352 of file BCEngineMCMC.h.

353  { fMCMCNIterationsUpdate = n; }
void BCEngineMCMC::MCMCSetNIterationsUpdateMax ( unsigned  n)
inline

Sets the maximum number of iterations in the pre-run after which an update on the statistics (convergence, efficiency, etc.) is done. If set to 0 no maximum is set.

Parameters
nmaximum number of iterations.

Definition at line 360 of file BCEngineMCMC.h.

361  { fMCMCNIterationsUpdateMax = n; }
void BCEngineMCMC::MCMCSetNLag ( unsigned  n)
inline

Sets the lag of the Markov chains

Definition at line 330 of file BCEngineMCMC.h.

331  { fMCMCNLag = n; }
void BCEngineMCMC::MCMCSetPrecision ( BCEngineMCMC::Precision  precision)

Set the precision for the MCMC run.

Definition at line 99 of file BCEngineMCMC.cxx.

100 {
101  switch(precision) {
102  case BCEngineMCMC::kLow:
103  fMCMCNChains = 1;
104  fMCMCNLag = 1;
105  fMCMCNIterationsMax = 10000;
106  fMCMCNIterationsRun = 10000;
107  fMCMCNIterationsPreRunMin = 100;
108  fMCMCNIterationsUpdate = 1000;
109  fMCMCNIterationsUpdateMax = 10000;
110  fMCMCRValueCriterion = 0.1;
111  fMCMCRValueParametersCriterion = 0.1;
112  fMCMCRValue = 100;
113  break;
114  case BCEngineMCMC::kMedium:
115  fMCMCNChains = 5;
116  fMCMCNLag = 1;
117  fMCMCNIterationsMax = 100000;
118  fMCMCNIterationsRun = 100000;
119  fMCMCNIterationsPreRunMin = 100;
120  fMCMCNIterationsUpdate = 1000;
121  fMCMCNIterationsUpdateMax = 10000;
122  fMCMCRValueCriterion = 0.1;
123  fMCMCRValueParametersCriterion = 0.1;
124  fMCMCRValue = 100;
125  break;
126  case BCEngineMCMC::kHigh:
127  fMCMCNChains = 10;
128  fMCMCNLag = 10;
129  fMCMCNIterationsMax = 1000000;
130  fMCMCNIterationsRun = 1000000;
131  fMCMCNIterationsPreRunMin = 100;
132  fMCMCNIterationsUpdate = 1000;
133  fMCMCNIterationsUpdateMax = 10000;
134  fMCMCRValueCriterion = 0.1;
135  fMCMCRValueParametersCriterion = 0.1;
136  fMCMCRValue = 100;
137  break;
138  case BCEngineMCMC::kVeryHigh:
139  fMCMCNChains = 10;
140  fMCMCNLag = 10;
141  fMCMCNIterationsMax = 10000000;
142  fMCMCNIterationsRun = 10000000;
143  fMCMCNIterationsPreRunMin = 100;
144  fMCMCNIterationsUpdate = 1000;
145  fMCMCNIterationsUpdateMax = 10000;
146  fMCMCRValueCriterion = 0.1;
147  fMCMCRValueParametersCriterion = 0.1;
148  fMCMCRValue = 100;
149  break;
150  }
151 
152  // re-initialize
153  MCMCInitialize();
154 }
void BCEngineMCMC::MCMCSetRandomSeed ( unsigned  seed)

Set the random number seed

Definition at line 434 of file BCEngineMCMC.cxx.

435 {
436  if (!fRandom)
437  fRandom = new TRandom3();
438 
439  // set main generator
440  fRandom->SetSeed(seed);
441 
442  // call once so return value of GetSeed() fixed
443  fRandom->Rndm();
444 
445  SyncThreadStorage();
446 
447  // type conversion to avoid compiler warnings
448  if (size_t(fMCMCNChains) != fMCMCThreadLocalStorage.size())
449  BCLog::OutError(Form("#chains does not match #(thread local storages): %d vs %u",
450  fMCMCNChains, unsigned(fMCMCThreadLocalStorage.size())));
451 
452  // set all single chain generators
453  for (unsigned i = 0; i < fMCMCNChains ; ++i){
454  // call once so return value of GetSeed() fixed
455  fMCMCThreadLocalStorage[i].rng->SetSeed(fRandom->GetSeed() + i);
456  fMCMCThreadLocalStorage[i].rng->Rndm();
457  }
458 }
void BCEngineMCMC::MCMCSetRValueCriterion ( double  r)
inline

Sets the R-value criterion for convergence of all chains.

Definition at line 417 of file BCEngineMCMC.h.

418  { fMCMCRValueCriterion = r; }
void BCEngineMCMC::MCMCSetRValueParametersCriterion ( double  r)
inline

Sets the parameter R-value criterion for convergence of all chains

Definition at line 422 of file BCEngineMCMC.h.

423  { fMCMCRValueParametersCriterion = r; }
void BCEngineMCMC::MCMCSetRValueStrict ( bool  strict = true)
inline

Use strict or relaxed rule for Gelman/Rubin R-value

Definition at line 426 of file BCEngineMCMC.h.

427  { fMCMCRValueUseStrict = strict; }
void BCEngineMCMC::MCMCSetTrialFunctionScaleFactor ( std::vector< double >  scale)
inline

Set the scale factors for the trial functions

Parameters
scalea vector of doubles containing the scale factors

Definition at line 321 of file BCEngineMCMC.h.

322  { fMCMCTrialFunctionScaleFactorStart = scale; }
void BCEngineMCMC::MCMCSetValuesDefault ( )

Set the default values for the MCMC chain.

Definition at line 39 of file BCEngineMCMC.cxx.

40 {
41  fMCMCFlagWriteChainToFile = false;
42  fMCMCFlagWritePreRunToFile = false;
43  fMCMCFlagPreRun = true;
44  fMCMCFlagRun = false;
45  fMCMCEfficiencyMin = 0.15;
46  fMCMCEfficiencyMax = 0.50;
47  fMCMCFlagInitialPosition = 1;
48  fMCMCNLag = 1;
49  fMCMCCurrentIteration = -1;
50  fMCMCCurrentChain = -1;
51  fMCMCLogMaximum = -std::numeric_limits<double>::max();
52 
53  MCMCSetValuesDetail();
54 }
void BCEngineMCMC::MCMCSetValuesDetail ( )

Set the values for a detailed MCMC run.

Definition at line 78 of file BCEngineMCMC.cxx.

79 {
80  fMCMCNChains = 5;
81  fMCMCNIterationsMax = 1000000;
82  fMCMCNIterationsRun = 100000;
83  fMCMCNIterationsPreRunMin = 500;
84  fMCMCFlagInitialPosition = 1;
85  fMCMCRValueUseStrict = false;
86  fMCMCRValueCriterion = 0.1;
87  fMCMCRValueParametersCriterion = 0.1;
88  fMCMCNIterationsConvergenceGlobal = -1;
89  fMCMCFlagConvergenceGlobal = false;
90  fMCMCRValue = 100;
91  fMCMCNIterationsUpdate = 1000;
92  fMCMCNIterationsUpdateMax = 10000;
93  fMCMCFlagOrderParameters = true;
94  fMCMCCurrentIteration = -1;
95  fMCMCCurrentChain = -1;
96 }
void BCEngineMCMC::MCMCSetValuesQuick ( )

Set the values for a quick MCMC run.

Definition at line 57 of file BCEngineMCMC.cxx.

58 {
59  fMCMCNChains = 2;
60  fMCMCNIterationsMax = 10000;
61  fMCMCNIterationsRun = 10000;
62  fMCMCNIterationsPreRunMin = 500;
63  fMCMCFlagInitialPosition = 1;
64  fMCMCRValueUseStrict = false;
65  fMCMCRValueCriterion = 0.1;
66  fMCMCRValueParametersCriterion = 0.1;
67  fMCMCNIterationsConvergenceGlobal = -1;
68  fMCMCFlagConvergenceGlobal = false;
69  fMCMCRValue = 100;
70  fMCMCNIterationsUpdate = 1000;
71  fMCMCNIterationsUpdateMax = 10000;
72  fMCMCFlagOrderParameters = true;
73  fMCMCCurrentIteration = -1;
74  fMCMCCurrentChain = -1;
75 }
void BCEngineMCMC::MCMCSetWriteChainToFile ( bool  flag)
inline

Sets flag to write Markov chains to file.

Definition at line 379 of file BCEngineMCMC.h.

380  { fMCMCFlagWriteChainToFile = flag; }
void BCEngineMCMC::MCMCSetWritePreRunToFile ( bool  flag)
inline

Sets flag to write pre run to file.

Definition at line 384 of file BCEngineMCMC.h.

385  { fMCMCFlagWritePreRunToFile = flag; }
void BCEngineMCMC::MCMCTrialFunction ( unsigned  ichain,
std::vector< double > &  x 
)
virtual

Random walk trial function. The default trial function is a Breit-Wigner. It can be overloaded by the user to set the trial function.

Parameters
ichainthe chain index
xpoint with the dimension fMCMCNParameters

Definition at line 483 of file BCEngineMCMC.cxx.

484 {
485  // call MCMCTrialFunctionSingle() for all parameters by default
486  for (unsigned i = 0; i < fParameters.Size(); ++i)
487  x[i] = MCMCTrialFunctionSingle(ichain, i);
488 }
double BCEngineMCMC::MCMCTrialFunctionSingle ( unsigned  ichain,
unsigned  ipar 
)
virtual

Random walk trial function. The default trial function is a Breit-Wigner. It can be overloaded by the user to set the trial function.

Parameters
ichainthe chain index
iparthe parameter index
Returns
the unscaled proposal point

Definition at line 491 of file BCEngineMCMC.cxx.

492 {
493  // no check of range for performance reasons
494 
495  // use uniform distribution
496  // return = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fRandom->Rndm());
497 
498  // Breit-Wigner width adjustable width
499  return fMCMCThreadLocalStorage[ichain].rng->BreitWigner(0.0,
500  fMCMCTrialFunctionScaleFactor[ichain * fParameters.Size() + iparameter]);
501 }
virtual void BCEngineMCMC::MCMCUserIterationInterface ( )
inlinevirtual

Method executed for every iteration of the MCMC. User's code should be provided via overloading in the derived class

Reimplemented in BCMTF, BCGoFTest, and BCMVCDataModel.

Definition at line 646 of file BCEngineMCMC.h.

647  {}
BCEngineMCMC & BCEngineMCMC::operator= ( const BCEngineMCMC &  engineMCMC)

Defaut assignment operator

Definition at line 257 of file BCEngineMCMC.cxx.

258 {
259  Copy(enginemcmc);
260  return *this;
261 }
void BCEngineMCMC::ResetResults ( )
virtual

Reset the MCMC variables.

Reimplemented in BCIntegrate.

Definition at line 1522 of file BCEngineMCMC.cxx.

1523 {
1524  // reset variables
1525  fMCMCNIterations.clear();
1526  fMCMCNTrialsTrue.clear();
1527  fMCMCNTrials = 0;
1528  fMCMCTrialFunctionScaleFactor.clear();
1529  fMCMCprobMean.clear();
1530  fMCMCprobVar.clear();
1531  fMCMCxMean.clear();
1532  fMCMCxVar.clear();
1533  fMCMCx.clear();
1534  fMCMCprob.clear();
1535  fMCMCxMax.clear();
1536  fMCMCprobMax.clear();
1537  fMCMCNIterationsConvergenceGlobal = -1;
1538  fMCMCRValueParameters.clear();
1539 
1540  for (unsigned i = 0; i < fMCMCH1Marginalized.size(); ++i)
1541  if (fMCMCH1Marginalized[i])
1542  delete fMCMCH1Marginalized[i];
1543 
1544  for (unsigned i = 0; i < fMCMCH2Marginalized.size(); ++i)
1545  if (fMCMCH2Marginalized[i])
1546  delete fMCMCH2Marginalized[i];
1547 
1548  // clear plots
1549  fMCMCH1Marginalized.clear();
1550  fMCMCH2Marginalized.clear();
1551 
1552  // reset flags
1553  fMCMCFlagPreRun = true;
1554  fMCMCFlagRun = false;
1555  fMCMCFlagConvergenceGlobal = false;
1556 
1557  fMCMCBestFitParameters.clear();
1558  fMCMCLogMaximum = -std::numeric_limits<double>::max();
1559  fMarginalModes.clear();
1560 }
int BCEngineMCMC::SetMarginalized ( unsigned  index,
TH1D *  h 
)

Sets the histogram with 1D marginalized distributions for parameter.

Parameters
iindex of the parameter
hpointer to an existing histogram

Definition at line 1711 of file BCEngineMCMC.cxx.

1712 {
1713  if(fMCMCH1Marginalized.size() <= index)
1714  return 0;
1715 
1716  if(h==0)
1717  return 0;
1718 
1719  if(fMCMCH1Marginalized.size() == index)
1720  fMCMCH1Marginalized.push_back(h);
1721  else
1722  fMCMCH1Marginalized[index]=h;
1723 
1724  return index;
1725 }
int BCEngineMCMC::SetMarginalized ( unsigned  index1,
unsigned  index2,
TH2D *  h 
)

Sets the histogram with 2D marginalized distributions for two parameters.

Parameters
index1index of the first parameter
index2index of the second parameter
hpointer to an existing histogram

Definition at line 1728 of file BCEngineMCMC.cxx.

1729 {
1730  unsigned counter = 0;
1731  unsigned index = 0;
1732 
1733  // search for correct combination
1734  for(unsigned i = 0; i < fParameters.Size(); i++)
1735  for (unsigned j = 0; j < i; ++j)
1736  {
1737  if(j == index1 && i == index2)
1738  index = counter;
1739  counter++;
1740  }
1741 
1742  if(fMCMCH2Marginalized.size()<=index)
1743  return 0;
1744 
1745  if(h==0)
1746  return 0;
1747 
1748  if(fMCMCH2Marginalized.size()==index)
1749  fMCMCH2Marginalized.push_back(h);
1750  else
1751  fMCMCH2Marginalized[index]=h;
1752 
1753  return index;
1754 }
void BCEngineMCMC::SetNbins ( unsigned int  nbins)

Set the number of bins for the marginalized distribution of all parameters

Parameters
nbinsNumber of bins

Definition at line 282 of file BCEngineMCMC.cxx.

283 {
284  for (unsigned i = 0 ; i < fParameters.Size() ; ++i)
285  fParameters[i]->SetNbins(nbins);
286 }
void BCEngineMCMC::SyncThreadStorage ( )
private

Ensure that there are as many storages as chains

Definition at line 1791 of file BCEngineMCMC.cxx.

1792 {
1793  // always need as many local storage as #chains
1794  const int n = fMCMCThreadLocalStorage.size() - fMCMCNChains;
1795  if (n < 0)
1796  {
1797  // fix return value of GetSeed()
1798  fRandom->Rndm();
1799 
1800  for (int i = 0; i < -n; ++i){
1801  // append one new storage
1802  fMCMCThreadLocalStorage.push_back(MCMCThreadLocalStorage(fParameters.Size()));
1803  // each chains gets a different seed
1804  // We assume that fRandom always returns same seed, presumably as it has generated at least one random number
1805  fMCMCThreadLocalStorage.back().rng->SetSeed(fRandom->GetSeed() + fMCMCThreadLocalStorage.size());
1806  }
1807  }
1808  else if (n > 0)
1809  {
1810  for (int i = 0; i < n; ++i){
1811  fMCMCThreadLocalStorage.pop_back();
1812  }
1813  }
1814 
1815  // update parameter size for each chain
1816  for (unsigned i = 0 ; i < fMCMCThreadLocalStorage.size(); ++i)
1817  fMCMCThreadLocalStorage[i].xLocal.assign(fParameters.Size(), 0.0);
1818 }
void BCEngineMCMC::WriteMarkovChain ( bool  flag)
inline

Flag for writing Markov chain to ROOT file (true) or not (false)

Definition at line 487 of file BCEngineMCMC.h.

488  {
489  fMCMCFlagWriteChainToFile = flag;
490  fMCMCFlagWritePreRunToFile = flag;
491  }

Member Data Documentation

bool BCEngineMCMC::fFlagMarginalized
private

Definition at line 711 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMarginalModes
protected

A vector of best fit parameters estimated from the marginalized probability

Definition at line 926 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCBestFitParameters
protected

A vector of best fit parameters found by MCMC

Definition at line 918 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCCurrentChain
protected

The current chain index. If not called within the running of the algorithm, return -1.

Definition at line 739 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCCurrentIteration
protected

The current iteration number. If not called within the running of the algorithm, return -1.

Definition at line 734 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCEfficiencyMax
protected

The maximum allowed efficiency for MCMC

Definition at line 817 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCEfficiencyMin
protected

The minimum required efficiency for MCMC

Definition at line 813 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagConvergenceGlobal
protected

Flag for convergence

Definition at line 756 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCFlagInitialPosition
protected

Variable which defines the initial position. 0 (default) center of the allowed region, (1) random initial position (2) pre-defined intial position.

Definition at line 823 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagOrderParameters
protected

Flag which controls the sequence parameters during the running of the MCMC.

Definition at line 828 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagPreRun
protected

Defines if a prerun should be performed or not

Definition at line 798 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagRun
protected

Defines if MCMC has been performed or not

Definition at line 802 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagWriteChainToFile
protected

Flag to write Markov chains to file

Definition at line 781 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagWritePreRunToFile
protected

Flag to write pre run to file

Definition at line 785 of file BCEngineMCMC.h.

std::vector<TH1D *> BCEngineMCMC::fMCMCH1Marginalized
protected

An array of marginalized distributions

Definition at line 908 of file BCEngineMCMC.h.

std::vector<TH2D *> BCEngineMCMC::fMCMCH2Marginalized
protected

Definition at line 909 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCInitialPosition
protected

The intial position of each Markov chain. The length of the vectors is equal to fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on

Definition at line 809 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCLogMaximum
protected

The function value at the mode on the log scale

Definition at line 922 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNChains
protected

Number of Markov chains ran in parallel

Definition at line 720 of file BCEngineMCMC.h.

std::vector<unsigned> BCEngineMCMC::fMCMCNIterations
protected

Number of total iterations of the Markov chains. The length of the vector is equal to fMCMCNChains.

Definition at line 729 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsConvergenceGlobal
protected

Number of iterations needed for all chains to convergence simultaneously

Definition at line 752 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNIterationsMax
protected

Maximum number of iterations for a Markov chain prerun

Definition at line 760 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNIterationsPreRunMin
protected

Minimum number of iterations for the pre-run

Definition at line 768 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNIterationsRun
protected

Number of iterations for a Markov chain run

Definition at line 764 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNIterationsUpdate
protected

Number of iterations for updating scale factors

Definition at line 743 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNIterationsUpdateMax
protected

Maximum number of iterations for updating scale factors

Definition at line 747 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNLag
protected

The lag for the Markov Chain

Definition at line 724 of file BCEngineMCMC.h.

unsigned BCEngineMCMC::fMCMCNTrials
protected

Number of trials

Definition at line 777 of file BCEngineMCMC.h.

std::vector<int> BCEngineMCMC::fMCMCNTrialsTrue
protected

Number of accepted trials for each chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.

Definition at line 773 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCPhase
protected

The phase of the run. 1: pre-run, 2: main run.

Definition at line 834 of file BCEngineMCMC.h.

MCMCPointerToGetProposalPoint BCEngineMCMC::fMCMCPointerToGetProposalPoint
private

Pointer to a member function

Definition at line 671 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCprob
protected

The log of the probability of the current points of each Markov chain. The length of the vectors is fMCMCNChains.

Definition at line 862 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCprobMax
protected

The maximum (log) probability of each Markov chain. The length of the vector is fMCMCNChains.

Definition at line 867 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCprobMean
protected

The mean of all log prob values of each Markov chain. The length of the vector is equal to fMCMCNChains.

Definition at line 872 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCprobVar
protected

The variance of all log prob values of each Markov chain. The length of the vector is equal to fMCMCNChains.

Definition at line 877 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValue
protected

The R-value at which the chains did converge

Definition at line 897 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValueCriterion
protected

The R-value criterion for convergence of log-likelihood

Definition at line 889 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCRValueParameters
protected

The R-values for each parameter

Definition at line 900 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValueParametersCriterion
protected

The R-value criterion for convergence of parameters

Definition at line 893 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCRValueUseStrict
protected

flag: use exactly the R-value definition of Gelman/Rubin (R_strict) or a relaxed, simplified version (R_relax) [default]. Note that R_relax <= R_strict, and in some cases even R_relax < 1, but we always have R_strict >= 1.

Definition at line 885 of file BCEngineMCMC.h.

std::vector<MCMCThreadLocalStorage> BCEngineMCMC::fMCMCThreadLocalStorage
private

Keep thread local variables private.

Definition at line 704 of file BCEngineMCMC.h.

std::vector<TTree *> BCEngineMCMC::fMCMCTrees
protected

The trees containing the Markov chains. The length of the vector is fMCMCNChains.

Definition at line 914 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCTrialFunctionScaleFactor
protected

Scales the width of the trial functions by a scale factor for each parameter and chain

Definition at line 790 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCTrialFunctionScaleFactorStart
protected

Start values of the scale factors for the trial functions.

Definition at line 794 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCx
protected

The current points of each Markov chain. The length of the vectors is equal to fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on.

Definition at line 841 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCxMax
protected

The maximum points of each Markov chain. The length of the vector is fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on.

Definition at line 847 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCxMean
protected

The mean of all parameters of each Markov chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.

Definition at line 852 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCxVar
protected

The variance of all parameters of each Markov chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.

Definition at line 857 of file BCEngineMCMC.h.

BCParameterSet BCEngineMCMC::fParameters
protected

Parameter settings

Definition at line 716 of file BCEngineMCMC.h.

TRandom3* BCEngineMCMC::fRandom
protected

Random number generator

Definition at line 904 of file BCEngineMCMC.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.4