Protected Attributes | Private Types | Private Attributes

BCEngineMCMC Class Reference

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

#include <BCEngineMCMC.h>

Inherited by BCIntegrate.

List of all members.

Public Types

Enumerators

enum  Precision { kLow, kMedium, kHigh, kVeryHigh }

Public Member Functions

Constructors and destructors

 BCEngineMCMC ()
 BCEngineMCMC (int n)
 BCEngineMCMC (const BCEngineMCMC &enginemcmc)
virtual ~BCEngineMCMC ()
Assignment operators

BCEngineMCMCoperator= (const BCEngineMCMC &engineMCMC)
Getters

int MCMCGetNParameters ()
int MCMCGetNChains ()
int MCMCGetNLag ()
std::vector< int > MCMCGetNIterations ()
int MCMCGetCurrentIteration ()
int MCMCGetCurrentChain ()
int MCMCGetNIterationsConvergenceGlobal ()
bool MCMCGetFlagConvergenceGlobal ()
int MCMCGetNIterationsMax ()
int MCMCGetNIterationsRun ()
int MCMCGetNIterationsPreRunMin ()
int MCMCGetNIterationsUpdate ()
int MCMCGetNIterationsUpdateMax ()
std::vector< int > MCMCGetNTrialsTrue ()
std::vector< int > MCMCGetNTrialsFalse ()
std::vector< double > MCMCGetprobMean ()
std::vector< double > MCMCGetVariance ()
std::vector< double > MCMCGetTrialFunctionScaleFactor ()
std::vector< double > MCMCGetTrialFunctionScaleFactor (int ichain)
double MCMCGetTrialFunctionScaleFactor (int ichain, int ipar)
std::vector< double > MCMCGetx ()
std::vector< double > MCMCGetx (int ichain)
double MCMCGetx (int ichain, int ipar)
std::vector< double > MCMCGetLogProbx ()
double MCMCGetLogProbx (int ichain)
int MCMCGetPhase ()
int MCMCGetCycle ()
std::vector< double > MCMCGetMaximumPoints ()
std::vector< double > MCMCGetMaximumPoint (int i)
std::vector< double > MCMCGetMaximumLogProb ()
int MCMCGetFlagInitialPosition ()
double MCMCGetRValueCriterion ()
double MCMCGetRValueParametersCriterion ()
double MCMCGetRValue ()
double MCMCGetRValueParameters (int i)
bool MCMCGetRValueStrict ()
bool MCMCGetFlagRun ()
TTree * MCMCGetMarkovChainTree (int i)
TH1D * MCMCGetH1Marginalized (int i)
TH2D * MCMCGetH2Marginalized (int i, int j)
TRandom3 * MCMCGetTRandom3 ()
Setters

void MCMCSetTrialFunctionScaleFactor (std::vector< double > scale)
void MCMCSetNChains (int n)
void MCMCSetNLag (int n)
void MCMCSetNIterationsMax (int n)
void MCMCSetNIterationsRun (int n)
void MCMCSetNIterationsPreRunMin (int n)
void MCMCSetNIterationsUpdate (int n)
void MCMCSetNIterationsUpdateMax (int n)
void MCMCSetMinimumEfficiency (double efficiency)
void MCMCSetMaximumEfficiency (double efficiency)
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 MCMCSetFlagFillHistograms (int index, bool flag)
void MCMCSetFlagPreRun (bool flag)
void MCMCSetRValueCriterion (double r)
void MCMCSetRValueParametersCriterion (double r)
void MCMCSetRValueStrict (bool strict=true)
void MCMCSetMarkovChainTrees (std::vector< TTree * > trees)
void MCMCInitializeMarkovChainTrees ()
int SetMarginalized (int index, TH1D *h)
int SetMarginalized (int index1, int index2, TH2D *h)
void MCMCSetValuesDefault ()
void MCMCSetValuesQuick ()
void MCMCSetValuesDetail ()
void MCMCSetPrecision (BCEngineMCMC::Precision precision)
Miscellaneous methods

int MCMCAddParameter (double min, double max)
virtual void MCMCTrialFunction (int ichain, std::vector< double > &x)
virtual double MCMCTrialFunctionSingle (int ichain, int ipar)
bool MCMCGetProposalPointMetropolis (int chain, std::vector< double > &x)
bool MCMCGetProposalPointMetropolis (int chain, int parameter, std::vector< double > &x)
bool MCMCGetNewPointMetropolis (int chain=0)
bool MCMCGetNewPointMetropolis (int chain, int 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 ()
int MCMCResetResults ()
virtual void MCMCIterationInterface ()
virtual void MCMCCurrentPointInterface (std::vector< double > &point, int ichain, bool accepted)

Protected Attributes

int fMCMCNParameters
std::vector< double > fMCMCBoundaryMin
std::vector< double > fMCMCBoundaryMax
std::vector< bool > fMCMCFlagsFillHistograms
int fMCMCNChains
int fMCMCNLag
std::vector< int > fMCMCNIterations
int fMCMCCurrentIteration
int fMCMCCurrentChain
int fMCMCNIterationsUpdate
int fMCMCNIterationsUpdateMax
int fMCMCNIterationsConvergenceGlobal
bool fMCMCFlagConvergenceGlobal
int fMCMCNIterationsMax
int fMCMCNIterationsRun
int fMCMCNIterationsPreRunMin
std::vector< int > fMCMCNTrialsTrue
std::vector< int > fMCMCNTrialsFalse
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
bool fMCMCFlagFillHistograms
int fMCMCPhase
int fMCMCCycle
std::vector< double > fMCMCx
std::vector< double > fMCMCxMax
std::vector< double > fMCMCxMean
std::vector< double > fMCMCxVar
std::vector< double > fMCMCxLocal
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< int > fMCMCH1NBins
std::vector< TH1D * > fMCMCH1Marginalized
std::vector< TH2D * > fMCMCH2Marginalized
std::vector< TTree * > fMCMCTrees

Private Types

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

Private Attributes

MCMCPointerToGetProposalPoint fMCMCPointerToGetProposalPoint

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 35 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 579 of file BCEngineMCMC.h.


Member Enumeration Documentation

An enumerator for the status of a test.

Enumerator:
kLow 
kMedium 
kHigh 
kVeryHigh 

Definition at line 44 of file BCEngineMCMC.h.

{ kLow, kMedium, kHigh, kVeryHigh };


Constructor & Destructor Documentation

BCEngineMCMC::BCEngineMCMC (  ) 

Default constructor.

Definition at line 24 of file BCEngineMCMC.cxx.

{
   // set default parameters for the mcmc
   MCMCSetValuesDefault();

   // initialize random number generator
   fRandom = new TRandom3(0);
}

BCEngineMCMC::BCEngineMCMC ( int  n  ) 

Constructor.

Parameters:
n number of chains

Definition at line 34 of file BCEngineMCMC.cxx.

{
   // set number of chains to n
   fMCMCNChains = n;

   // call default constructor
   BCEngineMCMC();
}

BCEngineMCMC::BCEngineMCMC ( const BCEngineMCMC enginemcmc  ) 

Default copy constructor.

Definition at line 183 of file BCEngineMCMC.cxx.

{
   fMCMCPointerToGetProposalPoint = enginemcmc.fMCMCPointerToGetProposalPoint;
   fMCMCNParameters           = enginemcmc.fMCMCNParameters;
   fMCMCBoundaryMin           = enginemcmc.fMCMCBoundaryMin;
   fMCMCBoundaryMax           = enginemcmc.fMCMCBoundaryMax;
   fMCMCFlagsFillHistograms   = enginemcmc.fMCMCFlagsFillHistograms;
   fMCMCNChains               = enginemcmc.fMCMCNChains;
   fMCMCNLag                  = enginemcmc.fMCMCNLag;
   fMCMCNIterations           = enginemcmc.fMCMCNIterations;
   fMCMCCurrentIteration      = enginemcmc.fMCMCCurrentIteration;
   fMCMCCurrentChain          = enginemcmc.fMCMCCurrentChain;
   fMCMCNIterationsUpdate     = enginemcmc.fMCMCNIterationsUpdate;
   fMCMCNIterationsUpdateMax  = enginemcmc.fMCMCNIterationsUpdateMax;
   fMCMCNIterationsConvergenceGlobal = enginemcmc.fMCMCNIterationsConvergenceGlobal;
   fMCMCFlagConvergenceGlobal = enginemcmc.fMCMCFlagConvergenceGlobal;
   fMCMCNIterationsMax        = enginemcmc.fMCMCNIterationsMax;
   fMCMCNIterationsRun        = enginemcmc.fMCMCNIterationsRun;
   fMCMCNIterationsPreRunMin  = enginemcmc.fMCMCNIterationsPreRunMin;
   fMCMCNTrialsTrue           = enginemcmc.fMCMCNTrialsTrue;
   fMCMCNTrialsFalse          = enginemcmc.fMCMCNTrialsFalse;
   fMCMCFlagWriteChainToFile  = enginemcmc.fMCMCFlagWriteChainToFile;
   fMCMCFlagWritePreRunToFile = enginemcmc.fMCMCFlagWritePreRunToFile;
   fMCMCTrialFunctionScaleFactor = enginemcmc.fMCMCTrialFunctionScaleFactor;
   fMCMCTrialFunctionScaleFactorStart = enginemcmc.fMCMCTrialFunctionScaleFactorStart;
   fMCMCFlagPreRun            = enginemcmc.fMCMCFlagPreRun;
   fMCMCFlagRun               = enginemcmc.fMCMCFlagRun;
   fMCMCInitialPosition       = enginemcmc.fMCMCInitialPosition;
   fMCMCEfficiencyMin         = enginemcmc.fMCMCEfficiencyMin;
   fMCMCEfficiencyMax         = enginemcmc.fMCMCEfficiencyMax;
   fMCMCFlagInitialPosition   = enginemcmc.fMCMCFlagInitialPosition;
   fMCMCFlagOrderParameters   = enginemcmc.fMCMCFlagOrderParameters;
   fMCMCFlagFillHistograms    = enginemcmc.fMCMCFlagFillHistograms;
   fMCMCPhase                 = enginemcmc.fMCMCPhase;
   fMCMCCycle                 = enginemcmc.fMCMCCycle;
   fMCMCx                     = enginemcmc.fMCMCx;
   fMCMCxMax                  = enginemcmc.fMCMCxMax;
   fMCMCxMean                 = enginemcmc.fMCMCxMean;
   fMCMCxVar                  = enginemcmc.fMCMCxVar;
   fMCMCxLocal                = enginemcmc.fMCMCxLocal;
   fMCMCprob                  = enginemcmc.fMCMCprob;
   fMCMCprobMax               = enginemcmc.fMCMCprobMax;
   fMCMCprobMean              = enginemcmc.fMCMCprobMean;
   fMCMCprobVar               = enginemcmc.fMCMCprobVar;
   fMCMCRValueUseStrict       = enginemcmc.fMCMCRValueUseStrict;
   fMCMCRValueCriterion       = enginemcmc.fMCMCRValueCriterion ;
   fMCMCRValueParametersCriterion = enginemcmc.fMCMCRValueParametersCriterion;
   fMCMCRValue                = enginemcmc.fMCMCRValue;
   fMCMCRValueParameters      = enginemcmc.fMCMCRValueParameters;
   if (enginemcmc.fRandom)
      fRandom = new TRandom3(*enginemcmc.fRandom);
   else
      fRandom = 0;
   fMCMCH1NBins               = enginemcmc.fMCMCH1NBins;

   for (int i = 0; i < int(enginemcmc.fMCMCH1Marginalized.size()); ++i) {
      if (enginemcmc.fMCMCH1Marginalized.at(i))
         fMCMCH1Marginalized.push_back(new TH1D(*(enginemcmc.fMCMCH1Marginalized.at(i))));
      else
         fMCMCH1Marginalized.push_back(0);
   }

   for (int i = 0; i < int(enginemcmc.fMCMCH2Marginalized.size()); ++i) {
      if (enginemcmc.fMCMCH2Marginalized.at(i))
         fMCMCH2Marginalized.push_back(new TH2D(*(enginemcmc.fMCMCH2Marginalized.at(i))));
      else
         fMCMCH2Marginalized.push_back(0);
   }

   for (int i = 0; i < int(enginemcmc.fMCMCTrees.size()); ++i) {
      // debugKK
      fMCMCTrees.push_back(0);
   }
}

BCEngineMCMC::~BCEngineMCMC (  )  [virtual]

Default destructor.

Definition at line 163 of file BCEngineMCMC.cxx.

{
   // delete random number generator
   if (fRandom)
      delete fRandom;

   // delete 1-d marginalized distributions
   for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
      if (fMCMCH1Marginalized[i])
         delete fMCMCH1Marginalized[i];
   fMCMCH1Marginalized.clear();

   // delete 2-d marginalized distributions
   for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
      if (fMCMCH2Marginalized[i])
         delete fMCMCH2Marginalized[i];
   fMCMCH2Marginalized.clear();
}


Member Function Documentation

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 900 of file BCEngineMCMC.cxx.

{
   // test function for now
   // this will be overloaded by the user
   return 0.0;
}

int BCEngineMCMC::MCMCAddParameter ( double  min,
double  max 
)

Adds a parameter.

Parameters:
min minimum value of the parameter
max maximum value of the parameter
Returns:
number of parameters after adding

Definition at line 1474 of file BCEngineMCMC.cxx.

{
   // add the boundaries to the corresponding vectors
   fMCMCBoundaryMin.push_back(min);
   fMCMCBoundaryMax.push_back(max);

   // set flag for individual parameters
   fMCMCFlagsFillHistograms.push_back(true);

   // increase the number of parameters by one
   fMCMCNParameters++;

   // return the number of parameters
   return fMCMCNParameters;
}

virtual void BCEngineMCMC::MCMCCurrentPointInterface ( std::vector< double > &  point,
int  ichain,
bool  accepted 
) [inline, virtual]

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

Parameters:
point point that was generated and checked
ichain index of the chain
accepted flag whether or not the point was accepted for the chain

Definition at line 564 of file BCEngineMCMC.h.

         {
            // suppress warnings for unused parameters
            // with optimization, no code should be generated
            (void)point;
            (void)ichain;
            (void)accepted;
         }

int BCEngineMCMC::MCMCGetCurrentChain (  )  [inline]
Returns:
current chain index

Definition at line 106 of file BCEngineMCMC.h.

         { return fMCMCCurrentChain; }

int BCEngineMCMC::MCMCGetCurrentIteration (  )  [inline]
Returns:
current iterations

Definition at line 101 of file BCEngineMCMC.h.

         { return fMCMCCurrentIteration; }

int BCEngineMCMC::MCMCGetCycle (  )  [inline]
Returns:
pointer to the cycle of a pre-run.

Definition at line 216 of file BCEngineMCMC.h.

         { return fMCMCCycle; }

bool BCEngineMCMC::MCMCGetFlagConvergenceGlobal (  )  [inline]
Returns:
flag if converged or not

Definition at line 117 of file BCEngineMCMC.h.

int BCEngineMCMC::MCMCGetFlagInitialPosition (  )  [inline]
Returns:
flag which defined initial position

Definition at line 236 of file BCEngineMCMC.h.

bool BCEngineMCMC::MCMCGetFlagRun (  )  [inline]
Returns:
the flag if MCMC has been performed or not

Definition at line 266 of file BCEngineMCMC.h.

      { return fMCMCFlagRun; }

TH1D * BCEngineMCMC::MCMCGetH1Marginalized ( int  i  ) 

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

Parameters:
i index of the parameter
Returns:
pointer to the histogram

Definition at line 339 of file BCEngineMCMC.cxx.

{
      // check index
   if (index < 0 || index >= int(fMCMCH1Marginalized.size()))
   {
      BCLog::OutWarning("BCEngineMCMC::MCMCGetH1Marginalized. Index out of range.");
      return 0;
   }

   return fMCMCH1Marginalized[index];
}

TH2D * BCEngineMCMC::MCMCGetH2Marginalized ( int  i,
int  j 
)

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

Parameters:
i index of the first parameter
j index of the second parameter
Returns:
pointer to the histogram

Definition at line 352 of file BCEngineMCMC.cxx.

{
   int counter = 0;
   int index = 0;

   // search for correct combination
   for(int i = 0; i < fMCMCNParameters; i++)
      for (int j = 0; j < i; ++j)
      {
         if(j == index1 && i == index2)
            index = counter;
         counter++;
      }

   // check index
   if (index < 0 || index >= int(fMCMCH2Marginalized.size()))
   {
      BCLog::OutWarning("BCEngineMCMC::MCMCGetH2Marginalized. Index out of range.");
      return 0;
   }

   return fMCMCH2Marginalized[index];
}

std::vector<double> BCEngineMCMC::MCMCGetLogProbx (  )  [inline]
Returns:
log of the probability of the current points of each Markov chain

Definition at line 201 of file BCEngineMCMC.h.

         { return fMCMCprob; }

double BCEngineMCMC::MCMCGetLogProbx ( int  ichain  ) 
Returns:
log of the probability of the current points of the Markov chain.
Parameters:
ichain chain index

Definition at line 573 of file BCEngineMCMC.cxx.

{
   // check if ichain is in range
   if (ichain < 0 || ichain >= fMCMCNChains)
      return -1;

   // return log of the probability at the current point in the ith chain
   return fMCMCprob.at(ichain);
}

TTree* BCEngineMCMC::MCMCGetMarkovChainTree ( int  i  )  [inline]

Rtrieve the tree containing the Markov chain.

Parameters:
i index of the Markov chain
Returns:
pointer to the tree

Definition at line 273 of file BCEngineMCMC.h.

         { return fMCMCTrees.at(i); }

std::vector<double> BCEngineMCMC::MCMCGetMaximumLogProb (  )  [inline]
Returns:
maximum (log) probability of each Markov chain

Definition at line 231 of file BCEngineMCMC.h.

         { return fMCMCprobMax; }

std::vector< double > BCEngineMCMC::MCMCGetMaximumPoint ( int  i  ) 
Returns:
maximum point of Markov chain
Parameters:
i The index of the Markov chain

Definition at line 377 of file BCEngineMCMC.cxx.

{
   // create a new vector with the lenght of fMCMCNParameters
   std::vector<double> x;

   // check if i is in range
   if (i < 0 || i >= fMCMCNChains)
      return x;

   // copy the point in the ith chain into the temporary vector
   for (int j = 0; j < fMCMCNParameters; ++j)
   x.push_back(fMCMCxMax.at(i * fMCMCNParameters + j));

   return x;
}

std::vector<double> BCEngineMCMC::MCMCGetMaximumPoints (  )  [inline]
Returns:
maximum points of each Markov chain

Definition at line 221 of file BCEngineMCMC.h.

         { return fMCMCxMax; }

int BCEngineMCMC::MCMCGetNChains (  )  [inline]
Returns:
number of Markov chains

Definition at line 86 of file BCEngineMCMC.h.

         { return fMCMCNChains; }

bool BCEngineMCMC::MCMCGetNewPointMetropolis ( int  chain = 0  ) 

Generates a new point using the Metropolis algorithm.

Parameters:
chain chain index

Definition at line 693 of file BCEngineMCMC.cxx.

{
   // calculate index
   int index = chain * fMCMCNParameters;

   fMCMCCurrentChain = chain;

   // increase counter
   fMCMCNIterations[chain]++;

   // get proposal point
   if (!MCMCGetProposalPointMetropolis(chain, fMCMCxLocal))
   {
      // increase counter
      for (int i = 0; i < fMCMCNParameters; ++i)
         fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;

      // execute user code for every point
      MCMCCurrentPointInterface(fMCMCxLocal, chain, false);

      return false;
   }

   // calculate probabilities of the old and new points
   double p0 = fMCMCprob[chain];
   double p1 = LogEval(fMCMCxLocal);

   // flag for accept
   bool accept = false;

   // if the new point is more probable, keep it ...
   if (p1 >= p0)
      accept = true;

   // ... or else throw dice.
   else
   {
      double r = log(fRandom->Rndm());

      if(r < p1 - p0)
         accept = true;
   }

   // fill the new point
   if(accept)
   {
      // increase counter
      for (int i = 0; i < fMCMCNParameters; ++i)
         fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;

      // copy the point
      for(int i = 0; i < fMCMCNParameters; ++i)
      {
         // save the point
         fMCMCx[index + i] = fMCMCxLocal[i];

         // save the probability of the point
         fMCMCprob[chain] = p1;
      }
   }
   else
   {
      // increase counter
      for (int i = 0; i < fMCMCNParameters; ++i)
         fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
   }

   // execute user code for every point
   MCMCCurrentPointInterface(fMCMCxLocal, chain, accept);

   return accept;
}

bool BCEngineMCMC::MCMCGetNewPointMetropolis ( int  chain,
int  parameter 
)

Definition at line 623 of file BCEngineMCMC.cxx.

{
   // calculate index
   int index = chain * fMCMCNParameters;

   fMCMCCurrentChain = chain;

   // increase counter
   fMCMCNIterations[chain]++;

   // get proposal point
   if (!MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal))
   {
      // increase counter
      fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;

      // execute user code for every point
      MCMCCurrentPointInterface(fMCMCxLocal, chain, false);

      return false;
   }

   // calculate probabilities of the old and new points
   double p0 = fMCMCprob[chain];
   double p1 = LogEval(fMCMCxLocal);

   // flag for accept
   bool accept = false;

   // if the new point is more probable, keep it ...
   if (p1 >= p0)
      accept = true;
   // ... or else throw dice.
   else
   {
      double r = log(fRandom->Rndm());

      if(r < p1 - p0)
         accept = true;
   }

   // fill the new point
   if(accept)
   {
      // increase counter
      fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++;

      // copy the point
      for(int i = 0; i < fMCMCNParameters; ++i)
      {
         // save the point
         fMCMCx[index + i] = fMCMCxLocal[i];

         // save the probability of the point
         fMCMCprob[chain] = p1;
      }
   }
   else
   {
      // increase counter
      fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
   }

   // execute user code for every point
   MCMCCurrentPointInterface(fMCMCxLocal, chain, accept);

   return accept;
}

std::vector<int> BCEngineMCMC::MCMCGetNIterations (  )  [inline]
Returns:
number of iterations

Definition at line 96 of file BCEngineMCMC.h.

         { return fMCMCNIterations; }

int BCEngineMCMC::MCMCGetNIterationsConvergenceGlobal (  )  [inline]
Returns:
number of iterations needed for all chains to converge simultaneously

Definition at line 112 of file BCEngineMCMC.h.

int BCEngineMCMC::MCMCGetNIterationsMax (  )  [inline]
Returns:
maximum number of iterations for a Markov chain

Definition at line 122 of file BCEngineMCMC.h.

         { return fMCMCNIterationsMax; }

int BCEngineMCMC::MCMCGetNIterationsPreRunMin (  )  [inline]
Returns:
minimum number of pre-run iterations for a Markov chain

Definition at line 132 of file BCEngineMCMC.h.

int BCEngineMCMC::MCMCGetNIterationsRun (  )  [inline]
Returns:
number of iterations for a Markov chain

Definition at line 127 of file BCEngineMCMC.h.

         { return fMCMCNIterationsRun; }

int BCEngineMCMC::MCMCGetNIterationsUpdate (  )  [inline]
Returns:
number of iterations after statistics update.

Definition at line 137 of file BCEngineMCMC.h.

int BCEngineMCMC::MCMCGetNIterationsUpdateMax (  )  [inline]
Returns:
maximum number of iterations after statistics update.

Definition at line 142 of file BCEngineMCMC.h.

int BCEngineMCMC::MCMCGetNLag (  )  [inline]
Returns:
lag of the Markov chains

Definition at line 91 of file BCEngineMCMC.h.

         { return fMCMCNLag; }

int BCEngineMCMC::MCMCGetNParameters (  )  [inline]
Returns:
number of parameters of the Markov chain

Definition at line 81 of file BCEngineMCMC.h.

         { return fMCMCNParameters; }

std::vector<int> BCEngineMCMC::MCMCGetNTrialsFalse (  )  [inline]
Returns:
number of not-accepted trials for each chain

Definition at line 152 of file BCEngineMCMC.h.

         { return fMCMCNTrialsFalse; }

std::vector<int> BCEngineMCMC::MCMCGetNTrialsTrue (  )  [inline]
Returns:
number of accepted trials for each chain

Definition at line 147 of file BCEngineMCMC.h.

         { return fMCMCNTrialsTrue; }

int BCEngineMCMC::MCMCGetPhase (  )  [inline]
Returns:
pointer to the phase of a run.

Definition at line 211 of file BCEngineMCMC.h.

         { return fMCMCPhase; }

std::vector<double> BCEngineMCMC::MCMCGetprobMean (  )  [inline]
Returns:
mean value of the probability for each chain up to the current iteration

Definition at line 158 of file BCEngineMCMC.h.

         { return fMCMCprobMean; }

bool BCEngineMCMC::MCMCGetProposalPointMetropolis ( int  chain,
std::vector< double > &  x 
)

Returns a trial point for the Metropolis algorithm.

Parameters:
chain chain index
x proposal point
Returns:
flag indicating whether the new point lies within the allowed range

Definition at line 584 of file BCEngineMCMC.cxx.

{
   // get unscaled random point. this point might not be in the correct volume.
   MCMCTrialFunction(chain, x);

   // get a proposal point from the trial function and scale it
   for (int i = 0; i < fMCMCNParameters; ++i)
      x[i] = fMCMCx[chain * fMCMCNParameters + i] + x[i] * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));

   // check if the point is in the correct volume.
   for (int i = 0; i < fMCMCNParameters; ++i)
      if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
         return false;

   return true;
}

bool BCEngineMCMC::MCMCGetProposalPointMetropolis ( int  chain,
int  parameter,
std::vector< double > &  x 
)

Returns a trial point for the Metropolis algorithm.

Parameters:
chain chain index
x proposal point
Returns:
flag indicating whether the new point lies within the allowed range

Definition at line 602 of file BCEngineMCMC.cxx.

{
   // get unscaled random point in the dimension of the chosen
   // parameter. this point might not be in the correct volume.
   double proposal = MCMCTrialFunctionSingle(ichain, ipar);

   // copy the old point into the new
   for (int i = 0; i < fMCMCNParameters; ++i)
      x[i] = fMCMCx[ichain * fMCMCNParameters + i];

   // modify the parameter under study
   x[ipar] += proposal * (fMCMCBoundaryMax[ipar] - fMCMCBoundaryMin[ipar]);

   // check if the point is in the correct volume.
   if ((x[ipar] < fMCMCBoundaryMin[ipar]) || (x[ipar] > fMCMCBoundaryMax[ipar]))
      return false;

   return true;
}

double BCEngineMCMC::MCMCGetRValue (  )  [inline]
Returns:
R-value

Definition at line 251 of file BCEngineMCMC.h.

         { return fMCMCRValue; }

double BCEngineMCMC::MCMCGetRValueCriterion (  )  [inline]
Returns:
R-value criterion

Definition at line 241 of file BCEngineMCMC.h.

         { return fMCMCRValueCriterion; }

double BCEngineMCMC::MCMCGetRValueParameters ( int  i  )  [inline]
Returns:
R-value for a parameter
Parameters:
i parameter index

Definition at line 257 of file BCEngineMCMC.h.

         { return fMCMCRValueParameters.at(i); }

double BCEngineMCMC::MCMCGetRValueParametersCriterion (  )  [inline]
Returns:
R-value criterion for parameters

Definition at line 246 of file BCEngineMCMC.h.

bool BCEngineMCMC::MCMCGetRValueStrict (  )  [inline]

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

Definition at line 261 of file BCEngineMCMC.h.

         { return fMCMCRValueUseStrict; }

TRandom3* BCEngineMCMC::MCMCGetTRandom3 (  )  [inline]

Return the random number generator. Any non-zero seed gives reproducible behavior,e.g. m->MCMCGetTRandom3()->SetSeed(21340)

Definition at line 293 of file BCEngineMCMC.h.

         { return fRandom; }

std::vector<double> BCEngineMCMC::MCMCGetTrialFunctionScaleFactor (  )  [inline]
Returns:
scale factor for all parameters and chains

Definition at line 169 of file BCEngineMCMC.h.

std::vector< double > BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( int  ichain  ) 
Returns:
scale factor for all parameters and achain.
Parameters:
ichain chain index

Definition at line 509 of file BCEngineMCMC.cxx.

{
   // create a new vector with the length of fMCMCNParameters
   std::vector<double> x;

   // check if ichain is in range
   if (ichain < 0 || ichain >= fMCMCNChains)
      return x;

   // copy the scale factors into the temporary vector
   for (int j = 0; j < fMCMCNParameters; ++j)
      x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNParameters + j));

   return x;
}

double BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( int  ichain,
int  ipar 
)
Returns:
scale factor for a parameter and a chain.
Parameters:
ichain chain index
ipar parameter index

Definition at line 526 of file BCEngineMCMC.cxx.

{
   // check if ichain is in range
   if (ichain < 0 || ichain >= fMCMCNChains)
      return 0;

   // check if ipar is in range
   if (ipar < 0 || ipar >= fMCMCNParameters)
      return 0;

   // return component of ipar point in the ichain chain
   return fMCMCTrialFunctionScaleFactor.at(ichain *  fMCMCNChains + ipar);
}

std::vector<double> BCEngineMCMC::MCMCGetVariance (  )  [inline]
Returns:
mean value of the probability for each chain up to the current iteration

Definition at line 164 of file BCEngineMCMC.h.

         { return fMCMCprobVar; }

std::vector<double> BCEngineMCMC::MCMCGetx (  )  [inline]
Returns:
current point of each Markov chain

Definition at line 185 of file BCEngineMCMC.h.

         { return fMCMCx; }

std::vector< double > BCEngineMCMC::MCMCGetx ( int  ichain  ) 
Parameters:
ichain index of the Markov chain
Returns:
current point of the Markov chain

Definition at line 541 of file BCEngineMCMC.cxx.

{
   // create a new vector with the length of fMCMCNParameters
   std::vector<double> x;

   // check if ichain is in range
   if (ichain < 0 || ichain >= fMCMCNChains)
      return x;

   // copy the point in the ichain chain into the temporary vector
   for (int j = 0; j < fMCMCNParameters; ++j)
      x.push_back(fMCMCx.at(ichain * fMCMCNParameters + j));

   return x;
}

double BCEngineMCMC::MCMCGetx ( int  ichain,
int  ipar 
)
Parameters:
ichain chain index
ipar parameter index
Returns:
parameter of the Markov chain

Definition at line 558 of file BCEngineMCMC.cxx.

{
   // check if ichain is in range
   if (ichain < 0 || ichain >= fMCMCNChains)
      return 0;

   // check if ipar is in range
   if (ipar < 0 || ipar >= fMCMCNParameters)
      return 0;

   // return component of jth point in the ith chain
   return fMCMCx.at(ichain *  fMCMCNParameters + ipar);
}

void BCEngineMCMC::MCMCInChainCheckMaximum (  ) 

Updates statistics: find new maximum

Definition at line 767 of file BCEngineMCMC.cxx.

{
   // loop over all chains
   for (int i = 0; i < fMCMCNChains; ++i)
   {
      // check if new maximum is found or chain is at the beginning
      if (fMCMCprob[i] > fMCMCprobMax[i] || fMCMCNIterations[i] == 1)
      {
         // copy maximum value
         fMCMCprobMax[i] = fMCMCprob[i];

         // copy mode of chain
         for (int j = 0; j < fMCMCNParameters; ++j)
            fMCMCxMax[i * fMCMCNParameters + j] = fMCMCx[i * fMCMCNParameters + j];
      }
   }
}

void BCEngineMCMC::MCMCInChainFillHistograms (  ) 

Updates statistics: fill marginalized distributions

Definition at line 819 of file BCEngineMCMC.cxx.

{
   // check if histograms are supposed to be filled
   if (!fMCMCFlagFillHistograms)
      return;

   // loop over chains
   for (int i = 0; i < fMCMCNChains; ++i)
   {
      // fill each 1-dimensional histogram (if supposed to be filled)
      for (int j = 0; j < fMCMCNParameters; ++j)
         if (fMCMCFlagsFillHistograms.at(j))
            fMCMCH1Marginalized[j]->Fill(fMCMCx[i * fMCMCNParameters + j]);

      // fill each 2-dimensional histogram (if supposed to be filled)
      int counter = 0;

      for (int j = 0; j < fMCMCNParameters; ++j)
         for (int k = 0; k < j; ++k)
         {
            if (fMCMCFlagsFillHistograms.at(j) && fMCMCFlagsFillHistograms.at(k))
               fMCMCH2Marginalized[counter]->Fill(fMCMCx[i*fMCMCNParameters+k],fMCMCx[i* fMCMCNParameters+j]);
            counter ++;
         }
   }
}

void BCEngineMCMC::MCMCInChainTestConvergenceAllChains (  ) 

Updates statistics: check convergence

Definition at line 847 of file BCEngineMCMC.cxx.

{
   // calculate number of entries in this part of the chain
   int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0];

   if (fMCMCNChains > 1 && npoints > 1)
   {
      // define flag for convergence
      bool flag_convergence = true;

      // loop over parameters
      for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters){
         // extract means and variances
         std::vector<double> means(fMCMCNChains);
         std::vector<double> variances(fMCMCNChains);

         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains) {
            // get parameter index
            int index = ichains * fMCMCNParameters + iparameters;
            means[ichains] = fMCMCxMean[index];
            variances[ichains] = fMCMCxVar[index];
         }

         fMCMCRValueParameters[iparameters] = BCMath::Rvalue(means, variances, npoints, fMCMCRValueUseStrict);

         // set flag to false if convergence criterion is not fulfilled for the parameter
         if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
            flag_convergence = false;

         // else: leave convergence flag true for that parameter
      }

      fMCMCRValue = BCMath::Rvalue(fMCMCprobMean, fMCMCprobVar, npoints, fMCMCRValueUseStrict);

      // set flag to false if convergence criterion is not fulfilled for the log-likelihood
      if (!((fMCMCRValue - 1.0) < fMCMCRValueCriterion))
         flag_convergence = false;

      // remember number of iterations needed to converge
      if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
         fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters;
   }
}

void BCEngineMCMC::MCMCInChainUpdateStatistics (  ) 

Updates statistics:

Definition at line 786 of file BCEngineMCMC.cxx.

{
   // calculate number of entries in this part of the chain
   int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0];

   // length of vectors
   int nentries = fMCMCNParameters * fMCMCNChains;

   // loop over all parameters of all chains
   for (int i = 0; i < nentries; ++i) {
      // calculate mean value of each parameter in the chain for this part
      fMCMCxMean[i] += (fMCMCx[i] - fMCMCxMean[i]) / double(npoints);

      // calculate variance of each chain for this part
      if (npoints > 1)
         fMCMCxVar[i] = (1.0 - 1./double(npoints)) * fMCMCxVar[i]
            + (fMCMCx[i] - fMCMCxMean[i]) * (fMCMCx[i] - fMCMCxMean[i]) / double(npoints - 1);
   }

   // loop over chains
   for (int i = 0; i < fMCMCNChains; ++i) {
      // calculate mean value of each chain for this part
      fMCMCprobMean[i] += (fMCMCprob[i] - fMCMCprobMean[i]) / double(npoints);

      // calculate variance of each chain for this part
      if (npoints > 1)
         fMCMCprobVar[i] = (1.0 - 1/double(npoints)) * fMCMCprobVar[i]
            + (fMCMCprob[i] - fMCMCprobMean[i]) * (fMCMCprob[i] - fMCMCprobMean[i]) / double(npoints - 1);
   }

}

void BCEngineMCMC::MCMCInChainWriteChains (  ) 

Updates statistics: write chains to file

Definition at line 892 of file BCEngineMCMC.cxx.

{
   // loop over all chains
   for (int i = 0; i < fMCMCNChains; ++i)
      fMCMCTrees[i]->Fill();
}

int BCEngineMCMC::MCMCInitialize (  ) 

Initializes the engine.

Returns:
An error code

Definition at line 1548 of file BCEngineMCMC.cxx.

{
   // reset values
   MCMCResetResults();

   // free memory for vectors
   fMCMCNIterations.assign(fMCMCNChains, 0);
   fMCMCprobMean.assign(fMCMCNChains, 0);
   fMCMCprobVar.assign(fMCMCNChains, 0);
   fMCMCprob.assign(fMCMCNChains, -1.0);
   fMCMCprobMax.assign(fMCMCNChains, -1.0);

   fMCMCNTrialsTrue.assign(fMCMCNChains * fMCMCNParameters, 0);
   fMCMCNTrialsFalse.assign(fMCMCNChains * fMCMCNParameters, 0);
   fMCMCxMax.assign(fMCMCNChains * fMCMCNParameters, 0.);
   fMCMCxMean.assign(fMCMCNChains * fMCMCNParameters, 0);
   fMCMCxVar.assign(fMCMCNChains * fMCMCNParameters, 0);

   fMCMCRValueParameters.assign(fMCMCNParameters, 0.);

   if (fMCMCTrialFunctionScaleFactorStart.size() == 0)
      fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fMCMCNParameters, 1.0);
   else
      for (int i = 0; i < fMCMCNChains; ++i)
         for (int j = 0; j < fMCMCNParameters; ++j)
            fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j));

   // set initial position
   if (fMCMCFlagInitialPosition == 2) // user defined points
   {
      // define flag
      bool flag = true;

      // check the length of the array of initial positions
      if (int(fMCMCInitialPosition.size()) != (fMCMCNChains * fMCMCNParameters))
      {
         BCLog::OutError("BCEngine::MCMCInitialize : Length of vector containing initial positions does not have required length.");
         flag = false;
      }

      // check the boundaries
      if (flag)
      {
         for (int j = 0; j < fMCMCNChains; ++j)
            for (int i = 0; i < fMCMCNParameters; ++i)
               if (fMCMCInitialPosition[j * fMCMCNParameters + i] < fMCMCBoundaryMin[i] ||
                     fMCMCInitialPosition[j * fMCMCNParameters + i] > fMCMCBoundaryMax[i])
               {
                  BCLog::OutError("BCEngine::MCMCInitialize : Initial position out of boundaries.");
                  flag = false;
               }
      }

      // check flag
      if (!flag)
         fMCMCFlagInitialPosition = 1;
   }

   if (fMCMCFlagInitialPosition == 0) // center of the interval
      for (int j = 0; j < fMCMCNChains; ++j)
         for (int i = 0; i < fMCMCNParameters; ++i)
            fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));

   else if (fMCMCFlagInitialPosition == 2) // user defined
   {
      for (int j = 0; j < fMCMCNChains; ++j)
         for (int i = 0; i < fMCMCNParameters; ++i)
            fMCMCx.push_back(fMCMCInitialPosition.at(j * fMCMCNParameters + i));
   }

   else
      for (int j = 0; j < fMCMCNChains; ++j) // random number (default)
         for (int i = 0; i < fMCMCNParameters; ++i)
            fMCMCx.push_back(fMCMCBoundaryMin[i] + fRandom->Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));

   // copy the point of the first chain
   fMCMCxLocal.clear();
   for (int i = 0; i < fMCMCNParameters; ++i)
      fMCMCxLocal.push_back(fMCMCx[i]);

   // define 1-dimensional histograms for marginalization
   for(int i = 0; i < fMCMCNParameters; ++i)
   {
      TH1D * h1 = 0;
      if (fMCMCFlagsFillHistograms.at(i))
         h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i), "",
               fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i]);
      fMCMCH1Marginalized.push_back(h1);
   }

   for(int i = 0; i < fMCMCNParameters; ++i)
      for (int k = 0; k < i; ++k)
      {
         TH2D * h2 = 0;
         if (fMCMCFlagsFillHistograms.at(i) && fMCMCFlagsFillHistograms.at(k))
            h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, k), "",
                  fMCMCH1NBins[k], fMCMCBoundaryMin[k], fMCMCBoundaryMax[k],
                  fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i] );
         fMCMCH2Marginalized.push_back(h2);
      }

   return 1;
}

void BCEngineMCMC::MCMCInitializeMarkovChains (  ) 

Initializes Markov chains.

Definition at line 1491 of file BCEngineMCMC.cxx.

{
   // evaluate function at the starting point
   std::vector<double> x0;

   for (int j = 0; j < fMCMCNChains; ++j)
   {
      x0.clear();
      for (int i = 0; i < fMCMCNParameters; ++i)
         x0.push_back(fMCMCx[j * fMCMCNParameters + i]);
      fMCMCprob[j] = LogEval(x0);
   }

   x0.clear();
}

void BCEngineMCMC::MCMCInitializeMarkovChainTrees (  ) 

Initialize trees containing the Markov chains.

Definition at line 467 of file BCEngineMCMC.cxx.

{
// clear vector
   fMCMCTrees.clear();

// create new trees
   for (int i = 0; i < fMCMCNChains; ++i) {
      fMCMCTrees.push_back(new TTree(TString::Format("MarkovChainTree_%i", i), "MarkovChainTree"));
      fMCMCTrees[i]->Branch("Iteration",       &fMCMCNIterations[i],  "iteration/I");
      fMCMCTrees[i]->Branch("NParameters",     &fMCMCNParameters,     "parameters/I");
      fMCMCTrees[i]->Branch("LogProbability",  &fMCMCprob[i],         "log(probability)/D");
      fMCMCTrees[i]->Branch("Phase",           &fMCMCPhase,           "phase/I");
      fMCMCTrees[i]->Branch("Cycle",           &fMCMCCycle,           "cycle/I");

      for (int j = 0; j < fMCMCNParameters; ++j)
         fMCMCTrees[i]->Branch(TString::Format("Parameter%i", j),
               &fMCMCx[i * fMCMCNParameters + j],
               TString::Format("parameter %i/D", j));
   }
}

virtual void BCEngineMCMC::MCMCIterationInterface (  )  [inline, 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 BCIntegrate, and BCRooInterface.

Definition at line 553 of file BCEngineMCMC.h.

         {}

int BCEngineMCMC::MCMCMetropolis (  ) 

Runs Metropolis algorithm.

Definition at line 1296 of file BCEngineMCMC.cxx.

{
   // check if prerun should be performed
   if (fMCMCFlagPreRun)
      MCMCMetropolisPreRun();
   else
      BCLog::OutWarning("BCEngineMCMC::MCMCMetropolis. Not running prerun. This can cause trouble if the data have changed.");

   // print to screen
   BCLog::OutSummary( "Run Metropolis MCMC...");

   // reset run statistics
   MCMCResetRunStatistics();

   // set phase and cycle number
   fMCMCPhase = 2;
   fMCMCCycle = 0;

   // perform run
   BCLog::OutSummary(Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));

//   int counterupdate = 0;
//   bool convergence = false;
//   bool flagefficiency = false;

//   std::vector<double> efficiency;

//   for (int i = 0; i < fMCMCNParameters; ++i)
//      for (int j = 0; j < fMCMCNChains; ++j)
//         efficiency.push_back(0.0);

   int nwrite = fMCMCNIterationsRun/10;
   if(nwrite < 100)
      nwrite=100;
   else if(nwrite < 500)
      nwrite=1000;
   else if(nwrite < 10000)
      nwrite=1000;
   else
      nwrite=10000;

   // start the run
   for (fMCMCCurrentIteration = 1; fMCMCCurrentIteration <= fMCMCNIterationsRun; ++fMCMCCurrentIteration)
   {
      if ( (fMCMCCurrentIteration)%nwrite == 0 )
         BCLog::OutDetail(Form(" --> iteration number %i (%.2f%%)", fMCMCCurrentIteration, (double)(fMCMCCurrentIteration)/(double)fMCMCNIterationsRun*100.));

      // if the flag is set then run over the parameters one after the other.
      if (fMCMCFlagOrderParameters)
      {
         // loop over parameters
         for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
         {
            // loop over chains
            for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
               MCMCGetNewPointMetropolis(ichains, iparameters);

            // reset current chain
            fMCMCCurrentChain = -1;

            // update search for maximum
            MCMCInChainCheckMaximum();

         } // end loop over all parameters

         // check if the current iteration is consistent with the lag
         if ( fMCMCCurrentIteration % fMCMCNLag == 0)
         {
            // fill histograms
            MCMCInChainFillHistograms();

            // write chain to file
            if (fMCMCFlagWriteChainToFile)
               MCMCInChainWriteChains();

            // do anything interface
            MCMCIterationInterface();
         }
      }
      // if the flag is not set then run over the parameters at the same time.
      else
      {
         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
            // get new point
            MCMCGetNewPointMetropolis(ichains);

         // reset current chain
         fMCMCCurrentChain = -1;

         // update search for maximum
         MCMCInChainCheckMaximum();

         // check if the current iteration is consistent with the lag
         if (fMCMCCurrentIteration % fMCMCNLag == 0)
         {
            // fill histograms
            MCMCInChainFillHistograms();

            // write chain to file
            if (fMCMCFlagWriteChainToFile)
               MCMCInChainWriteChains();

            // do anything interface
            MCMCIterationInterface();
         }
      }

   } // end run

   // print convergence status
   BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));

   // print modes

   // find global maximum
   double probmax = fMCMCprobMax.at(0);
   int probmaxindex = 0;

   // loop over all chains and find the maximum point
   for (int i = 1; i < fMCMCNChains; ++i) {
      if (fMCMCprobMax.at(i) > probmax)
      {
         probmax = fMCMCprobMax.at(i);
         probmaxindex = i;
      }
   }

   BCLog::OutDetail(" --> Global mode from MCMC:");
   BCLog::OutDebug(Form(" --> Posterior value: %g", probmax));
   int ndigits = (int) log10(fMCMCNParameters);
   for (int i = 0; i < fMCMCNParameters; ++i)
      BCLog::OutDetail(Form( TString::Format(" -->      parameter %%%di:   %%.4g", ndigits+1),
            i, fMCMCxMax[probmaxindex * fMCMCNParameters + i]));

   // reset counter
   fMCMCCurrentIteration = -1;

   // reset current chain
   fMCMCCurrentChain = -1;

   // set flags
   fMCMCFlagRun = true;

   return 1;
}

int BCEngineMCMC::MCMCMetropolisPreRun (  ) 

Runs a pre run for the Metropolis algorithm.

Definition at line 908 of file BCEngineMCMC.cxx.

{
   // print on screen
   BCLog::OutSummary("Pre-run Metropolis MCMC...");

   // initialize Markov chain
   MCMCInitialize();
   MCMCInitializeMarkovChains();

   // helper variable containing number of digits in the number of parameters
   int ndigits = (int)log10(fMCMCNParameters) +1;
   if(ndigits<4)
      ndigits=4;

   // reset run statistics
   MCMCResetRunStatistics();
   fMCMCNIterationsConvergenceGlobal = -1;

   // perform run
   BCLog::OutSummary(Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax));


   // don't write to file during pre run
   bool tempflag_writetofile = fMCMCFlagWriteChainToFile;
   fMCMCFlagWriteChainToFile = false;

   // initialize counter variables and flags
   fMCMCCurrentIteration = 1;   // counts the number of iterations
   int counterupdate = 1;        // after how many iterations is an update needed?
   bool convergence = false;     // convergence reached?
   bool flagefficiency = false;  // efficiency reached?

   // array of efficiencies
   std::vector<double> efficiency;
   efficiency.assign(fMCMCNParameters * fMCMCNChains, 0.0);

   // how often to check convergence and efficiencies?
   // it's either every fMCMCNParameters*nMCMCNIterationsUpdate (for 5 parameters the default would be 5000)
   // or it's fMCMCNIterationsUpdateMax (10000 by default)
   // whichever of the two is smaller
   int updateLimit = ( fMCMCNIterationsUpdateMax<fMCMCNIterationsUpdate*(fMCMCNParameters)  && fMCMCNIterationsUpdateMax>0 ) ?
         fMCMCNIterationsUpdateMax : fMCMCNIterationsUpdate*(fMCMCNParameters);

   // loop over chains
   for (int ichains = 0; ichains < fMCMCNChains; ++ichains) {
      // loop over parameters
      for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter){
         // global index of the parameter (throughout all the chains)
         int index = ichains * fMCMCNParameters + iparameter;
         // reset counters
         fMCMCNTrialsTrue[index] = 0;
         fMCMCNTrialsFalse[index] = 0;
         fMCMCxMean[index] = fMCMCx[index];
         fMCMCxVar[index] = 0; 
      }
      fMCMCprobMean[ichains] = fMCMCprob[ichains];
      fMCMCprobVar[ichains] = 0;
   }

   // set phase and cycle number
   fMCMCPhase = 1; 
   fMCMCCycle = 0;

   // run chain ...
   // (a) for at least a minimum number of iterations,
   // (b) until a maximum number of iterations is reached,
   // (c) or until convergence is reached and the efficiency is in the
   //     specified region
   while (fMCMCCurrentIteration < fMCMCNIterationsPreRunMin || (fMCMCCurrentIteration < fMCMCNIterationsMax && !(convergence && flagefficiency)))
   {
      //-------------------------------------------
      // reset flags and counters
      //-------------------------------------------

      // set convergence to false by default
      convergence = false;

      // set number of iterations needed to converge to negative
      fMCMCNIterationsConvergenceGlobal = -1;

      //-------------------------------------------
      // get new point in n-dim space
      //-------------------------------------------

      // if the flag is set then run over the parameters one after the other.
      if (fMCMCFlagOrderParameters)
      {
         // loop over parameters
         for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
         {
            // loop over chains
            for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
               MCMCGetNewPointMetropolis(ichains, iparameters);

            // search for global maximum
            MCMCInChainCheckMaximum();
         }
      }

      // if the flag is not set then run over the parameters at the same time.
      else
      {
         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
            // get new point
            MCMCGetNewPointMetropolis(ichains);

         // search for global maximum
         MCMCInChainCheckMaximum();
      }

      //-------------------------------------------
      // print out message to log
      //-------------------------------------------

      // progress printout
      if ( fMCMCCurrentIteration > 0 && fMCMCCurrentIteration % fMCMCNIterationsUpdate == 0 )
         BCLog::OutDetail(Form(" --> Iteration %i", fMCMCNIterations[0]/fMCMCNParameters));

      //-------------------------------------------
      // update statistics
      //-------------------------------------------

      if (counterupdate > 1)
         MCMCInChainUpdateStatistics();

      //-------------------------------------------
      // update scale factors and check convergence
      //-------------------------------------------

      // debugKK
      // check if this line makes sense
      if ( counterupdate % updateLimit == 0 && counterupdate > 0 && fMCMCCurrentIteration >= fMCMCNIterationsPreRunMin)
//      if ( fMCMCCurrentIteration % fMCMCNIterationsUpdate == 0 && counterupdate > 1 && fMCMCCurrentIteration >= fMCMCNIterationsPreRunMin)
      {
         // -----------------------------
         // reset flags and counters
         // -----------------------------

         bool rvalues_ok = true;

         static bool has_converged = false;

         // reset the number of iterations needed for convergence to
         // negative
         fMCMCNIterationsConvergenceGlobal = -1;

         // -----------------------------
         // check convergence status
         // -----------------------------

         // test convergence
         MCMCInChainTestConvergenceAllChains();

         // set convergence flag
         if (fMCMCNIterationsConvergenceGlobal > 0)
            convergence = true;

         // print convergence status:
         if (convergence && fMCMCNChains > 1)
            BCLog::OutDetail(Form("     * Convergence status: Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
         else if (!convergence && fMCMCNChains > 1)
         {
            BCLog::OutDetail(Form("     * Convergence status: Set of %i Markov chains did not converge after %i iterations.", fMCMCNChains, fMCMCCurrentIteration));

            BCLog::OutDetail("       - R-Values:");
            for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
            {
               if(fabs(fMCMCRValueParameters[iparameter]-1.) < fMCMCRValueParametersCriterion)
                  BCLog::OutDetail(Form( TString::Format("         parameter %%%di :  %%.06f",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
               else
               {
                  BCLog::OutDetail(Form( TString::Format("         parameter %%%di :  %%.06f <--",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
                  rvalues_ok = false;
               }
            }
            if(fabs(fMCMCRValue-1.) < fMCMCRValueCriterion)
               BCLog::OutDetail(Form("         log-likelihood :  %.06f", fMCMCRValue));
            else
            {
               BCLog::OutDetail(Form("         log-likelihood :  %.06f <--", fMCMCRValue));
               rvalues_ok = false;
            }
         }

         // set convergence flag
         if(!has_converged)
            if(rvalues_ok)
               has_converged = true;

         // -----------------------------
         // check efficiency status
         // -----------------------------

         // set flag
         flagefficiency = true;

         bool flagprintefficiency = true;

         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
         {
            // loop over parameters
            for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
            {
               // global index of the parameter (throughout all the chains)
               int index = ichains * fMCMCNParameters + iparameter;

               // calculate efficiency
               efficiency[index] = double(fMCMCNTrialsTrue[index]) / double(fMCMCNTrialsTrue[index] + fMCMCNTrialsFalse[index]);

               // adjust scale factors if efficiency is too low
               if (efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01)
               {
                  if (flagprintefficiency)
                  {
                     BCLog::OutDetail(Form("     * Efficiency status: Efficiencies not within pre-defined range."));
                     BCLog::OutDetail(Form("       - Efficiencies:"));
                     flagprintefficiency = false;
                  }

                  double fscale=2.;
                  if(has_converged && fMCMCEfficiencyMin/efficiency[index] > 2.)
                     fscale = 4.;
                  fMCMCTrialFunctionScaleFactor[index] /= fscale;

                  BCLog::OutDetail(Form("         Efficiency of parameter %i dropped below %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
                        iparameter, 100. * fMCMCEfficiencyMin, 100. * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index]));
               }

               // adjust scale factors if efficiency is too high
               else if (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.0)
               {
                  if (flagprintefficiency)
                  {
                     BCLog::OutDetail(Form("   * Efficiency status: Efficiencies not within pre-defined ranges."));
                     BCLog::OutDetail(Form("     - Efficiencies:"));
                     flagprintefficiency = false;
                  }

                  fMCMCTrialFunctionScaleFactor[index] *= 2.;

                  BCLog::OutDetail(Form("         Efficiency of parameter %i above %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
                        iparameter, 100.0 * fMCMCEfficiencyMax, 100.0 * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index]));
               }

               // check flag
               if ((efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01)
                     || (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.))
                  flagefficiency = false;
            } // end of running over all parameters
         } // end of running over all chains

         // print to screen
         if (flagefficiency)
            BCLog::OutDetail(Form("     * Efficiency status: Efficiencies within pre-defined ranges."));

         // -----------------------------
         // reset counters
         // -----------------------------

         counterupdate = 0;

         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains) {
            // loop over parameters
            for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter){
               // global index of the parameter (throughout all the chains)
               int index = ichains * fMCMCNParameters + iparameter;
               // reset counters
               fMCMCNTrialsTrue[index] = 0;
               fMCMCNTrialsFalse[index] = 0;
               fMCMCxMean[index] = fMCMCx[index];
               fMCMCxVar[index] = 0;
            }
            fMCMCprobMean[ichains] = fMCMCprob[ichains];
            fMCMCprobVar[ichains] = 0;
         }
      } // end if update scale factors and check convergence

      //-------------------------------------------
      // write chain to file
      //-------------------------------------------

      // write chain to file
      if (fMCMCFlagWritePreRunToFile)
         MCMCInChainWriteChains();

      //-------------------------------------------
      // increase counters
      //-------------------------------------------

      if (counterupdate == 0 && fMCMCCurrentIteration > 1)
        fMCMCCycle++;
      fMCMCCurrentIteration++;
      counterupdate++;

   } // end of running

   // decrease counter by one since it didn't really run that long
   //   fMCMCCurrentIteration--;
   //   counterupdate--;

   // did we check convergence at least once ?
   if (fMCMCCurrentIteration<updateLimit)
   {
      BCLog::OutWarning(" Convergence never checked !");
      BCLog::OutWarning("   Increase maximum number of iterations in the pre-run /MCMCSetNIterationsMax()/");
      BCLog::OutWarning("   or decrease maximum number of iterations for update  /MCMCSetNIterationsUpdateMax()/");
   }

   // ---------------
   // after chain run
   // ---------------

   // define convergence status
   if (fMCMCNIterationsConvergenceGlobal > 0)
      fMCMCFlagConvergenceGlobal = true;
   else
      fMCMCFlagConvergenceGlobal = false;

   // print convergence status
   if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && !flagefficiency)
      BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations but could not adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));

   else if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && flagefficiency)
      BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations and all scales are adjusted.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));

   else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && flagefficiency)
      BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));

   else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && !flagefficiency)
      BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations and could not adjust scales.", fMCMCNChains, fMCMCNIterationsMax));

   else if(fMCMCNChains == 1)
      BCLog::OutSummary(" --> No convergence criterion for a single chain defined.");

   else
      BCLog::OutSummary(" --> Only one Markov chain. No global convergence criterion defined.");

   BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCCurrentIteration));


   // print efficiencies
   std::vector<double> efficiencies;

   for (int i = 0; i < fMCMCNParameters; ++i)
      efficiencies.push_back(0.);

   BCLog::OutDetail(" --> Average efficiencies:");
   for (int i = 0; i < fMCMCNParameters; ++i)
   {
      for (int j = 0; j < fMCMCNChains; ++j)
         efficiencies[i] += efficiency[j * fMCMCNParameters + i] / double(fMCMCNChains);

      BCLog::OutDetail(Form(TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * efficiencies[i]));
   }


   // print scale factors
   std::vector<double> scalefactors;

   for (int i = 0; i < fMCMCNParameters; ++i)
      scalefactors.push_back(0.0);

   BCLog::OutDetail(" --> Average scale factors:");
   for (int i = 0; i < fMCMCNParameters; ++i)
   {
      for (int j = 0; j < fMCMCNChains; ++j)
         scalefactors[i] += fMCMCTrialFunctionScaleFactor[j * fMCMCNParameters + i] / double(fMCMCNChains);

      BCLog::OutDetail(Form( TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * scalefactors[i]));
   }

   // reset flag
   fMCMCFlagWriteChainToFile = tempflag_writetofile;

   // reset current iteration
   fMCMCCurrentIteration = -1;

   // reset current chain
   fMCMCCurrentChain = -1;

   // no error
   return 1;
}

int BCEngineMCMC::MCMCResetResults (  ) 

Reset the MCMC variables.

Returns:
An error code

Definition at line 1508 of file BCEngineMCMC.cxx.

{
   // reset variables
   fMCMCNIterations.clear();
   fMCMCNTrialsTrue.clear();
   fMCMCNTrialsFalse.clear();
   fMCMCTrialFunctionScaleFactor.clear();
   fMCMCprobMean.clear();
   fMCMCprobVar.clear();
   fMCMCxMean.clear();
   fMCMCxVar.clear();
   fMCMCx.clear();
   fMCMCprob.clear();
   fMCMCxMax.clear();
   fMCMCprobMax.clear();
   fMCMCNIterationsConvergenceGlobal = -1;
   fMCMCRValueParameters.clear();

   for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
      if (fMCMCH1Marginalized[i])
         delete fMCMCH1Marginalized[i];

   for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
      if (fMCMCH2Marginalized[i])
         delete fMCMCH2Marginalized[i];

   // clear plots
   fMCMCH1Marginalized.clear();
   fMCMCH2Marginalized.clear();

   // reset flags
   fMCMCFlagPreRun = true;
   fMCMCFlagRun = false;
   fMCMCFlagConvergenceGlobal = false;

   // no errors
   return 1;
}

void BCEngineMCMC::MCMCResetRunStatistics (  ) 

Resets the run statistics.

Definition at line 1444 of file BCEngineMCMC.cxx.

{
   for (int j = 0; j < fMCMCNChains; ++j)
   {
      fMCMCNIterations[j]  = 0;
      fMCMCNTrialsTrue[j]  = 0;
      fMCMCNTrialsFalse[j] = 0;
      fMCMCprobMean[j]         = 0;
      fMCMCprobVar[j]     = 0;

      for (int k = 0; k < fMCMCNParameters; ++k)
      {
         fMCMCNTrialsTrue[j * fMCMCNParameters + k]  = 0;
         fMCMCNTrialsFalse[j * fMCMCNParameters + k] = 0;
      }
   }

   // reset marginalized distributions
   for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
      if (fMCMCH1Marginalized[i])
         fMCMCH1Marginalized[i]->Reset();

   for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
      if (fMCMCH2Marginalized[i])
         fMCMCH2Marginalized[i]->Reset();

   fMCMCRValue = 100;
}

void BCEngineMCMC::MCMCSetFlagFillHistograms ( bool  flag  ) 

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

Definition at line 433 of file BCEngineMCMC.cxx.

{
   fMCMCFlagFillHistograms = flag;

   for (int i = 0; i < fMCMCNParameters; ++i)
      fMCMCFlagsFillHistograms[i] = flag;
}

void BCEngineMCMC::MCMCSetFlagFillHistograms ( int  index,
bool  flag 
)

Sets the flag for a single parameter to either fill histograms or not.

Definition at line 442 of file BCEngineMCMC.cxx.

{
   // check if index is within range
   if (index < 0 || index > fMCMCNParameters)
   {
      BCLog::OutWarning("BCEngineMCMC :MCMCSetFlagFillHistograms. Index out of range.");
      return;
   }

   // set flag
   fMCMCFlagsFillHistograms[index] = flag;
}

void BCEngineMCMC::MCMCSetFlagInitialPosition ( int  flag  )  [inline]

Sets flag which defines initial position.

Definition at line 377 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetFlagOrderParameters ( bool  flag  )  [inline]

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

Definition at line 383 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetFlagPreRun ( bool  flag  )  [inline]

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

Definition at line 393 of file BCEngineMCMC.h.

      { fMCMCFlagPreRun = flag; }

void BCEngineMCMC::MCMCSetInitialPositions ( std::vector< std::vector< double > >  x0s  ) 

Sets the initial positions for all chains.

Parameters:
x0s initial positions for all chains.

Definition at line 419 of file BCEngineMCMC.cxx.

{
   // create new vector
   std::vector<double> y0s;

   // loop over vector elements
   for (int i = 0; i < int(x0s.size()); ++i)
      for (int j = 0; j < int((x0s.at(i)).size()); ++j)
         y0s.push_back((x0s.at(i)).at(j));

   MCMCSetInitialPositions(y0s);
}

void BCEngineMCMC::MCMCSetInitialPositions ( const std::vector< double > &  x0s  ) 

Sets the initial positions for all chains.

Parameters:
x0s initial positions for all chains.

Definition at line 403 of file BCEngineMCMC.cxx.

{
   // clear the existing initial position vector
   fMCMCInitialPosition.clear();

   // copy the initial positions
   int n = int(x0s.size());

   for (int i = 0; i < n; ++i)
      fMCMCInitialPosition.push_back(x0s.at(i));

   // use these intial positions for the Markov chain
   MCMCSetFlagInitialPosition(2);
}

void BCEngineMCMC::MCMCSetMarkovChainTrees ( std::vector< TTree * >  trees  ) 

Sets the tree containing the Markov chains.

Definition at line 456 of file BCEngineMCMC.cxx.

{
// clear vector
   fMCMCTrees.clear();

   // copy tree
   for (int i = 0; i < int(trees.size()); ++i)
      fMCMCTrees.push_back(trees[i]);
}

void BCEngineMCMC::MCMCSetMaximumEfficiency ( double  efficiency  )  [inline]

Sets the maximum efficiency required for a chain.

Definition at line 352 of file BCEngineMCMC.h.

         { fMCMCEfficiencyMax = efficiency; }

void BCEngineMCMC::MCMCSetMinimumEfficiency ( double  efficiency  )  [inline]

Sets the minimum efficiency required for a chain.

Definition at line 347 of file BCEngineMCMC.h.

         { fMCMCEfficiencyMin = efficiency; }

void BCEngineMCMC::MCMCSetNChains ( int  n  ) 

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

Definition at line 394 of file BCEngineMCMC.cxx.

{
   fMCMCNChains = n;

   // re-initialize
   MCMCInitialize();
}

void BCEngineMCMC::MCMCSetNIterationsMax ( int  n  )  [inline]

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

Definition at line 317 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetNIterationsPreRunMin ( int  n  )  [inline]

Sets the minimum number of iterations in the pre-run

Definition at line 327 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetNIterationsRun ( int  n  )  [inline]

Sets the number of iterations.

Definition at line 322 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetNIterationsUpdate ( int  n  )  [inline]

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

Parameters:
n The number of iterations.

Definition at line 334 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetNIterationsUpdateMax ( int  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:
n maximum number of iterations.

Definition at line 342 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetNLag ( int  n  )  [inline]

Sets the lag of the Markov chains

Definition at line 312 of file BCEngineMCMC.h.

         { fMCMCNLag = n; }

void BCEngineMCMC::MCMCSetPrecision ( BCEngineMCMC::Precision  precision  ) 
void BCEngineMCMC::MCMCSetRValueCriterion ( double  r  )  [inline]

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

Definition at line 398 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetRValueParametersCriterion ( double  r  )  [inline]

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

Definition at line 403 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetRValueStrict ( bool  strict = true  )  [inline]

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

Definition at line 407 of file BCEngineMCMC.h.

      { fMCMCRValueUseStrict = strict; }

void BCEngineMCMC::MCMCSetTrialFunctionScaleFactor ( std::vector< double >  scale  )  [inline]

Set the scale factors for the trial functions

Parameters:
scale a vector of doubles containing the scale factors

Definition at line 303 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetValuesDefault (  ) 

Set the default values for the MCMC chain.

Definition at line 44 of file BCEngineMCMC.cxx.

void BCEngineMCMC::MCMCSetValuesDetail (  ) 
void BCEngineMCMC::MCMCSetValuesQuick (  ) 
void BCEngineMCMC::MCMCSetWriteChainToFile ( bool  flag  )  [inline]

Sets flag to write Markov chains to file.

Definition at line 357 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCSetWritePreRunToFile ( bool  flag  )  [inline]

Sets flag to write pre run to file.

Definition at line 362 of file BCEngineMCMC.h.

void BCEngineMCMC::MCMCTrialFunction ( int  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:
ichain the chain index
x point with the dimension fMCMCNParameters

Definition at line 489 of file BCEngineMCMC.cxx.

{
   // call MCMCTrialFunctionSingle() for all parameters by default
   for (int i = 0; i < fMCMCNParameters; ++i)
      x[i] = MCMCTrialFunctionSingle(ichain, i);
}

double BCEngineMCMC::MCMCTrialFunctionSingle ( int  ichain,
int  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:
ichain the chain index
ipar the parameter index
Returns:
the unscaled proposal point

Definition at line 497 of file BCEngineMCMC.cxx.

{
   // no check of range for performance reasons

   // use uniform distribution
//   return = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fRandom->Rndm());

   // Breit-Wigner width adjustable width
   return fRandom->BreitWigner(0.0, fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter]);
}

BCEngineMCMC & BCEngineMCMC::operator= ( const BCEngineMCMC engineMCMC  ) 

Defaut assignment operator

Definition at line 259 of file BCEngineMCMC.cxx.

{
   fMCMCPointerToGetProposalPoint = enginemcmc.fMCMCPointerToGetProposalPoint;
   fMCMCNParameters           = enginemcmc.fMCMCNParameters;
   fMCMCBoundaryMin           = enginemcmc.fMCMCBoundaryMin;
   fMCMCBoundaryMax           = enginemcmc.fMCMCBoundaryMax;
   fMCMCFlagsFillHistograms   = enginemcmc.fMCMCFlagsFillHistograms;
   fMCMCNChains               = enginemcmc.fMCMCNChains;
   fMCMCNLag                  = enginemcmc.fMCMCNLag;
   fMCMCNIterations           = enginemcmc.fMCMCNIterations;
   fMCMCCurrentIteration      = enginemcmc.fMCMCCurrentIteration;
   fMCMCCurrentChain          = enginemcmc.fMCMCCurrentChain;
   fMCMCNIterationsUpdate     = enginemcmc.fMCMCNIterationsUpdate;
   fMCMCNIterationsUpdateMax  = enginemcmc.fMCMCNIterationsUpdateMax;
   fMCMCNIterationsConvergenceGlobal = enginemcmc.fMCMCNIterationsConvergenceGlobal;
   fMCMCFlagConvergenceGlobal = enginemcmc.fMCMCFlagConvergenceGlobal;
   fMCMCNIterationsMax        = enginemcmc.fMCMCNIterationsMax;
   fMCMCNIterationsRun        = enginemcmc.fMCMCNIterationsRun;
   fMCMCNIterationsPreRunMin  = enginemcmc.fMCMCNIterationsPreRunMin;
   fMCMCNTrialsTrue           = enginemcmc.fMCMCNTrialsTrue;
   fMCMCNTrialsFalse          = enginemcmc.fMCMCNTrialsFalse;
   fMCMCFlagWriteChainToFile  = enginemcmc.fMCMCFlagWriteChainToFile;
   fMCMCFlagWritePreRunToFile = enginemcmc.fMCMCFlagWritePreRunToFile;
   fMCMCTrialFunctionScaleFactor = enginemcmc.fMCMCTrialFunctionScaleFactor;
   fMCMCTrialFunctionScaleFactorStart = enginemcmc.fMCMCTrialFunctionScaleFactorStart;
   fMCMCFlagPreRun            = enginemcmc.fMCMCFlagPreRun;
   fMCMCFlagRun               = enginemcmc.fMCMCFlagRun;
   fMCMCInitialPosition       = enginemcmc.fMCMCInitialPosition;
   fMCMCEfficiencyMin         = enginemcmc.fMCMCEfficiencyMin;
   fMCMCEfficiencyMax         = enginemcmc.fMCMCEfficiencyMax;
   fMCMCFlagInitialPosition   = enginemcmc.fMCMCFlagInitialPosition;
   fMCMCFlagOrderParameters   = enginemcmc.fMCMCFlagOrderParameters;
   fMCMCFlagFillHistograms    = enginemcmc.fMCMCFlagFillHistograms;
   fMCMCPhase                 = enginemcmc.fMCMCPhase;
   fMCMCCycle                 = enginemcmc.fMCMCCycle;
   fMCMCx                     = enginemcmc.fMCMCx;
   fMCMCxMax                  = enginemcmc.fMCMCxMax;
   fMCMCxMean                 = enginemcmc.fMCMCxMean;
   fMCMCxVar                  = enginemcmc.fMCMCxVar;
   fMCMCxLocal                = enginemcmc.fMCMCxLocal;
   fMCMCprob                  = enginemcmc.fMCMCprob;
   fMCMCprobMax               = enginemcmc.fMCMCprobMax;
   fMCMCprobMean              = enginemcmc.fMCMCprobMean;
   fMCMCprobVar               = enginemcmc.fMCMCprobVar;
   fMCMCRValueUseStrict       = enginemcmc.fMCMCRValueUseStrict;
   fMCMCRValueCriterion       = enginemcmc.fMCMCRValueCriterion ;
   fMCMCRValueParametersCriterion = enginemcmc.fMCMCRValueParametersCriterion;
   fMCMCRValue                = enginemcmc.fMCMCRValue;
   fMCMCRValueParameters      = enginemcmc.fMCMCRValueParameters;
   if (enginemcmc.fRandom)
      fRandom = new TRandom3(*enginemcmc.fRandom);
   else
      fRandom = 0;
   fMCMCH1NBins               = enginemcmc.fMCMCH1NBins;

   for (int i = 0; i < int(enginemcmc.fMCMCH1Marginalized.size()); ++i) {
      if (enginemcmc.fMCMCH1Marginalized.at(i))
         fMCMCH1Marginalized.push_back(new TH1D(*(enginemcmc.fMCMCH1Marginalized.at(i))));
      else
         fMCMCH1Marginalized.push_back(0);
   }

   for (int i = 0; i < int(enginemcmc.fMCMCH2Marginalized.size()); ++i) {
      if (enginemcmc.fMCMCH2Marginalized.at(i))
         fMCMCH2Marginalized.push_back(new TH2D(*(enginemcmc.fMCMCH2Marginalized.at(i))));
      else
         fMCMCH2Marginalized.push_back(0);
   }

   for (int i = 0; i < int(enginemcmc.fMCMCTrees.size()); ++i) {
      // debugKK
      fMCMCTrees.push_back(0);
   }

   // return this
   return *this;
}

int BCEngineMCMC::SetMarginalized ( int  index,
TH1D *  h 
)

Sets the histogram with 1D marginalized distributions for parameter.

Parameters:
i index of the parameter
h pointer to an existing histogram

Definition at line 1653 of file BCEngineMCMC.cxx.

{
   if((int)fMCMCH1Marginalized.size()<=index)
      return 0;

   if(h==0)
      return 0;

   if((int)fMCMCH1Marginalized.size()==index)
      fMCMCH1Marginalized.push_back(h);
   else
      fMCMCH1Marginalized[index]=h;

   return index;
}

int BCEngineMCMC::SetMarginalized ( int  index1,
int  index2,
TH2D *  h 
)

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

Parameters:
index1 index of the first parameter
index2 index of the second parameter
h pointer to an existing histogram

Definition at line 1670 of file BCEngineMCMC.cxx.

{
   int counter = 0;
   int index = 0;

   // search for correct combination
   for(int i = 0; i < fMCMCNParameters; i++)
      for (int j = 0; j < i; ++j)
      {
         if(j == index1 && i == index2)
            index = counter;
         counter++;
      }

   if((int)fMCMCH2Marginalized.size()<=index)
      return 0;

   if(h==0)
      return 0;

   if((int)fMCMCH2Marginalized.size()==index)
      fMCMCH2Marginalized.push_back(h);
   else
      fMCMCH2Marginalized[index]=h;

   return index;
}


Member Data Documentation

std::vector<double> BCEngineMCMC::fMCMCBoundaryMax [protected]

Definition at line 594 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCBoundaryMin [protected]

Parameter boundaries

Definition at line 593 of file BCEngineMCMC.h.

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

Definition at line 621 of file BCEngineMCMC.h.

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

Definition at line 616 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCCycle [protected]

The cycle of the pre-run

Definition at line 727 of file BCEngineMCMC.h.

The maximum allowed efficiency for MCMC

Definition at line 701 of file BCEngineMCMC.h.

The minimum required efficiency for MCMC

Definition at line 697 of file BCEngineMCMC.h.

Flag for convergence

Definition at line 638 of file BCEngineMCMC.h.

Flag which controls fill histograms during main run.

Definition at line 716 of file BCEngineMCMC.h.

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 707 of file BCEngineMCMC.h.

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

Definition at line 712 of file BCEngineMCMC.h.

Defines if a prerun should be performed or not

Definition at line 682 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagRun [protected]

Defines if MCMC has been performed or not

Definition at line 686 of file BCEngineMCMC.h.

std::vector<bool> BCEngineMCMC::fMCMCFlagsFillHistograms [protected]

Parameter flags for marginalization

Definition at line 598 of file BCEngineMCMC.h.

Flag to write Markov chains to file

Definition at line 664 of file BCEngineMCMC.h.

Flag to write pre run to file

Definition at line 668 of file BCEngineMCMC.h.

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

An array of marginalized distributions

Definition at line 809 of file BCEngineMCMC.h.

std::vector<int> BCEngineMCMC::fMCMCH1NBins [protected]

Number of bins per dimension for the marginalized distributions.

Definition at line 805 of file BCEngineMCMC.h.

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

Definition at line 810 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 693 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNChains [protected]

Number of Markov chains ran in parallel

Definition at line 602 of file BCEngineMCMC.h.

std::vector<int> BCEngineMCMC::fMCMCNIterations [protected]

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

Definition at line 611 of file BCEngineMCMC.h.

Number of iterations needed for all chains to convergence simulaneously

Definition at line 634 of file BCEngineMCMC.h.

Maximum number of iterations for a Markov chain prerun

Definition at line 642 of file BCEngineMCMC.h.

Minimum number of iterations for the pre-run

Definition at line 650 of file BCEngineMCMC.h.

Number of iterations for a Markov chain run

Definition at line 646 of file BCEngineMCMC.h.

Number of iterations for updating scale factors

Definition at line 625 of file BCEngineMCMC.h.

Maximum number of iterations for updating scale factors

Definition at line 629 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNLag [protected]

The lag for the Markov Chain

Definition at line 606 of file BCEngineMCMC.h.

Number of parameters

Definition at line 589 of file BCEngineMCMC.h.

std::vector<int> BCEngineMCMC::fMCMCNTrialsFalse [protected]

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

Definition at line 660 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 655 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCPhase [protected]

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

Definition at line 722 of file BCEngineMCMC.h.

Pointer to a member function

Definition at line 583 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 759 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 764 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 769 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 774 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValue [protected]

The R-value at which the chains did converge

Definition at line 794 of file BCEngineMCMC.h.

The R-value criterion for convergence of log-likelihood

Definition at line 786 of file BCEngineMCMC.h.

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

The R-values for each parameter

Definition at line 797 of file BCEngineMCMC.h.

The R-value criterion for convergence of parameters

Definition at line 790 of file BCEngineMCMC.h.

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 782 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 815 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 673 of file BCEngineMCMC.h.

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

Start values of the scale factors for the trial functions.

Definition at line 678 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 734 of file BCEngineMCMC.h.

std::vector<double> BCEngineMCMC::fMCMCxLocal [protected]

A temporary vector for a single Markov chain

Definition at line 754 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 740 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 745 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 750 of file BCEngineMCMC.h.

TRandom3* BCEngineMCMC::fRandom [protected]

Random number generator

Definition at line 801 of file BCEngineMCMC.h.


The documentation for this class was generated from the following files: