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

#include <BCEngineMCMC.h>

Inheritance diagram for BCEngineMCMC:

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

BCEngineMCMC &	operator= (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 ()
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	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 (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	MCMCSetRValueCriterion (double r)
void	MCMCSetRValueParametersCriterion (double r)
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 (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
std::vector< double >	fMCMCBoundaryMax
std::vector< double >	fMCMCBoundaryMin
int	fMCMCCurrentChain
int	fMCMCCurrentIteration
int	fMCMCCycle
double	fMCMCEfficiencyMax
double	fMCMCEfficiencyMin
bool	fMCMCFlagConvergenceGlobal
bool	fMCMCFlagFillHistograms
int	fMCMCFlagInitialPosition
bool	fMCMCFlagOrderParameters
bool	fMCMCFlagPreRun
bool	fMCMCFlagRun
std::vector< bool >	fMCMCFlagsFillHistograms
bool	fMCMCFlagWriteChainToFile
bool	fMCMCFlagWritePreRunToFile
std::vector< TH1D * >	fMCMCH1Marginalized
std::vector< int >	fMCMCH1NBins
std::vector< TH2D * >	fMCMCH2Marginalized
std::vector< double >	fMCMCInitialPosition
int	fMCMCNChains
std::vector< int >	fMCMCNIterations
int	fMCMCNIterationsConvergenceGlobal
int	fMCMCNIterationsMax
int	fMCMCNIterationsPreRunMin
int	fMCMCNIterationsRun
int	fMCMCNIterationsUpdate
int	fMCMCNIterationsUpdateMax
int	fMCMCNLag
int	fMCMCNParameters
std::vector< int >	fMCMCNTrialsFalse
std::vector< int >	fMCMCNTrialsTrue
int	fMCMCPhase
std::vector< double >	fMCMCprob
std::vector< double >	fMCMCprobMax
std::vector< double >	fMCMCprobMean
std::vector< double >	fMCMCprobVar
double	fMCMCRValue
double	fMCMCRValueCriterion
std::vector< double >	fMCMCRValueParameters
double	fMCMCRValueParametersCriterion
std::vector< TTree * >	fMCMCTrees
std::vector< double >	fMCMCTrialFunctionScaleFactor
std::vector< double >	fMCMCTrialFunctionScaleFactorStart
std::vector< double >	fMCMCx
std::vector< double >	fMCMCxLocal
std::vector< double >	fMCMCxMax
std::vector< double >	fMCMCxMean
std::vector< double >	fMCMCxVar
TRandom3 *	fRandom
Private Types
typedef bool(BCEngineMCMC::*	MCMCPointerToGetProposalPoint )(int chain, std::vector< double > xnew, std::vector< double > xold) const
Private Member Functions
void	Copy (BCEngineMCMC &enginemcmc) 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]

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

{ kLow, kMedium, kHigh, kVeryHigh };

Constructor & Destructor Documentation

BCEngineMCMC::BCEngineMCMC ( )

Default constructor.

Definition at line 25 of file BCEngineMCMC.cxx.

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

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

BCEngineMCMC::BCEngineMCMC ( int n )

Constructor.

Parameters:

n

number of chains

Definition at line 35 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 182 of file BCEngineMCMC.cxx.

{
   enginemcmc.Copy(*this);
}

BCEngineMCMC::~BCEngineMCMC ( ) [virtual]

Default destructor.

Definition at line 162 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

void BCEngineMCMC::Copy ( BCEngineMCMC & enginemcmc ) const [private]

Definition at line 347 of file BCEngineMCMC.cxx.

{}

double BCEngineMCMC::LogEval ( std::vector< double > parameters ) [virtual]

Reimplemented in BCIntegrate, and BCModel.

Definition at line 800 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
	)

Definition at line 1374 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]`

Definition at line 534 of file BCEngineMCMC.h.

{};

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]

Definition at line 201 of file BCEngineMCMC.h.

         { return fMCMCCycle; };

bool BCEngineMCMC::MCMCGetFlagConvergenceGlobal ( ) [inline]

Definition at line 117 of file BCEngineMCMC.h.

         { return fMCMCFlagConvergenceGlobal; };

int BCEngineMCMC::MCMCGetFlagInitialPosition ( ) [inline]

Definition at line 221 of file BCEngineMCMC.h.

         { return fMCMCFlagInitialPosition; };

bool BCEngineMCMC::MCMCGetFlagRun ( ) [inline]

Definition at line 247 of file BCEngineMCMC.h.

      { return fMCMCFlagRun; };

TH1D * BCEngineMCMC::MCMCGetH1Marginalized ( int i )

Definition at line 197 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
	)

Definition at line 210 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]

Definition at line 186 of file BCEngineMCMC.h.

         { return fMCMCprob; };

double BCEngineMCMC::MCMCGetLogProbx ( int ichain )

Definition at line 435 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]

Definition at line 254 of file BCEngineMCMC.h.

         { return fMCMCTrees.at(i); };

std::vector<double> BCEngineMCMC::MCMCGetMaximumLogProb ( ) [inline]

Definition at line 216 of file BCEngineMCMC.h.

         { return fMCMCprobMax; };

std::vector< double > BCEngineMCMC::MCMCGetMaximumPoint ( int i )

Definition at line 235 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]

Definition at line 206 of file BCEngineMCMC.h.

         { return fMCMCxMax; };

int BCEngineMCMC::MCMCGetNChains ( ) [inline]

Definition at line 86 of file BCEngineMCMC.h.

         { return fMCMCNChains; };

bool BCEngineMCMC::MCMCGetNewPointMetropolis ( int chain = 0 )

Definition at line 555 of file BCEngineMCMC.cxx.

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

   fMCMCCurrentChain = chain;

   // increase counter
   fMCMCNIterations[chain]++;

   // get proposal point
   if (!this->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 = this->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 485 of file BCEngineMCMC.cxx.

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

   fMCMCCurrentChain = chain;

   // increase counter
   fMCMCNIterations[chain]++;

   // get proposal point
   if (!this->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 = this->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]

Definition at line 96 of file BCEngineMCMC.h.

         { return fMCMCNIterations; };

int BCEngineMCMC::MCMCGetNIterationsConvergenceGlobal ( ) [inline]

Definition at line 112 of file BCEngineMCMC.h.

         { return fMCMCNIterationsConvergenceGlobal; };

int BCEngineMCMC::MCMCGetNIterationsMax ( ) [inline]

Definition at line 122 of file BCEngineMCMC.h.

         { return fMCMCNIterationsMax; };

int BCEngineMCMC::MCMCGetNIterationsRun ( ) [inline]

Definition at line 127 of file BCEngineMCMC.h.

         { return fMCMCNIterationsRun; };

int BCEngineMCMC::MCMCGetNLag ( ) [inline]

Definition at line 91 of file BCEngineMCMC.h.

         { return fMCMCNLag; };

int BCEngineMCMC::MCMCGetNParameters ( ) [inline]

Definition at line 81 of file BCEngineMCMC.h.

         { return fMCMCNParameters; };

std::vector<int> BCEngineMCMC::MCMCGetNTrialsFalse ( ) [inline]

Definition at line 137 of file BCEngineMCMC.h.

         { return fMCMCNTrialsFalse; };

std::vector<int> BCEngineMCMC::MCMCGetNTrialsTrue ( ) [inline]

Definition at line 132 of file BCEngineMCMC.h.

         { return fMCMCNTrialsTrue; };

int BCEngineMCMC::MCMCGetPhase ( ) [inline]

Definition at line 196 of file BCEngineMCMC.h.

         { return fMCMCPhase; };

std::vector<double> BCEngineMCMC::MCMCGetprobMean ( ) [inline]

Definition at line 143 of file BCEngineMCMC.h.

         { return fMCMCprobMean; };

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

Definition at line 446 of file BCEngineMCMC.cxx.

{
   // get unscaled random point. this point might not be in the correct volume.
   this->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
	)

Definition at line 464 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]

Definition at line 236 of file BCEngineMCMC.h.

         { return fMCMCRValue; };

double BCEngineMCMC::MCMCGetRValueCriterion ( ) [inline]

Definition at line 226 of file BCEngineMCMC.h.

         { return fMCMCRValueCriterion; };

double BCEngineMCMC::MCMCGetRValueParameters ( int i ) [inline]

Definition at line 242 of file BCEngineMCMC.h.

         { return fMCMCRValueParameters.at(i); };

double BCEngineMCMC::MCMCGetRValueParametersCriterion ( ) [inline]

Definition at line 231 of file BCEngineMCMC.h.

         { return fMCMCRValueParametersCriterion; };

TRandom3* BCEngineMCMC::MCMCGetTRandom3 ( ) [inline]

Definition at line 272 of file BCEngineMCMC.h.

         { return fRandom; };

std::vector<double> BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( ) [inline]

Definition at line 154 of file BCEngineMCMC.h.

         { return fMCMCTrialFunctionScaleFactor; };

std::vector< double > BCEngineMCMC::MCMCGetTrialFunctionScaleFactor ( int ichain )

Definition at line 371 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
	)

Definition at line 388 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]

Definition at line 149 of file BCEngineMCMC.h.

         { return fMCMCprobVar; };

std::vector<double> BCEngineMCMC::MCMCGetx ( ) [inline]

Definition at line 170 of file BCEngineMCMC.h.

         { return fMCMCx; };

std::vector< double > BCEngineMCMC::MCMCGetx ( int ichain )

Definition at line 403 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
	)

Definition at line 420 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 ( )

Definition at line 629 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 ( )

Definition at line 681 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 ( )

Definition at line 709 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)
      {
         double sum = 0;
         double sum2 = 0;
         double sumv = 0;

         // loop over chains
         for (int ichains = 0; ichains < fMCMCNChains; ++ichains) {
            // get parameter index
            int index = ichains * fMCMCNParameters + iparameters;

            // add to sums
            sum  += fMCMCxMean[index];
            sum2 += fMCMCxMean[index] * fMCMCxMean[index];
            sumv += fMCMCxVar[index];
         }

         // calculate r-value for each parameter
         double mean = sum / double(fMCMCNChains);
         double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints);
         double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);

         double r = 100.0;

         // check denominator and convergence
         if (W > 0) {
            r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W);
            fMCMCRValueParameters[iparameters] = r;
            
            // 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
      }
      // convergence criterion applied on the log-likelihood
      double sum = 0;
      double sum2 = 0;
      double sumv = 0;

      // loop over chains
      for (int i = 0; i < fMCMCNChains; ++i)
      {
         sum  += fMCMCprobMean[i];
         sum2 += fMCMCprobMean[i] * fMCMCprobMean[i]; ;
         sumv += fMCMCprobVar[i];
      }

      // calculate r-value for log-likelihood
      double mean = sum / double(fMCMCNChains);
      double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints);
      double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);
      double r = 100.0;

      if (W > 0)
      {
         r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W);
         fMCMCRValue = r;

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

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

void BCEngineMCMC::MCMCInChainUpdateStatistics ( )

Definition at line 648 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 ( )

Definition at line 792 of file BCEngineMCMC.cxx.

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

int BCEngineMCMC::MCMCInitialize ( )

Definition at line 1448 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 ( )

Definition at line 1391 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] = this->LogEval(x0);
   }

   x0.clear();
}

void BCEngineMCMC::MCMCInitializeMarkovChainTrees ( )

Definition at line 325 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]

Reimplemented in BCIntegrate.

Definition at line 523 of file BCEngineMCMC.h.

{};

int BCEngineMCMC::MCMCMetropolis ( )

Definition at line 1199 of file BCEngineMCMC.cxx.

{
   // check if prerun has been performed
   if (!fMCMCFlagPreRun)
      this->MCMCMetropolisPreRun();

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

   // reset run statistics
   this->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)
               this->MCMCGetNewPointMetropolis(ichains, iparameters);

            // reset current chain
            fMCMCCurrentChain = -1;

            // update search for maximum
            this->MCMCInChainCheckMaximum();

         } // end loop over all parameters

         // check if the current iteration is consistent with the lag
         if ( fMCMCCurrentIteration % fMCMCNLag == 0)
         {
            // fill histograms
            this->MCMCInChainFillHistograms();
               
            // write chain to file
            if (fMCMCFlagWriteChainToFile)
               this->MCMCInChainWriteChains();

            // do anything interface
            this->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
            this->MCMCGetNewPointMetropolis(ichains);

         // reset current chain
         fMCMCCurrentChain = -1;
               
         // update search for maximum
         this->MCMCInChainCheckMaximum();

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

            // write chain to file
            if (fMCMCFlagWriteChainToFile)
               this->MCMCInChainWriteChains();

            // do anything interface
            this->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:");
   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 coutner
   fMCMCCurrentIteration = -1;

   // reset current chain
   fMCMCCurrentChain = -1;

   // set flags
// fMCMCFlagPreRun = false;
   fMCMCFlagRun = true;

   return 1;
}

int BCEngineMCMC::MCMCMetropolisPreRun ( )

Definition at line 808 of file BCEngineMCMC.cxx.

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

   // initialize Markov chain
   this->MCMCInitialize();
   this->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
   this->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)
               this->MCMCGetNewPointMetropolis(ichains, iparameters);

            // search for global maximum
            this->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
            this->MCMCGetNewPointMetropolis(ichains);

         // search for global maximum
         this->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 
         this->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)
         this->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;

   // set pre-run flag
   fMCMCFlagPreRun = true;

   // reset current iteration
   fMCMCCurrentIteration = -1;

   // reset current chain
   fMCMCCurrentChain = -1;

   // no error
   return 1;
}

int BCEngineMCMC::MCMCResetResults ( )

Definition at line 1408 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 = false;
   fMCMCFlagRun = false;
   fMCMCFlagConvergenceGlobal = false;

   // no errors
   return 1;
}

void BCEngineMCMC::MCMCResetRunStatistics ( )

Definition at line 1344 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 )

Definition at line 291 of file BCEngineMCMC.cxx.

{
   fMCMCFlagFillHistograms = flag;

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

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

Definition at line 300 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]

Definition at line 356 of file BCEngineMCMC.h.

         { fMCMCFlagInitialPosition = flag; };

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

Definition at line 362 of file BCEngineMCMC.h.

         { fMCMCFlagOrderParameters = flag; };

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

Definition at line 277 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));

   this->MCMCSetInitialPositions(y0s);
}

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

Definition at line 261 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
   this->MCMCSetFlagInitialPosition(2);
}

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

Definition at line 314 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]

Definition at line 331 of file BCEngineMCMC.h.

         { fMCMCEfficiencyMax = efficiency; };

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

Definition at line 326 of file BCEngineMCMC.h.

         { fMCMCEfficiencyMin = efficiency; };

void BCEngineMCMC::MCMCSetNChains ( int n )

Definition at line 252 of file BCEngineMCMC.cxx.

{
   fMCMCNChains = n;

   // re-initialize
   this->MCMCInitialize();
}

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

Definition at line 296 of file BCEngineMCMC.h.

         { fMCMCNIterationsMax = n; };

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

Definition at line 306 of file BCEngineMCMC.h.

         { fMCMCNIterationsPreRunMin = n; };

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

Definition at line 301 of file BCEngineMCMC.h.

         { fMCMCNIterationsRun = n; };

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

Definition at line 313 of file BCEngineMCMC.h.

         { fMCMCNIterationsUpdate = n; };

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

Definition at line 321 of file BCEngineMCMC.h.

         { fMCMCNIterationsUpdateMax = n; };

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

Definition at line 291 of file BCEngineMCMC.h.

         { fMCMCNLag = n; };

void BCEngineMCMC::MCMCSetPrecision ( BCEngineMCMC::Precision precision )

Set the precision for the MCMC run.

Definition at line 104 of file BCEngineMCMC.cxx.

{
   switch(precision) {
   case BCEngineMCMC::kLow:
      fMCMCNChains              = 1;
      fMCMCNLag                 = 1;
      fMCMCNIterationsMax       = 10000;
      fMCMCNIterationsRun       = 10000;
      fMCMCNIterationsPreRunMin = 100;
      fMCMCNIterationsUpdate    = 1000;
      fMCMCNIterationsUpdateMax = 10000;
      fMCMCRValueCriterion      = 0.1;
      fMCMCRValueParametersCriterion = 0.1;
      fMCMCRValue               = 100;
      break;
   case  BCEngineMCMC::kMedium:
      fMCMCNChains              = 5;
      fMCMCNLag                 = 1;
      fMCMCNIterationsMax       = 100000;
      fMCMCNIterationsRun       = 100000;
      fMCMCNIterationsPreRunMin = 100;
      fMCMCNIterationsUpdate    = 1000;
      fMCMCNIterationsUpdateMax = 10000;
      fMCMCRValueCriterion      = 0.1;
      fMCMCRValueParametersCriterion = 0.1;
      fMCMCRValue               = 100;
      break;
   case  BCEngineMCMC::kHigh:
      fMCMCNChains              = 10;
      fMCMCNLag                 = 10;
      fMCMCNIterationsMax       = 1000000;
      fMCMCNIterationsRun       = 1000000;
      fMCMCNIterationsPreRunMin = 100;
      fMCMCNIterationsUpdate    = 1000;
      fMCMCNIterationsUpdateMax = 10000;
      fMCMCRValueCriterion      = 0.1;
      fMCMCRValueParametersCriterion = 0.1;
      fMCMCRValue               = 100;
      break;
   case  BCEngineMCMC::kVeryHigh:
      fMCMCNChains              = 10;
      fMCMCNLag                 = 10;
      fMCMCNIterationsMax       = 10000000;
      fMCMCNIterationsRun       = 10000000;
      fMCMCNIterationsPreRunMin = 100;
      fMCMCNIterationsUpdate    = 1000;
      fMCMCNIterationsUpdateMax = 10000;
      fMCMCRValueCriterion      = 0.1;
      fMCMCRValueParametersCriterion = 0.1;
      fMCMCRValue               = 100;    
      break;
   }

   // re-initialize
   MCMCInitialize();
}

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

Definition at line 373 of file BCEngineMCMC.h.

         { fMCMCRValueCriterion = r; };

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

Definition at line 378 of file BCEngineMCMC.h.

         { fMCMCRValueParametersCriterion = r; };

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

Definition at line 282 of file BCEngineMCMC.h.

         { fMCMCTrialFunctionScaleFactorStart = scale; };

void BCEngineMCMC::MCMCSetValuesDefault ( )

Definition at line 45 of file BCEngineMCMC.cxx.

{
   fMCMCNParameters          = 0;
   fMCMCFlagWriteChainToFile = false;
   fMCMCFlagWritePreRunToFile = false;
   fMCMCFlagPreRun           = false;
   fMCMCFlagRun              = false;
   fMCMCFlagFillHistograms   = true;
   fMCMCEfficiencyMin        = 0.15;
   fMCMCEfficiencyMax        = 0.50;
   fMCMCFlagInitialPosition  = 1;
   fMCMCNLag                 = 1;
   fMCMCCurrentIteration     = -1;
   fMCMCCurrentChain         = -1;

   this->MCMCSetValuesDetail();
}

void BCEngineMCMC::MCMCSetValuesDetail ( )

Definition at line 84 of file BCEngineMCMC.cxx.

{
   fMCMCNChains              = 5;
   fMCMCNIterationsMax       = 1000000;
   fMCMCNIterationsRun       = 100000;
   fMCMCNIterationsPreRunMin = 500;
   fMCMCFlagInitialPosition  = 1;
   fMCMCRValueCriterion      = 0.1;
   fMCMCRValueParametersCriterion = 0.1;
   fMCMCNIterationsConvergenceGlobal = -1;
   fMCMCFlagConvergenceGlobal = false;
   fMCMCRValue               = 100;
   fMCMCNIterationsUpdate    = 1000;
   fMCMCNIterationsUpdateMax = 10000;
   fMCMCFlagOrderParameters  = true;
   fMCMCCurrentIteration     = -1;
   fMCMCCurrentChain         = -1;
}

void BCEngineMCMC::MCMCSetValuesQuick ( )

Definition at line 64 of file BCEngineMCMC.cxx.

{
   fMCMCNChains              = 1;
   fMCMCNIterationsMax       = 1000;
   fMCMCNIterationsRun       = 10000;
   fMCMCNIterationsPreRunMin = 0;
   fMCMCFlagInitialPosition  = 1;
   fMCMCRValueCriterion      = 0.1;
   fMCMCRValueParametersCriterion = 0.1;
   fMCMCNIterationsConvergenceGlobal = -1;
   fMCMCFlagConvergenceGlobal = false;
   fMCMCRValue               = 100;
   fMCMCNIterationsUpdate    = 1000;
   fMCMCNIterationsUpdateMax = 10000;
   fMCMCFlagOrderParameters  = true;
   fMCMCCurrentIteration     = -1;
   fMCMCCurrentChain         = -1;
}

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

Definition at line 336 of file BCEngineMCMC.h.

         { fMCMCFlagWriteChainToFile = flag; };

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

Definition at line 341 of file BCEngineMCMC.h.

         { fMCMCFlagWritePreRunToFile = flag; };

void BCEngineMCMC::MCMCTrialFunction	(	int	ichain,
		std::vector< double > &	x
	)			`[virtual]`

Definition at line 351 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]`

Definition at line 359 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 188 of file BCEngineMCMC.cxx.

{
   if (this != &enginemcmc)
      enginemcmc.Copy(* this);

   return * this;
}

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

Definition at line 1553 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
	)

Definition at line 1570 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;
}