BCEngineMCMC Class Reference

#include <BCEngineMCMC.h>

Inheritance diagram for BCEngineMCMC:

[legend]

List of all members.

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


Public Member Functions
Constructors and destructors

	BCEngineMCMC (const BCEngineMCMC &enginemcmc)
	BCEngineMCMC (int n)
	BCEngineMCMC ()
virtual	~BCEngineMCMC ()
Miscellaneous methods

virtual double	LogEval (std::vector< double > parameters)
int	MCMCAddParameter (double min, double max)
bool	MCMCGetNewPointMetropolis (int chain, int parameter, bool pca=false)
bool	MCMCGetNewPointMetropolis (int chain=0, bool pca=false)
bool	MCMCGetNewPointMetropolisHastings (int chain=0)
void	MCMCGetNewPointPCA ()
bool	MCMCGetProposalPointMetropolis (int chain, int parameter, std::vector< double > &x, bool pca)
bool	MCMCGetProposalPointMetropolis (int chain, std::vector< double > &x, bool pca)
bool	MCMCGetProposalPointMetropolisHastings (int chain, std::vector< double > &xnew, std::vector< double > &xold)
int	MCMCInitialize ()
void	MCMCInitializeMarkovChains ()
virtual void	MCMCIterationInterface ()
int	MCMCMetropolis ()
int	MCMCMetropolisHastings ()
int	MCMCMetropolisPreRun ()
void	MCMCResetRunStatistics ()
void	MCMCSetValuesDefault ()
void	MCMCSetValuesDetail ()
void	MCMCSetValuesQuick ()
void	MCMCTrialFunction (std::vector< double > &x)
void	MCMCTrialFunctionAuto (std::vector< double > &x)
double	MCMCTrialFunctionIndependent (std::vector< double > &xnew, std::vector< double > &xold, bool newpoint)
void	MCMCTrialFunctionSingle (int ichain, int iparameter, std::vector< double > &x)
void	MCMCUpdateStatistics ()
void	MCMCUpdateStatisticsCheckMaximum ()
void	MCMCUpdateStatisticsFillHistograms ()
void	MCMCUpdateStatisticsModeConvergence ()
void	MCMCUpdateStatisticsTestConvergenceAllChains ()
void	MCMCUpdateStatisticsWriteChains ()
int	SetMarginalized (int index1, int index2, TH2D *h)
int	SetMarginalized (int index, TH1D *h)
Getters

bool	MCMCGetFlagConvergenceGlobal ()
int	MCMCGetFlagInitialPosition ()
bool	MCMCGetFlagIterationsAuto ()
bool	MCMCGetFlagPCA ()
TH1D *	MCMCGetH1Marginalized (int i)
TH2D *	MCMCGetH2Marginalized (int i, int j)
double	MCMCGetLogProbx (int ichain)
std::vector< double >	MCMCGetLogProbx ()
TTree *	MCMCGetMarkovChainTree (int i)
std::vector< double >	MCMCGetMaximumLogProb ()
std::vector< double >	MCMCGetMaximumPoint (int i)
std::vector< double >	MCMCGetMaximumPoints ()
std::vector< double >	MCMCGetMean ()
int	MCMCGetNChains ()
std::vector< int >	MCMCGetNIterations ()
int	MCMCGetNIterationsBurnIn ()
int	MCMCGetNIterationsConvergenceGlobal ()
std::vector< int >	MCMCGetNIterationsConvergenceLocal ()
int	MCMCGetNIterationsMax ()
int	MCMCGetNIterationsPCA ()
int	MCMCGetNParameters ()
std::vector< int >	MCMCGetNTrialsFalse ()
std::vector< int >	MCMCGetNTrialsTrue ()
std::vector< double > *	MCMCGetP2LogProbx ()
std::vector< int > *	MCMCGetP2NIterations ()
std::vector< double > *	MCMCGetP2x ()
double	MCMCGetRValue ()
double	MCMCGetRValueCriterion ()
double	MCMCGetRValueParameters (int i)
double	MCMCGetRValueParametersCriterion ()
double	MCMCGetTrialFunctionScale ()
double	MCMCGetTrialFunctionScaleFactor (int ichain, int ipar)
std::vector< double >	MCMCGetTrialFunctionScaleFactor (int ichain)
std::vector< double >	MCMCGetTrialFunctionScaleFactor ()
std::vector< double >	MCMCGetVariance ()
double	MCMCGetx (int ichain, int ipar)
std::vector< double >	MCMCGetx (int ichain)
std::vector< double >	MCMCGetx ()
Setters

void	MCMCSetFlagInitialPosition (int flag)
void	MCMCSetFlagOrderParameters (bool flag)
void	MCMCSetFlagPCA (bool flag)
void	MCMCSetFlagPCATruncate (bool flag)
void	MCMCSetInitialPosition (int chain, std::vector< double > x0)
void	MCMCSetInitialPosition (std::vector< double > x0)
void	MCMCSetInitialPositions (std::vector< double > x0s)
void	MCMCSetIterationsAuto (bool flag)
void	MCMCSetMarkovChainTrees (std::vector< TTree * > trees)
void	MCMCSetMaximumEfficiency (double efficiency)
void	MCMCSetMinimumEfficiency (double efficiency)
void	MCMCSetNChains (int n)
void	MCMCSetNIterationsBurnIn (int n)
void	MCMCSetNIterationsMax (int n)
void	MCMCSetNIterationsPCA (int n)
void	MCMCSetNIterationsRun (int n)
void	MCMCSetNIterationsUpdate (int n)
void	MCMCSetPCAMinimumRatio (double ratio)
void	MCMCSetRValueCriterion (double r)
void	MCMCSetRValueParametersCriterion (double r)
void	MCMCSetTrialFunctionScale (double scale)
void	MCMCSetTrialFunctionScaleFactor (std::vector< double > scale)
void	MCMCSetWriteChainToFile (bool flag)
Assignment operators

BCEngineMCMC &	operator= (const BCEngineMCMC &engineMCMC)
Protected Attributes
std::vector< double >	fMCMCBoundaryMax
std::vector< double >	fMCMCBoundaryMin
double	fMCMCEfficiencyMax
double	fMCMCEfficiencyMin
bool	fMCMCFlagConvergenceGlobal
int	fMCMCFlagInitialPosition
bool	fMCMCFlagIterationsAuto
bool	fMCMCFlagOrderParameters
bool	fMCMCFlagPCA
bool	fMCMCFlagPCATruncate
bool	fMCMCFlagPreRun
bool	fMCMCFlagWriteChainToFile
std::vector< TH1D * >	fMCMCH1Marginalized
std::vector< int >	fMCMCH1NBins
std::vector< TH2D * >	fMCMCH2Marginalized
std::vector< double >	fMCMCInitialPosition
std::vector< double >	fMCMCLogProbx
std::vector< double >	fMCMCMaximumLogProb
std::vector< double >	fMCMCMaximumPoints
std::vector< double >	fMCMCMean
std::vector< double >	fMCMCMeanx
int	fMCMCNChains
int	fMCMCNDimensionsPCA
std::vector< int >	fMCMCNIterations
int	fMCMCNIterationsBurnIn
int	fMCMCNIterationsConvergenceGlobal
std::vector< int >	fMCMCNIterationsConvergenceLocal
int	fMCMCNIterationsMax
int	fMCMCNIterationsPCA
int	fMCMCNIterationsPreRunMin
int	fMCMCNIterationsRun
int	fMCMCNIterationsUpdate
int	fMCMCNParameters
std::vector< int >	fMCMCNTrialsFalse
std::vector< int >	fMCMCNTrialsTrue
std::vector< double >	fMCMCNumericalPrecisionModeValues
TPrincipal *	fMCMCPCA
std::vector< double >	fMCMCPCAMean
double	fMCMCPCAMinimumRatio
std::vector< double >	fMCMCPCAVariance
TRandom3 *	fMCMCRandom
double	fMCMCRValue
double	fMCMCRValueCriterion
std::vector< double >	fMCMCRValueParameters
double	fMCMCRValueParametersCriterion
std::vector< double >	fMCMCSum
std::vector< double >	fMCMCSum2
std::vector< TTree * >	fMCMCTrees
double	fMCMCTrialFunctionScale
std::vector< double >	fMCMCTrialFunctionScaleFactor
std::vector< double >	fMCMCTrialFunctionScaleFactorStart
std::vector< double >	fMCMCVariance
std::vector< double >	fMCMCVariancex
std::vector< double >	fMCMCx
std::vector< double >	fMCMCxLocal
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

Member Typedef Documentation

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

Constructor & Destructor Documentation

BCEngineMCMC::BCEngineMCMC ( )

Default constructor.

Definition at line 31 of file BCEngineMCMC.cxx.

00032 {
00033    // default settings
00034    fMCMCNParameters          = 0;
00035    fMCMCFlagWriteChainToFile = false;
00036    fMCMCFlagPreRun           = false;
00037    fMCMCEfficiencyMin        = 0.15;
00038    fMCMCEfficiencyMax        = 0.50;
00039 
00040    this -> MCMCSetValuesDefault();
00041 
00042 // fMCMCRelativePrecisionMode = 1e-3;
00043 
00044    // set pointer to control histograms to NULL
00045 // for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
00046 //    fMCMCH1Marginalized[i] = 0;
00047 
00048 // for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
00049 //    fMCMCH2Marginalized[i] = 0;
00050 
00051 // fMCMCH1RValue = 0;
00052 // fMCMCH1Efficiency = 0;
00053 
00054    // initialize random number generator
00055    fMCMCRandom = new TRandom3(0);
00056 
00057    // initialize
00058 // this -> MCMCInitialize();
00059 }

BCEngineMCMC::BCEngineMCMC ( int n )

Constructor.

Parameters:

n

number of chains

Definition at line 63 of file BCEngineMCMC.cxx.

00064 {
00065    // set number of chains to n
00066    fMCMCNChains = n;
00067 
00068    // call default constructor
00069    BCEngineMCMC();
00070 }

BCEngineMCMC::BCEngineMCMC ( const BCEngineMCMC & enginemcmc )

Default copy constructor.

Definition at line 101 of file BCEngineMCMC.cxx.

00102 {
00103    enginemcmc.Copy(*this);
00104 }

BCEngineMCMC::~BCEngineMCMC ( ) [virtual]

Default destructor.

Definition at line 74 of file BCEngineMCMC.cxx.

00075 {
00076    // delete random number generator
00077    if (fMCMCRandom) delete fMCMCRandom;
00078 
00079    // delete PCA object
00080    if (fMCMCPCA) delete fMCMCPCA;
00081 
00082    // delete constrol histograms and plots
00083    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
00084       if (fMCMCH1Marginalized[i])
00085          delete fMCMCH1Marginalized[i];
00086 
00087    fMCMCH1Marginalized.clear();
00088 
00089    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
00090       if (fMCMCH2Marginalized[i])
00091          delete fMCMCH2Marginalized[i];
00092 
00093    fMCMCH2Marginalized.clear();
00094 
00095 // if (fMCMCH1RValue) delete fMCMCH1RValue;
00096 // if (fMCMCH1Efficiency) delete fMCMCH1Efficiency;
00097 }

Member Function Documentation

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

Definition at line 245 of file BCEngineMCMC.cxx.

00246 {}

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

Reimplemented in BCIntegrate, and BCModel.

Definition at line 926 of file BCEngineMCMC.cxx.

00927 {
00928    // test function for now
00929    // this will be overloaded by the user
00930    return 0.0;
00931 }

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

Definition at line 1672 of file BCEngineMCMC.cxx.

01673 {
01674    // add the boundaries to the corresponding vectors
01675    fMCMCBoundaryMin.push_back(min);
01676    fMCMCBoundaryMax.push_back(max);
01677 
01678    // increase the number of parameters by one
01679    fMCMCNParameters++;
01680 
01681    // return the number of parameters
01682    return fMCMCNParameters;
01683 }

bool BCEngineMCMC::MCMCGetFlagConvergenceGlobal ( ) [inline]

Definition at line 108 of file BCEngineMCMC.h.

00109          { return fMCMCFlagConvergenceGlobal; };

int BCEngineMCMC::MCMCGetFlagInitialPosition ( ) [inline]

Definition at line 238 of file BCEngineMCMC.h.

00239          { return fMCMCFlagInitialPosition; };

bool BCEngineMCMC::MCMCGetFlagIterationsAuto ( ) [inline]

Definition at line 153 of file BCEngineMCMC.h.

00154          { return fMCMCFlagIterationsAuto; };

bool BCEngineMCMC::MCMCGetFlagPCA ( ) [inline]

Definition at line 269 of file BCEngineMCMC.h.

00270          { return fMCMCFlagPCA; };

TH1D* BCEngineMCMC::MCMCGetH1Marginalized ( int i ) [inline]

Definition at line 283 of file BCEngineMCMC.h.

00284          { return fMCMCH1Marginalized[i]; };

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

Definition at line 118 of file BCEngineMCMC.cxx.

00119 {
00120    int counter = 0;
00121    int index = 0;
00122 
00123    // search for correct combination
00124    for(int i = 0; i < fMCMCNParameters; i++)
00125       for (int j = 0; j < i; ++j)
00126       {
00127          if(j == index1 && i == index2)
00128             index = counter;
00129          counter++;
00130       }
00131 
00132    return fMCMCH2Marginalized[index];
00133 }

double BCEngineMCMC::MCMCGetLogProbx ( int ichain )

Definition at line 362 of file BCEngineMCMC.cxx.

00363 {
00364    // check if ichain is in range
00365    if (ichain < 0 || ichain >= fMCMCNChains)
00366       return -1;
00367 
00368    // return log of the probability at the current point in the ith chain
00369    return fMCMCLogProbx.at(ichain);
00370 }

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

Definition at line 200 of file BCEngineMCMC.h.

00201          { return fMCMCLogProbx; };

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

Definition at line 276 of file BCEngineMCMC.h.

00277          { return fMCMCTrees.at(i); };

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

Definition at line 225 of file BCEngineMCMC.h.

00226          { return fMCMCMaximumLogProb; };

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

Definition at line 137 of file BCEngineMCMC.cxx.

00138 {
00139    // create a new vector with the lenght of fMCMCNParameters
00140    std::vector <double> x;
00141 
00142    // check if i is in range
00143    if (i < 0 || i >= fMCMCNChains)
00144       return x;
00145 
00146    // copy the point in the ith chain into the temporary vector
00147    for (int j = 0; j < fMCMCNParameters; ++j)
00148    x.push_back(fMCMCMaximumPoints.at(i * fMCMCNParameters + j));
00149 
00150    return x;
00151 }

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

Definition at line 215 of file BCEngineMCMC.h.

00216          { return fMCMCMaximumPoints; };

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

Definition at line 141 of file BCEngineMCMC.h.

00142          { return fMCMCMean; };

int BCEngineMCMC::MCMCGetNChains ( ) [inline]

Definition at line 81 of file BCEngineMCMC.h.

00082          { return fMCMCNChains; };

bool BCEngineMCMC::MCMCGetNewPointMetropolis	(	int	chain,
		int	parameter,
		bool	pca = `false`
	)

Definition at line 518 of file BCEngineMCMC.cxx.

00519 {
00520    // calculate index
00521    int index = chain * fMCMCNParameters;
00522 
00523    // get proposal point
00524    int counter = 0;
00525 
00526    while (!this -> MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal, pca) && counter < 1000)
00527       counter++;
00528 
00529    // calculate probabilities of the old and new points
00530    double p0 = fMCMCLogProbx[chain];
00531    double p1 = this -> LogEval(fMCMCxLocal);
00532 
00533    // flag for accept
00534    bool accept = false;
00535 
00536    // if the new point is more probable, keep it ...
00537    if (p1 >= p0)
00538       accept = true;
00539    // ... or else throw dice.
00540    else
00541    {
00542       double r = log(fMCMCRandom -> Rndm());
00543 
00544       if(r < p1 - p0)
00545          accept = true;
00546    }
00547 
00548    // fill the new point
00549    if(accept)
00550    {
00551       // increase counter
00552       fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++;
00553 
00554       // copy the point
00555       for(int i = 0; i < fMCMCNParameters; ++i)
00556       {
00557          // save the point
00558          fMCMCx[index + i] = fMCMCxLocal[i];
00559 
00560          // save the probability of the point
00561          fMCMCLogProbx[chain] = p1;
00562       }
00563    }
00564    else
00565    {
00566       // increase counter
00567       fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00568    }
00569 
00570    // increase counter
00571    fMCMCNIterations[chain]++;
00572 
00573    return accept;
00574 }

bool BCEngineMCMC::MCMCGetNewPointMetropolis	(	int	chain = `0`,
		bool	pca = `false`
	)

Definition at line 578 of file BCEngineMCMC.cxx.

00579 {
00580    // calculate index
00581    int index = chain * fMCMCNParameters;
00582 
00583    // get proposal point
00584    int counter = 0;
00585 
00586    while (!this -> MCMCGetProposalPointMetropolis(chain, fMCMCxLocal, pca) && counter < 1000)
00587       counter++;
00588 
00589    // calculate probabilities of the old and new points
00590    double p0 = fMCMCLogProbx[chain];
00591    double p1 = this -> LogEval(fMCMCxLocal);
00592 
00593    // flag for accept
00594    bool accept = false;
00595 
00596    // if the new point is more probable, keep it ...
00597    if (p1 >= p0)
00598       accept = true;
00599 
00600    // ... or else throw dice.
00601    else
00602    {
00603       double r = log(fMCMCRandom -> Rndm());
00604 
00605       if(r < p1 - p0)
00606          accept = true;
00607    }
00608 
00609    // fill the new point
00610    if(accept)
00611    {
00612       // increase counter
00613       for (int i = 0; i < fMCMCNParameters; ++i)
00614          fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00615 
00616       // copy the point
00617       for(int i = 0; i < fMCMCNParameters; ++i)
00618       {
00619          // save the point
00620          fMCMCx[index + i] = fMCMCxLocal[i];
00621 
00622          // save the probability of the point
00623          fMCMCLogProbx[chain] = p1;
00624       }
00625    }
00626    else
00627    {
00628       // increase counter
00629       for (int i = 0; i < fMCMCNParameters; ++i)
00630          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00631    }
00632 
00633    // increase counter
00634    fMCMCNIterations[chain]++;
00635 
00636    return accept;
00637 }

bool BCEngineMCMC::MCMCGetNewPointMetropolisHastings ( int chain = 0 )

Definition at line 641 of file BCEngineMCMC.cxx.

00642 {
00643    // calculate index
00644    int index = chain * fMCMCNParameters;
00645 
00646    // save old point
00647    std::vector <double> xold;
00648    for (int i = 0; i < fMCMCNParameters; ++i)
00649       xold.push_back(fMCMCxLocal.at(i));
00650 
00651    // get proposal point
00652    int counter = 0;
00653    while (!this -> MCMCGetProposalPointMetropolisHastings(chain, fMCMCxLocal, xold) && counter < 1000)
00654       counter++;
00655 
00656    // calculate probabilities of the old and new points
00657    double p0 = fMCMCLogProbx[chain] + log(this -> MCMCTrialFunctionIndependent(xold, fMCMCxLocal, false));
00658    double p1 = this -> LogEval(fMCMCxLocal) + log(this -> MCMCTrialFunctionIndependent(xold, fMCMCxLocal, false));
00659 
00660    // flag for accept
00661    bool accept = false;
00662 
00663    // if the new point is more probable, keep it ...
00664    if (p1 >= p0)
00665       accept = true;
00666    // ... or else throw dice.
00667    else
00668    {
00669       if( log(fMCMCRandom -> Rndm()) < p1 - p0)
00670          accept = true;
00671    }
00672 
00673    // fill the new point
00674    if(accept)
00675    {
00676       // increase counter
00677       for (int i = 0; i < fMCMCNParameters; ++i)
00678          fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00679 
00680       // copy the point
00681       for(int i = 0; i < fMCMCNParameters; ++i)
00682       {
00683          // save the point
00684          fMCMCx[index + i] = fMCMCxLocal[i];
00685 
00686          // save the probability of the point
00687          fMCMCLogProbx[chain] = p1;
00688       }
00689    }
00690    else
00691    {
00692       // increase counter
00693       for (int i = 0; i < fMCMCNParameters; ++i)
00694          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00695    }
00696 
00697    // increase counter
00698    fMCMCNIterations[chain]++;
00699 
00700    return accept;
00701 }

void BCEngineMCMC::MCMCGetNewPointPCA ( )

Definition at line 509 of file BCEngineMCMC.cxx.

00510 {
00511    // get random point in allowed parameter space
00512    for (int i = 0; i < fMCMCNParameters; ++i)
00513       fMCMCx[i] = fMCMCBoundaryMin.at(i) + (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i)) * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00514 }

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

Definition at line 86 of file BCEngineMCMC.h.

00087          { return fMCMCNIterations; };

int BCEngineMCMC::MCMCGetNIterationsBurnIn ( ) [inline]

Definition at line 119 of file BCEngineMCMC.h.

00120          { return fMCMCNIterationsBurnIn; };

int BCEngineMCMC::MCMCGetNIterationsConvergenceGlobal ( ) [inline]

Definition at line 103 of file BCEngineMCMC.h.

00104          { return fMCMCNIterationsConvergenceGlobal; };

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

Definition at line 97 of file BCEngineMCMC.h.

00098          { return fMCMCNIterationsConvergenceLocal; };

int BCEngineMCMC::MCMCGetNIterationsMax ( ) [inline]

Definition at line 113 of file BCEngineMCMC.h.

00114          { return fMCMCNIterationsMax; };

int BCEngineMCMC::MCMCGetNIterationsPCA ( ) [inline]

Definition at line 125 of file BCEngineMCMC.h.

00126          { return fMCMCNIterationsPCA; };

int BCEngineMCMC::MCMCGetNParameters ( ) [inline]

Definition at line 76 of file BCEngineMCMC.h.

00077          { return fMCMCNParameters; };

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

Definition at line 135 of file BCEngineMCMC.h.

00136          { return fMCMCNTrialsFalse; };

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

Definition at line 130 of file BCEngineMCMC.h.

00131          { return fMCMCNTrialsTrue; };

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

Definition at line 210 of file BCEngineMCMC.h.

00211          { return &fMCMCLogProbx; };

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

Definition at line 91 of file BCEngineMCMC.h.

00092          { return &fMCMCNIterations; };

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

Definition at line 184 of file BCEngineMCMC.h.

00185          { return &fMCMCx; };

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

Definition at line 433 of file BCEngineMCMC.cxx.

00434 {
00435    // get unscaled random point in the dimension of the chosen
00436    // parameter. this point might not be in the correct volume.
00437    this -> MCMCTrialFunctionSingle(chain, parameter, x);
00438 
00439    // shift the point to the old point (x0) and scale it.
00440    if (pca == false)
00441    {
00442       // copy the old point into the new
00443       for (int i = 0; i < fMCMCNParameters; ++i)
00444          if (i != parameter)
00445             x[i] = fMCMCx[chain * fMCMCNParameters + i];
00446 
00447       // modify the parameter under study
00448       x[parameter] = fMCMCx[chain * fMCMCNParameters + parameter] + x[parameter] * (fMCMCBoundaryMax.at(parameter) - fMCMCBoundaryMin.at(parameter));
00449    }
00450    else
00451    {
00452       // create temporary points in x and p space
00453       double * newp = new double[fMCMCNParameters];
00454       double * newx = new double[fMCMCNParameters];
00455 
00456       for (int i = 0; i < fMCMCNParameters; i++)
00457       {
00458          newp[i] = 0.;
00459          newx[i] = 0.;
00460       }
00461 
00462       // get the old point in x space
00463       for (int i = 0; i < fMCMCNParameters; i++)
00464          newx[i] = fMCMCx[chain * fMCMCNParameters + i];
00465 
00466       // transform the old point into p space
00467       fMCMCPCA -> X2P(newx, newp);
00468 
00469       // add new trial point to old point
00470       newp[parameter] += x[parameter] * sqrt(fMCMCPCAVariance[parameter]);
00471 
00472       // transform new point back to x space
00473 //    fMCMCPCA -> P2X(newp, newx, fMCMCNParameters);
00474       fMCMCPCA -> P2X(newp, newx, fMCMCNDimensionsPCA);
00475 
00476       // copy point into vector
00477       for (int i = 0; i < fMCMCNParameters; ++i)
00478          x[i] = newx[i];
00479 
00480       delete [] newp;
00481       delete [] newx;
00482    }
00483 
00484    // check if the point is in the correct volume.
00485    for (int i = 0; i < fMCMCNParameters; ++i)
00486       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00487          return false;
00488 
00489    return true;
00490 }

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

Definition at line 374 of file BCEngineMCMC.cxx.

00375 {
00376    // get unscaled random point. this point might not be in the correct volume.
00377    this -> MCMCTrialFunction(x);
00378 
00379    // shift the point to the old point (x0) and scale it.
00380    if (pca == false)
00381    {
00382       // get a proposal point from the trial function and scale it
00383       for (int i = 0; i < fMCMCNParameters; ++i)
00384          x[i] = fMCMCx[chain * fMCMCNParameters + i] + x[i] * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));
00385    }
00386    else
00387    {
00388       // create temporary points in x and p space
00389       double * newp = new double[fMCMCNParameters];
00390       double * newx = new double[fMCMCNParameters];
00391 
00392       for (int i = 0; i < fMCMCNParameters; i++)
00393       {
00394          newp[i] = 0.;
00395          newx[i] = 0.;
00396       }
00397 
00398       // get a new trial point
00399       this -> MCMCTrialFunctionAuto(x);
00400 
00401       // get the old point in x space
00402       for (int i = 0; i < fMCMCNParameters; i++)
00403          newx[i] = fMCMCx[chain * fMCMCNParameters + i];
00404 
00405       // transform the old point into p space
00406       fMCMCPCA -> X2P(newx, newp);
00407 
00408       // add new trial point to old point
00409       for (int i = 0; i < fMCMCNParameters; i++)
00410          newp[i] += fMCMCTrialFunctionScale * x[i];
00411 
00412       // transform new point back to x space
00413 //    fMCMCPCA -> P2X(newp, newx, fMCMCNParameters);
00414       fMCMCPCA -> P2X(newp, newx, fMCMCNDimensionsPCA);
00415 
00416       // copy point into vector
00417       for (int i = 0; i < fMCMCNParameters; ++i)
00418          x[i] = newx[i];
00419 
00420       delete [] newp;
00421       delete [] newx;
00422    }
00423 
00424    // check if the point is in the correct volume.
00425    for (int i = 0; i < fMCMCNParameters; ++i)
00426       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00427          return false;
00428 
00429    return true;
00430 }

bool BCEngineMCMC::MCMCGetProposalPointMetropolisHastings	(	int	chain,
		std::vector< double > &	xnew,
		std::vector< double > &	xold
	)

Definition at line 494 of file BCEngineMCMC.cxx.

00495 {
00496    // get a scaled random point.
00497    this -> MCMCTrialFunctionIndependent(xnew, xold, true);
00498 
00499    // check if the point is in the correct volume.
00500    for (int i = 0; i < fMCMCNParameters; ++i)
00501       if ((xnew[i] < fMCMCBoundaryMin[i]) || (xnew[i] > fMCMCBoundaryMax[i]))
00502          return false;
00503 
00504    return true;
00505 }

double BCEngineMCMC::MCMCGetRValue ( ) [inline]

Definition at line 253 of file BCEngineMCMC.h.

00254          { return fMCMCRValue; };

double BCEngineMCMC::MCMCGetRValueCriterion ( ) [inline]

Definition at line 243 of file BCEngineMCMC.h.

00244          { return fMCMCRValueCriterion; };

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

Definition at line 259 of file BCEngineMCMC.h.

00260          { return fMCMCRValueParameters.at(i); };

double BCEngineMCMC::MCMCGetRValueParametersCriterion ( ) [inline]

Definition at line 248 of file BCEngineMCMC.h.

00249          { return fMCMCRValueParametersCriterion; };

double BCEngineMCMC::MCMCGetTrialFunctionScale ( ) [inline]

Definition at line 158 of file BCEngineMCMC.h.

00159          { return fMCMCTrialFunctionScale; };

double BCEngineMCMC::MCMCGetTrialFunctionScaleFactor	(	int	ichain,
		int	ipar
	)

Definition at line 312 of file BCEngineMCMC.cxx.

00313 {
00314    // check if ichain is in range
00315    if (ichain < 0 || ichain >= fMCMCNChains)
00316       return 0;
00317 
00318    // check if ipar is in range
00319    if (ipar < 0 || ipar >= fMCMCNParameters)
00320       return 0;
00321 
00322    // return component of ipar point in the ichain chain
00323    return fMCMCTrialFunctionScaleFactor.at(ichain *  fMCMCNChains + ipar);
00324 }

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

Definition at line 294 of file BCEngineMCMC.cxx.

00295 {
00296    // create a new vector with the length of fMCMCNParameters
00297    std::vector <double> x;
00298 
00299    // check if ichain is in range
00300    if (ichain < 0 || ichain >= fMCMCNChains)
00301       return x;
00302 
00303    // copy the scale factors into the temporary vector
00304    for (int j = 0; j < fMCMCNParameters; ++j)
00305       x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNParameters + j));
00306 
00307    return x;
00308 }

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

Definition at line 163 of file BCEngineMCMC.h.

00164          { return fMCMCTrialFunctionScaleFactor; };

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

Definition at line 147 of file BCEngineMCMC.h.

00148          { return fMCMCVariance; };

double BCEngineMCMC::MCMCGetx	(	int	ichain,
		int	ipar
	)

Definition at line 346 of file BCEngineMCMC.cxx.

00347 {
00348    // check if ichain is in range
00349    if (ichain < 0 || ichain >= fMCMCNChains)
00350       return 0;
00351 
00352    // check if ipar is in range
00353    if (ipar < 0 || ipar >= fMCMCNParameters)
00354       return 0;
00355 
00356    // return component of jth point in the ith chain
00357    return fMCMCx.at(ichain *  fMCMCNParameters + ipar);
00358 }

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

Definition at line 328 of file BCEngineMCMC.cxx.

00329 {
00330    // create a new vector with the length of fMCMCNParameters
00331    std::vector <double> x;
00332 
00333    // check if ichain is in range
00334    if (ichain < 0 || ichain >= fMCMCNChains)
00335       return x;
00336 
00337    // copy the point in the ichain chain into the temporary vector
00338    for (int j = 0; j < fMCMCNParameters; ++j)
00339       x.push_back(fMCMCx.at(ichain * fMCMCNParameters + j));
00340 
00341    return x;
00342 }

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

Definition at line 179 of file BCEngineMCMC.h.

00180          { return fMCMCx; };

int BCEngineMCMC::MCMCInitialize ( )

Definition at line 1705 of file BCEngineMCMC.cxx.

01706 {
01707    // reset variables
01708    fMCMCNIterations.clear();
01709    fMCMCNIterationsConvergenceLocal.clear();
01710    fMCMCNTrialsTrue.clear();
01711    fMCMCNTrialsFalse.clear();
01712    fMCMCTrialFunctionScaleFactor.clear();
01713    fMCMCMean.clear();
01714    fMCMCVariance.clear();
01715    fMCMCMeanx.clear();
01716    fMCMCVariancex.clear();
01717    fMCMCx.clear();
01718    fMCMCLogProbx.clear();
01719    fMCMCMaximumPoints.clear();
01720    fMCMCMaximumLogProb.clear();
01721 // fMCMCRelativePrecisionModeValues.clear();
01722    fMCMCNumericalPrecisionModeValues.clear();
01723    fMCMCNIterationsConvergenceGlobal = -1;
01724    fMCMCFlagConvergenceGlobal = false;
01725    fMCMCRValueParameters.clear();
01726 
01727    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01728       if (fMCMCH1Marginalized[i])
01729          delete fMCMCH1Marginalized[i];
01730 
01731    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01732       if (fMCMCH2Marginalized[i])
01733          delete fMCMCH2Marginalized[i];
01734 
01735    // clear plots
01736    fMCMCH1Marginalized.clear();
01737    fMCMCH2Marginalized.clear();
01738 
01739 // if (fMCMCH1RValue)
01740 //    delete fMCMCH1RValue;
01741 
01742 // if (fMCMCH1Efficiency)
01743 //    delete fMCMCH1Efficiency;
01744 
01745 // free memory for vectors
01746    fMCMCNIterationsConvergenceLocal.assign(fMCMCNChains, -1);
01747    fMCMCNIterations.assign(fMCMCNChains, 0);
01748    fMCMCMean.assign(fMCMCNChains, 0);
01749    fMCMCVariance.assign(fMCMCNChains, 0);
01750    fMCMCLogProbx.assign(fMCMCNChains, -1.0);
01751    fMCMCMaximumLogProb.assign(fMCMCNChains, -1.0);
01752 
01753    fMCMCNumericalPrecisionModeValues.assign(fMCMCNParameters, 0.0);
01754 
01755    fMCMCNTrialsTrue.assign(fMCMCNChains * fMCMCNParameters, 0);
01756    fMCMCNTrialsFalse.assign(fMCMCNChains * fMCMCNParameters, 0);
01757    fMCMCMaximumPoints.assign(fMCMCNChains * fMCMCNParameters, 0.);
01758    fMCMCMeanx.assign(fMCMCNChains * fMCMCNParameters, 0);
01759    fMCMCVariancex.assign(fMCMCNChains * fMCMCNParameters, 0);
01760 
01761    fMCMCRValueParameters.assign(fMCMCNParameters, 0.);
01762 
01763    if (fMCMCTrialFunctionScaleFactorStart.size() == 0)
01764       fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fMCMCNParameters, 1.0);
01765    else
01766       for (int i = 0; i < fMCMCNChains; ++i)
01767          for (int j = 0; j < fMCMCNParameters; ++j)
01768             fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j));
01769 
01770    // set initial position
01771    for (int j = 0; j < fMCMCNChains; ++j)
01772       for (int i = 0; i < fMCMCNParameters; ++i)
01773          switch(fMCMCFlagInitialPosition)
01774          {
01775             // use the center of the region
01776             case 0 :
01777                fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01778                break;
01779             // use a random variable in the valid region
01780             case 1 :
01781                fMCMCx.push_back(fMCMCBoundaryMin[i] + fMCMCRandom -> Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01782                break;
01783             // use user-defined value
01784             case 2 :
01785                if (int(fMCMCInitialPosition.size()) == fMCMCNParameters * fMCMCNChains)
01786                   fMCMCx.push_back(fMCMCInitialPosition[j * fMCMCNParameters + i]);
01787                else
01788                   fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01789                break;
01790             // use the center of the region as default
01791             default:
01792                fMCMCx.push_back(fMCMCBoundaryMin[i] + 0.5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01793                break;
01794          }
01795 
01796    // copy the point of the first chain
01797    fMCMCxLocal.clear();
01798    for (int i = 0; i < fMCMCNParameters; ++i)
01799       fMCMCxLocal.push_back(fMCMCx[i]);
01800 
01801    // define 1-dimensional histograms for marginalization
01802    for(int i = 0; i < fMCMCNParameters; ++i)
01803    {
01804       double hmin1 = fMCMCBoundaryMin.at(i);
01805       double hmax1 = fMCMCBoundaryMax.at(i);
01806 
01807       TH1D * h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i), "", fMCMCH1NBins[i], hmin1, hmax1);
01808       fMCMCH1Marginalized.push_back(h1);
01809    }
01810 
01811    for(int i = 0; i < fMCMCNParameters; ++i)
01812       for (int k = 0; k < i; ++k)
01813       {
01814          double hmin1 = fMCMCBoundaryMin.at(k);
01815          double hmax1 = fMCMCBoundaryMax.at(k);
01816          double hmin2 = fMCMCBoundaryMin.at(i);
01817          double hmax2 = fMCMCBoundaryMax.at(i);
01818 
01819          TH2D * h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, k), "",
01820                fMCMCH1NBins[i], hmin1, hmax1,
01821                fMCMCH1NBins[k], hmin2, hmax2);
01822          fMCMCH2Marginalized.push_back(h2);
01823       }
01824 
01825    // define plot for R-value
01826 // if (fMCMCNChains > 1)
01827 // {
01828 //    fMCMCH1RValue = new TH1D("h1_rvalue", ";Iteration;R-value",
01829 //           100, 0, double(fMCMCNIterationsMax));
01830 //    fMCMCH1RValue -> SetStats(false);
01831 // }
01832 
01833 // fMCMCH1Efficiency = new TH1D("h1_efficiency", ";Chain;Efficiency",
01834 //       fMCMCNParameters, -0.5, fMCMCNParameters - 0.5);
01835 // fMCMCH1Efficiency -> SetStats(false);
01836 
01837    return 1;
01838 }

void BCEngineMCMC::MCMCInitializeMarkovChains ( )

Definition at line 1687 of file BCEngineMCMC.cxx.

01688 {
01689    // evaluate function at the starting point
01690    std::vector <double> x0;
01691 
01692    for (int j = 0; j < fMCMCNChains; ++j)
01693    {
01694       x0.clear();
01695       for (int i = 0; i < fMCMCNParameters; ++i)
01696          x0.push_back(fMCMCx[j * fMCMCNParameters + i]);
01697       fMCMCLogProbx[j] = this -> LogEval(x0);
01698    }
01699 
01700    x0.clear();
01701 }

virtual void BCEngineMCMC::MCMCIterationInterface ( ) [inline, virtual]

Reimplemented in BCIntegrate.

Definition at line 551 of file BCEngineMCMC.h.

00552          {};

int BCEngineMCMC::MCMCMetropolis ( )

Definition at line 1467 of file BCEngineMCMC.cxx.

01468 {
01469    // check if prerun has been performed
01470    if (!fMCMCFlagPreRun)
01471       this -> MCMCMetropolisPreRun();
01472 
01473    // print to screen
01474    BCLog::Out(BCLog::summary, BCLog::summary, "Run Metropolis MCMC...");
01475 
01476    if (fMCMCFlagPCA)
01477       BCLog::Out(BCLog::detail, BCLog::detail, " --> PCA switched on");
01478 
01479    // reset run statistics
01480    this -> MCMCResetRunStatistics();
01481 
01482    // perform run
01483    BCLog::Out(BCLog::summary, BCLog::summary,
01484          Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));
01485 
01486 // int counterupdate = 0;
01487 // bool convergence = false;
01488 // bool flagefficiency = false;
01489 
01490 // std::vector <double> efficiency;
01491 
01492 // for (int i = 0; i < fMCMCNParameters; ++i)
01493 //    for (int j = 0; j < fMCMCNChains; ++j)
01494 //       efficiency.push_back(0.0);
01495 
01496    int nwrite = fMCMCNIterationsRun/10;
01497    if(nwrite < 100)
01498       nwrite=100;
01499    else if(nwrite < 500)
01500       nwrite=1000;
01501    else if(nwrite < 10000)
01502       nwrite=1000;
01503    else
01504       nwrite=10000;
01505 
01506    // start the run
01507    for (int iiterations = 0; iiterations < fMCMCNIterationsRun; ++iiterations)
01508    {
01509       if ( (iiterations+1)%nwrite == 0 )
01510          BCLog::Out(BCLog::detail, BCLog::detail,
01511                Form(" --> iteration number %i (%.2f%%)", iiterations+1, (double)(iiterations+1)/(double)fMCMCNIterationsRun*100.));
01512 
01513       // if the flag is set then run over the parameters one after the other.
01514       if (fMCMCFlagOrderParameters)
01515       {
01516          // loop over parameters
01517          for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
01518             {
01519                // loop over chains
01520                for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01521                   this -> MCMCGetNewPointMetropolis(ichains, iparameters, fMCMCFlagPCA);
01522 
01523                // update statistics
01524                this -> MCMCUpdateStatistics();
01525 
01526                // do anything interface
01527                this -> MCMCIterationInterface();
01528             } // end loop over all parameters
01529       }
01530       // if the flag is not set then run over the parameters at the same time.
01531       else
01532       {
01533          // loop over chains
01534          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01535             // get new point
01536             this -> MCMCGetNewPointMetropolis(ichains, false);
01537 
01538          // update statistics
01539          this -> MCMCUpdateStatistics();
01540 
01541          // do anything interface
01542          this -> MCMCIterationInterface();
01543       }
01544    } // end run
01545 
01546    // print convergence status
01547    BCLog::Out(BCLog::summary, BCLog::summary,
01548          Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));
01549 
01550    // print modes
01551 
01552    // find global maximum
01553    double probmax = fMCMCMaximumLogProb.at(0);
01554    int probmaxindex = 0;
01555 
01556    // loop over all chains and find the maximum point
01557    for (int i = 1; i < fMCMCNChains; ++i)
01558    if (fMCMCMaximumLogProb.at(i) > probmax)
01559    {
01560       probmax = fMCMCMaximumLogProb.at(i);
01561       probmaxindex = i;
01562    }
01563 
01564    BCLog::Out(BCLog::detail, BCLog::detail, " --> Global mode from MCMC:");
01565    int ndigits = (int) log10(fMCMCNParameters);
01566    for (int i = 0; i < fMCMCNParameters; ++i)
01567       BCLog::Out(BCLog::detail, BCLog::detail,
01568             Form( TString::Format(" -->      parameter %%%di:   %%.4g", ndigits+1),
01569                   i, fMCMCMaximumPoints[probmaxindex * fMCMCNParameters + i]));
01570 
01571    // set flag
01572    fMCMCFlagPreRun = false;
01573 
01574    return 1;
01575 }

int BCEngineMCMC::MCMCMetropolisHastings ( )

Definition at line 1579 of file BCEngineMCMC.cxx.

01580 {
01581    BCLog::Out(BCLog::summary, BCLog::summary, "Run Metropolis-Hastings MCMC.");
01582 
01583    // initialize Markov chain
01584    this -> MCMCInitializeMarkovChains();
01585 
01586    // perform burn-in run
01587    BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Start burn-in run with %i iterations.", fMCMCNIterationsBurnIn));
01588    for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
01589       for (int j = 0; j < fMCMCNChains; ++j)
01590          this -> MCMCGetNewPointMetropolisHastings(j);
01591 
01592    // reset run statistics
01593    this -> MCMCResetRunStatistics();
01594 
01595    // perform run
01596    BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Perform MCMC run with %i iterations.", fMCMCNIterationsMax));
01597    for (int i = 0; i < fMCMCNIterationsMax; ++i)
01598    {
01599       for (int j = 0; j < fMCMCNChains; ++j)
01600          // get new point and increase counters
01601          this -> MCMCGetNewPointMetropolisHastings(j);
01602 
01603       // update statistics
01604       this -> MCMCUpdateStatistics();
01605 
01606       // call user interface
01607       this -> MCMCIterationInterface();
01608    }
01609 
01610    if (fMCMCNIterationsConvergenceGlobal > 0)
01611       BCLog::Out(BCLog::detail, BCLog::detail,
01612             Form(" --> Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
01613    else
01614       BCLog::Out(BCLog::detail, BCLog::detail,
01615             Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));
01616 
01617    // debug
01618    if (DEBUG)
01619    {
01620       for (int i = 0; i < fMCMCNChains; ++i)
01621       {
01622          std::cout << i << " " << fMCMCMaximumLogProb[i] << std::endl;
01623          for (int j = 0; j < fMCMCNParameters; ++j)
01624             std::cout << fMCMCMaximumPoints[i * fMCMCNParameters + j] << " ";
01625          std::cout << std::endl;
01626       }
01627    }
01628 
01629    return 1;
01630 }

int BCEngineMCMC::MCMCMetropolisPreRun ( )

Definition at line 935 of file BCEngineMCMC.cxx.

00936 {
00937    // print on screen
00938    BCLog::Out(BCLog::summary, BCLog::summary, "Pre-run Metropolis MCMC...");
00939 
00940    // initialize Markov chain
00941    this -> MCMCInitialize();
00942    this -> MCMCInitializeMarkovChains();
00943 
00944    // helper variable containing number of digits in the number of parameters
00945    int ndigits = (int)log10(fMCMCNParameters) +1;
00946    if(ndigits<4)
00947       ndigits=4;
00948 
00949    // perform burn-in run
00950    if (fMCMCNIterationsBurnIn > 0)
00951       BCLog::Out(BCLog::detail, BCLog::detail,
00952             Form(" --> Start burn-in run with %i iterations", fMCMCNIterationsBurnIn));
00953 
00954    // if the flag is set then run over the parameters one after the other.
00955    if (fMCMCFlagOrderParameters)
00956    {
00957       for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
00958          for (int j = 0; j < fMCMCNChains; ++j)
00959             for (int k = 0; k < fMCMCNParameters; ++k)
00960                this -> MCMCGetNewPointMetropolis(j, k, false);
00961    }
00962    // if the flag is not set then run over the parameters at the same time.
00963    else
00964    {
00965       for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
00966       {
00967          // loop over chains
00968          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00969             // get new point
00970             this -> MCMCGetNewPointMetropolis(ichains, false);
00971       }
00972    }
00973 
00974    // reset run statistics
00975    this -> MCMCResetRunStatistics();
00976    fMCMCNIterationsConvergenceGlobal = -1;
00977 
00978    // define data buffer for pca run
00979    double * dataall = 0;
00980    double * sum = 0;
00981    double * sum2 = 0;
00982 
00983    // allocate memory if pca is switched on
00984    if (fMCMCFlagPCA)
00985    {
00986       BCLog::Out(BCLog::detail, BCLog::detail, " --> PCA switched on");
00987 
00988       // create new TPrincipal
00989       fMCMCPCA = new TPrincipal(fMCMCNParameters, "D");
00990 
00991       // create buffer of sampled points
00992       dataall = new double[fMCMCNParameters * fMCMCNIterationsPCA];
00993 
00994       // create buffer for the sum and the sum of the squares
00995       sum = new double[fMCMCNParameters];
00996       sum2 = new double[fMCMCNParameters];
00997 
00998       // reset the buffers
00999       for (int i = 0; i < fMCMCNParameters; ++i)
01000       {
01001          sum[i]  = 0.0;
01002          sum2[i] = 0.0;
01003       }
01004 
01005       if (fMCMCFlagPCATruncate == true)
01006       {
01007          fMCMCNDimensionsPCA = 0;
01008 
01009          const double * ma = fMCMCPCA -> GetEigenValues() -> GetMatrixArray();
01010 
01011          for (int i = 0; i < fMCMCNParameters; ++i)
01012             if (ma[i]/ma[0] > fMCMCPCAMinimumRatio)
01013                fMCMCNDimensionsPCA++;
01014 
01015          BCLog::Out(BCLog::detail, BCLog::detail,
01016                Form(" --> Use %i out of %i dimensions", fMCMCNDimensionsPCA, fMCMCNParameters));
01017       }
01018       else
01019          fMCMCNDimensionsPCA = fMCMCNParameters;
01020    }
01021    else
01022       BCLog::Out(BCLog::detail, BCLog::detail, " --> No PCA run");
01023 
01024    // reset run statistics
01025    this -> MCMCResetRunStatistics();
01026    fMCMCNIterationsConvergenceGlobal = -1;
01027 
01028    // perform run
01029    BCLog::Out(BCLog::summary, BCLog::summary,
01030          Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax));
01031 
01032    bool tempflag_writetofile = fMCMCFlagWriteChainToFile;
01033 
01034    // don't write to file during pre run
01035    fMCMCFlagWriteChainToFile = false;
01036 
01037    int counter = 0;
01038    int counterupdate = 0;
01039    bool convergence = false;
01040    bool flagefficiency = false;
01041 
01042    std::vector <double> efficiency;
01043 
01044    for (int i = 0; i < fMCMCNParameters; ++i)
01045       for (int j = 0; j < fMCMCNChains; ++j)
01046          efficiency.push_back(0.);
01047 
01048    if (fMCMCFlagPCA)
01049    {
01050       fMCMCNIterationsPreRunMin = fMCMCNIterationsPCA;
01051       fMCMCNIterationsMax = fMCMCNIterationsPCA;
01052    }
01053 
01054    // run chain
01055    while (counter < fMCMCNIterationsPreRunMin || (counter < fMCMCNIterationsMax && !(convergence && flagefficiency)))
01056    {
01057       // set convergence to false by default
01058       convergence = false;
01059 
01060       // debugKK
01061       fMCMCNIterationsConvergenceGlobal = -1; 
01062 
01063       // if the flag is set then run over the parameters one after the other.
01064       if (fMCMCFlagOrderParameters)
01065       {
01066          // loop over parameters
01067          for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
01068          {
01069             // loop over chains
01070             for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01071             {
01072                this -> MCMCGetNewPointMetropolis(ichains, iparameters, false);
01073 
01074                // save point for finding the eigenvalues
01075                if (fMCMCFlagPCA)
01076                {
01077                   // create buffer for the new point
01078                   double * data = new double[fMCMCNParameters];
01079 
01080                   // copy the current point
01081                   for (int j = 0; j < fMCMCNParameters; ++j)
01082                   {
01083                      data[j] = fMCMCx[j];
01084                      dataall[ichains * fMCMCNParameters + j] = fMCMCx[j];
01085                   }
01086 
01087                   // add point to PCA object
01088                   fMCMCPCA -> AddRow(data);
01089 
01090                   // delete buffer
01091                   delete [] data;
01092 
01093                } // end if fMCMCPCAflag
01094             } // end loop over chains
01095 
01096             // search for global maximum
01097             this -> MCMCUpdateStatisticsCheckMaximum();
01098 
01099             // check convergence status
01100             // debugKK
01101             //          this -> MCMCUpdateStatisticsTestConvergenceAllChains();
01102 
01103          } // end loop over parameters
01104       } // end if fMCMCFlagOrderParameters
01105 
01106       // if the flag is not set then run over the parameters at the same time.
01107       else
01108       {
01109          // loop over chains
01110          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01111             // get new point
01112             this -> MCMCGetNewPointMetropolis(ichains, false);
01113 
01114          // save point for finding the eigenvalues
01115          if (fMCMCFlagPCA)
01116          {
01117             // create buffer for the new point
01118             double * data = new double[fMCMCNParameters];
01119 
01120             // copy the current point
01121             for (int j = 0; j < fMCMCNParameters; ++j)
01122             {
01123                // loop over chains
01124                for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01125                {
01126                   data[j] = fMCMCx[j];
01127                   dataall[ichains * fMCMCNParameters + j] = fMCMCx[j];
01128                }
01129             }
01130 
01131             // add point to PCA object
01132             fMCMCPCA -> AddRow(data);
01133 
01134             // delete buffer
01135             delete [] data;
01136          }
01137 
01138          // search for global maximum
01139          this -> MCMCUpdateStatisticsCheckMaximum();
01140 
01141          // check convergence status
01142          // debugKK
01143          //       this -> MCMCUpdateStatisticsTestConvergenceAllChains();
01144       }
01145 
01146       //-------------------------------------------
01147       // update scale factors and check convergence 
01148       //-------------------------------------------
01149       if (counterupdate % (fMCMCNIterationsUpdate*(fMCMCNParameters+1)) == 0 && counterupdate > 0 && counter >= fMCMCNIterationsPreRunMin)
01150       {
01151          // prompt status
01152          BCLog::Out(BCLog::detail, BCLog::detail,
01153                Form(" --> Iteration %i", fMCMCNIterations[0] / fMCMCNParameters));
01154 
01155          bool rvalues_ok = true;
01156          static bool has_converged = false;
01157 
01158          // -----------------------------
01159          // check convergence status
01160          // -----------------------------
01161          // debugKK
01162          fMCMCNIterationsConvergenceGlobal = -1; 
01163 
01164          this -> MCMCUpdateStatisticsTestConvergenceAllChains();
01165 
01166          // debugKK
01167          //       std::cout << " HERE " << fMCMCRValueParameters[0] << fMCMCNTrialsTrue[0] << " " << fMCMCNTrialsFalse[0] << std::endl; 
01168 
01169          // check convergence
01170          if (fMCMCNIterationsConvergenceGlobal > 0)
01171             convergence = true;
01172          
01173          // print convergence status: 
01174          if (convergence && fMCMCNChains > 1)
01175             BCLog::Out(BCLog::detail, BCLog::detail,
01176                             Form("     * Convergence status: Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
01177          else if (!convergence && fMCMCNChains > 1)
01178             {
01179                BCLog::Out(BCLog::detail, BCLog::detail,
01180                             Form("     * Convergence status: Set of %i Markov chains did not converge after %i iterations.", fMCMCNChains, counter));
01181                
01182                BCLog::Out(BCLog::detail, BCLog::detail, "       - R-Values:");
01183                for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
01184                   {
01185                      if(fabs(fMCMCRValueParameters[iparameter]-1.) < fMCMCRValueParametersCriterion)
01186                         BCLog::Out(BCLog::detail, BCLog::detail,
01187                                         Form( TString::Format("         parameter %%%di :  %%.06f",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
01188                      else
01189                         {
01190                            BCLog::Out(BCLog::detail, BCLog::detail,
01191                                            Form( TString::Format("         parameter %%%di :  %%.06f <--",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
01192                            rvalues_ok = false;
01193                         }
01194                   }
01195                if(fabs(fMCMCRValue-1.) < fMCMCRValueCriterion)
01196                   BCLog::Out(BCLog::detail, BCLog::detail, Form("         log-likelihood :  %.06f", fMCMCRValue));
01197                else
01198                   {
01199                      BCLog::Out(BCLog::detail, BCLog::detail, Form("         log-likelihood :  %.06f <--", fMCMCRValue));
01200                      rvalues_ok = false;
01201                   }
01202             }
01203 
01204          if(!has_converged)
01205             if(rvalues_ok)
01206                has_converged = true;
01207 
01208          // -----------------------------
01209          // check efficiency status
01210          // -----------------------------
01211 
01212          // set flag
01213          flagefficiency = true;
01214 
01215          bool flagprintefficiency = true; 
01216 
01217          // loop over chains
01218          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01219          {
01220             // loop over parameters
01221             for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
01222             {
01223                // global index of the parameter (throughout all the chains)
01224                int ipc = ichains * fMCMCNParameters + iparameter;
01225 
01226                // calculate efficiency
01227                efficiency[ipc] = double(fMCMCNTrialsTrue[ipc]) / double(fMCMCNTrialsTrue[ipc] + fMCMCNTrialsFalse[ipc]);
01228 
01229                // adjust scale factors if efficiency is too low
01230                if (efficiency[ipc] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[ipc] > .01)
01231                {
01232                   if (flagprintefficiency)
01233                      {
01234                         BCLog::Out(BCLog::detail, BCLog::detail,
01235                                         Form("     * Efficiency status: Efficiencies not within pre-defined range."));
01236                         BCLog::Out(BCLog::detail, BCLog::detail, 
01237                                         Form("       - Efficiencies:"));
01238                         flagprintefficiency = false; 
01239                      }
01240 
01241                   double fscale=2.;
01242                   if(has_converged && fMCMCEfficiencyMin/efficiency[ipc] > 2.)
01243                      fscale = 4.;
01244                   fMCMCTrialFunctionScaleFactor[ipc] /= fscale;
01245 
01246                   BCLog::Out(BCLog::detail, BCLog::detail,
01247                         Form("         Efficiency of parameter %i dropped below %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
01248                               iparameter, 100. * fMCMCEfficiencyMin, 100. * efficiency[ipc], ichains, fMCMCTrialFunctionScaleFactor[ipc]));
01249                }
01250 
01251                // adjust scale factors if efficiency is too high
01252                else if (efficiency[ipc] > fMCMCEfficiencyMax && (fMCMCTrialFunctionScaleFactor[ipc] < 1. || (fMCMCFlagPCA && fMCMCTrialFunctionScaleFactor[ipc] < 10.)))
01253                {
01254                   if (flagprintefficiency)
01255                      {
01256                         BCLog::Out(BCLog::detail, BCLog::detail,
01257                                         Form("   * Efficiency status: Efficiencies not within pre-defined ranges."));
01258                         BCLog::Out(BCLog::detail, BCLog::detail, 
01259                                         Form("     - Efficiencies:"));
01260                         flagprintefficiency = false; 
01261                      }
01262 
01263                   fMCMCTrialFunctionScaleFactor[ipc] *= 2.;
01264 
01265                   BCLog::Out(BCLog::detail, BCLog::detail,
01266                                   Form("         Efficiency of parameter %i above %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
01267                                           iparameter, 100.0 * fMCMCEfficiencyMax, 100.0 * efficiency[ipc], ichains, fMCMCTrialFunctionScaleFactor[ipc]));
01268                }
01269 
01270                // check flag
01271                if ((efficiency[ipc] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[ipc] > .01)
01272                      || (efficiency[ipc] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[ipc] < 1.))
01273                   flagefficiency = false;
01274 
01275                // reset counters
01276                counterupdate = 0;
01277                fMCMCNTrialsTrue[ipc] = 0;
01278                fMCMCNTrialsFalse[ipc] = 0;
01279             } // end of running over all parameters 
01280 
01281             // reset counters
01282             fMCMCMean[ichains] = 0;
01283             fMCMCVariance[ichains] = 0;
01284          } // end of running over all chains 
01285          
01286          // print to scrren 
01287          if (flagefficiency)
01288             BCLog::Out(BCLog::detail, BCLog::detail,
01289                             Form("     * Efficiency status: Efficiencies within pre-defined ranges.")); 
01290          
01291       } // end if update scale factors and check convergence
01292 
01293       // increase counter
01294       counter++;
01295       counterupdate++;
01296 
01297       // debugKK
01298 //       // check convergence
01299 //       if (fMCMCNIterationsConvergenceGlobal > 0)
01300 //          convergence = true;
01301    } // end of running 
01302 
01303    // ---------------
01304    // after chain ran
01305    // ---------------
01306 
01307    // define convergence status
01308    if (fMCMCNIterationsConvergenceGlobal > 0)
01309       fMCMCFlagConvergenceGlobal = true;
01310    else
01311       fMCMCFlagConvergenceGlobal = false;
01312 
01313    // print convergence status
01314    if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && !flagefficiency)
01315       BCLog::Out(BCLog::summary, BCLog::summary,
01316             Form(" --> Set of %i Markov chains converged within %i iterations but could not adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
01317 
01318    else if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && flagefficiency)
01319       BCLog::Out(BCLog::summary, BCLog::summary,
01320             Form(" --> Set of %i Markov chains converged within %i iterations and could adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
01321 
01322    else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && flagefficiency)
01323       BCLog::Out(BCLog::summary, BCLog::summary,
01324              Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));
01325 
01326    else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && !flagefficiency)
01327       BCLog::Out(BCLog::summary, BCLog::summary,
01328              Form(" --> Set of %i Markov chains did not converge within %i iterations and could not adjust scales.", fMCMCNChains, fMCMCNIterationsMax));
01329 
01330    else if(fMCMCNChains == 1)
01331       BCLog::Out(BCLog::summary, BCLog::summary,
01332              " --> No convergence criterion for a single chain defined.");
01333 
01334    else
01335       BCLog::Out(BCLog::summary, BCLog::summary,
01336             " --> Only one Markov chain. No global convergence criterion defined.");
01337 
01338    BCLog::Out(BCLog::summary, BCLog::summary,
01339          Form(" --> Markov chains ran for %i iterations.", counter));
01340 
01341 
01342    // print efficiencies
01343    std::vector <double> efficiencies;
01344 
01345    for (int i = 0; i < fMCMCNParameters; ++i)
01346       efficiencies.push_back(0.);
01347 
01348    BCLog::Out(BCLog::detail, BCLog::detail, " --> Average efficiencies:");
01349    for (int i = 0; i < fMCMCNParameters; ++i)
01350    {
01351       for (int j = 0; j < fMCMCNChains; ++j)
01352          efficiencies[i] += efficiency[j * fMCMCNParameters + i] / double(fMCMCNChains);
01353 
01354       BCLog::Out(BCLog::detail, BCLog::detail,
01355             Form( TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * efficiencies[i]));
01356    }
01357 
01358 
01359    // print scale factors
01360    std::vector <double> scalefactors;
01361 
01362    for (int i = 0; i < fMCMCNParameters; ++i)
01363       scalefactors.push_back(0.0);
01364 
01365    BCLog::Out(BCLog::detail, BCLog::detail, " --> Average scale factors:");
01366    for (int i = 0; i < fMCMCNParameters; ++i)
01367    {
01368       for (int j = 0; j < fMCMCNChains; ++j)
01369          scalefactors[i] += fMCMCTrialFunctionScaleFactor[j * fMCMCNParameters + i] / double(fMCMCNChains);
01370 
01371       BCLog::Out(BCLog::detail, BCLog::detail,
01372             Form( TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * scalefactors[i]));
01373    }
01374 
01375 
01376    // perform PCA analysis
01377    if (fMCMCFlagPCA)
01378    {
01379 
01380       // calculate eigenvalues and vectors
01381       fMCMCPCA -> MakePrincipals();
01382 
01383       // re-run over data points to gain a measure for the spread of the variables
01384       for (int i = 0; i < fMCMCNIterationsPCA; ++i)
01385       {
01386          double * data = new double[fMCMCNParameters];
01387          double * p    = new double[fMCMCNParameters];
01388 
01389          for (int j = 0; j < fMCMCNParameters; ++j)
01390             data[j] = dataall[i * fMCMCNParameters + j];
01391 
01392          fMCMCPCA -> X2P(data, p);
01393 
01394          for (int j = 0; j < fMCMCNParameters; ++j)
01395          {
01396             sum[j]  += p[j];
01397             sum2[j] += p[j] * p[j];
01398          }
01399 
01400          delete [] data;
01401          delete [] p;
01402       }
01403 
01404       delete [] dataall;
01405 
01406       fMCMCPCAMean.clear();
01407       fMCMCPCAVariance.clear();
01408 
01409       // calculate mean and variance
01410       for (int j = 0; j < fMCMCNParameters; ++j)
01411       {
01412          fMCMCPCAMean.push_back(sum[j] / double(fMCMCNIterationsPCA) / double(fMCMCNChains));
01413          fMCMCPCAVariance.push_back(sum2[j] / double(fMCMCNIterationsPCA) / double(fMCMCNChains) - fMCMCPCAMean[j] * fMCMCPCAMean[j]);
01414       }
01415 
01416       // check if all eigenvalues are found
01417       int neigenvalues = fMCMCPCA -> GetEigenValues() -> GetNoElements();
01418 
01419       const double * eigenv = fMCMCPCA -> GetEigenValues() -> GetMatrixArray();
01420 
01421       bool flageigenvalues = true;
01422 
01423       for (int i = 0; i < neigenvalues; ++i)
01424          if (isnan(eigenv[i]))
01425             flageigenvalues = false;
01426 
01427       // print on screen
01428       if (flageigenvalues)
01429          BCLog::Out(BCLog::detail, BCLog::detail, " --> PCA : All eigenvalues ok.");
01430       else
01431       {
01432          BCLog::Out(BCLog::detail, BCLog::detail, " --> PCA : Not all eigenvalues ok. Don't use PCA.");
01433          fMCMCFlagPCA = false;
01434       }
01435 
01436       // reset scale factors
01437       for (int i = 0; i < fMCMCNParameters; ++i)
01438          for (int j = 0; j < fMCMCNChains; ++j)
01439             fMCMCTrialFunctionScaleFactor[j * fMCMCNParameters + i] = 1.0;
01440 
01441       if (DEBUG)
01442       {
01443          fMCMCPCA -> Print("MSEV");
01444 
01445          for (int j = 0; j < fMCMCNParameters; ++j)
01446             std::cout << fMCMCPCAMean.at(j) << " " << sqrt(fMCMCPCAVariance.at(j)) << std::endl;
01447       }
01448 
01449    }
01450 
01451    // fill efficiency plot
01452    // for (int i = 0; i < fMCMCNParameters; ++i)
01453    //    fMCMCH1Efficiency -> SetBinContent(i + 1, efficiencies[i]);
01454 
01455    // reset flag
01456    fMCMCFlagWriteChainToFile = tempflag_writetofile;
01457 
01458    // set flag
01459    fMCMCFlagPreRun = true;
01460 
01461    return 1;
01462 
01463 }

void BCEngineMCMC::MCMCResetRunStatistics ( )

Definition at line 1634 of file BCEngineMCMC.cxx.

01635 {
01636    // debugKK -> this is set explicitly
01637    // fMCMCNIterationsConvergenceGlobal = -1;
01638 
01639    for (int j = 0; j < fMCMCNChains; ++j)
01640    {
01641       fMCMCNIterations[j]  = 0;
01642       fMCMCNTrialsTrue[j]  = 0;
01643       fMCMCNTrialsFalse[j] = 0;
01644       fMCMCMean[j]         = 0;
01645       fMCMCVariance[j]     = 0;
01646 
01647       for (int k = 0; k < fMCMCNParameters; ++k)
01648       {
01649          fMCMCNTrialsTrue[j * fMCMCNParameters + k]  = 0;
01650          fMCMCNTrialsFalse[j * fMCMCNParameters + k] = 0;
01651       }
01652    }
01653 
01654    // reset marginalized distributions
01655    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01656       fMCMCH1Marginalized[i] -> Reset();
01657 
01658    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01659       fMCMCH2Marginalized[i] -> Reset();
01660 
01661 // if (fMCMCH1RValue)
01662 //   fMCMCH1RValue -> Reset();
01663 
01664 // if (fMCMCH1Efficiency)
01665 //   fMCMCH1Efficiency -> Reset();
01666 
01667    fMCMCRValue = 100;
01668 }

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

Definition at line 373 of file BCEngineMCMC.h.

00374          { fMCMCFlagInitialPosition = flag; };

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

Definition at line 379 of file BCEngineMCMC.h.

00380       { fMCMCFlagOrderParameters = flag; };

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

Definition at line 399 of file BCEngineMCMC.h.

00400          { fMCMCFlagPCA = flag; };

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

Definition at line 409 of file BCEngineMCMC.h.

00410          { fMCMCFlagPCATruncate = flag; };

void BCEngineMCMC::MCMCSetInitialPosition	(	int	chain,
		std::vector< double >	x0
	)

Definition at line 165 of file BCEngineMCMC.cxx.

00166 {
00167    // check consistency
00168    if (chain >= fMCMCNChains)
00169    {
00170       fMCMCInitialPosition.clear();
00171       return;
00172    }
00173 
00174    if (int(x0.size()) == fMCMCNParameters)
00175    {
00176       fMCMCInitialPosition.clear();
00177       for (int j = 0; j < fMCMCNChains; ++j)
00178          for (int i = 0; i < fMCMCNParameters; ++i)
00179             fMCMCInitialPosition.push_back(x0.at(i));
00180    }
00181    else
00182       return;
00183 
00184    bool flagin = true;
00185    for (int j = 0; j < fMCMCNChains; ++j)
00186       for (int i = 0; i < fMCMCNParameters; ++i)
00187          if (fMCMCInitialPosition[j * fMCMCNParameters + i] < fMCMCBoundaryMin[i] ||
00188                fMCMCInitialPosition[j * fMCMCNParameters + i] > fMCMCBoundaryMax[i])
00189             flagin = false;
00190 
00191    if (flagin == false)
00192    {
00193       fMCMCInitialPosition.clear();
00194       return;
00195    }
00196    // copy this intial position into the particular Markov chain
00197    else
00198       for (int i = 0; i < fMCMCNParameters; ++i)
00199          fMCMCInitialPosition[chain * fMCMCNParameters + i] = x0.at(i);
00200 
00201    // use this intial position for the Markov chain
00202    this -> MCMCSetFlagInitialPosition(2);
00203 }

void BCEngineMCMC::MCMCSetInitialPosition ( std::vector< double > x0 )

Definition at line 207 of file BCEngineMCMC.cxx.

00208 {
00209    for (int i = 0; i < fMCMCNChains; ++i)
00210       this -> MCMCSetInitialPosition(i, x0);
00211 }

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

Definition at line 215 of file BCEngineMCMC.cxx.

00216 {
00217    if (int(fMCMCInitialPosition.size()) != (fMCMCNChains * fMCMCNParameters))
00218       return;
00219 
00220    if (int(x0s.size()) != (fMCMCNChains * fMCMCNParameters))
00221       return;
00222 
00223    // copy these intial positions into the Markov chains
00224    for (int i = 0; i < (fMCMCNChains * fMCMCNParameters); ++i)
00225       fMCMCInitialPosition[i] = x0s.at(i);
00226 
00227    // use these intial positions for the Markov chain
00228    this -> MCMCSetFlagInitialPosition(2);
00229 }

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

Definition at line 339 of file BCEngineMCMC.h.

00340          { fMCMCFlagIterationsAuto = flag; };

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

Definition at line 233 of file BCEngineMCMC.cxx.

00234 {
00235 // clear vector
00236    fMCMCTrees.clear();
00237 
00238    // copy tree
00239    for (int i = 0; i < int(trees.size()); ++i)
00240       fMCMCTrees.push_back(trees[i]);
00241 }

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

Definition at line 349 of file BCEngineMCMC.h.

00350          { fMCMCEfficiencyMax = efficiency; };

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

Definition at line 344 of file BCEngineMCMC.h.

00345          { fMCMCEfficiencyMin = efficiency; };

void BCEngineMCMC::MCMCSetNChains ( int n )

Definition at line 155 of file BCEngineMCMC.cxx.

00156 {
00157    fMCMCNChains = n;
00158 
00159    // re-initialize
00160    this -> MCMCInitialize();
00161 }

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

Definition at line 325 of file BCEngineMCMC.h.

00326          { fMCMCNIterationsBurnIn = n; };

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

Definition at line 315 of file BCEngineMCMC.h.

00316          { fMCMCNIterationsMax = n; };

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

Definition at line 330 of file BCEngineMCMC.h.

00331          { fMCMCNIterationsPCA = n; };

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

Definition at line 320 of file BCEngineMCMC.h.

00321          { fMCMCNIterationsRun = n; };

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

Definition at line 333 of file BCEngineMCMC.h.

00334       { fMCMCNIterationsUpdate = n; };

void BCEngineMCMC::MCMCSetPCAMinimumRatio ( double ratio ) [inline]

Definition at line 415 of file BCEngineMCMC.h.

00416          { fMCMCPCAMinimumRatio = ratio; };

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

Definition at line 384 of file BCEngineMCMC.h.

00385          { fMCMCRValueCriterion = r; };

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

Definition at line 389 of file BCEngineMCMC.h.

00390          { fMCMCRValueParametersCriterion = r; };

void BCEngineMCMC::MCMCSetTrialFunctionScale ( double scale ) [inline]

Definition at line 299 of file BCEngineMCMC.h.

00300          { fMCMCTrialFunctionScale = scale; };

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

Definition at line 302 of file BCEngineMCMC.h.

00303          { fMCMCTrialFunctionScaleFactorStart = scale; };

void BCEngineMCMC::MCMCSetValuesDefault ( )

Set the default values for the MCMC chain.

Definition at line 1890 of file BCEngineMCMC.cxx.

01891 {
01892    this -> MCMCSetValuesDetail();
01893 }

void BCEngineMCMC::MCMCSetValuesDetail ( )

Set the values for a detailed MCMC run.

Definition at line 1923 of file BCEngineMCMC.cxx.

01924 {
01925    fMCMCNChains              = 5;
01926    fMCMCNIterationsMax       = 1000000;
01927    fMCMCNIterationsRun       = 100000;
01928    fMCMCNIterationsBurnIn    = 0;
01929    fMCMCNIterationsPCA       = 10000;
01930    fMCMCNIterationsPreRunMin = 500;
01931    fMCMCFlagIterationsAuto   = false;
01932    fMCMCTrialFunctionScale   = 1.0;
01933    fMCMCFlagInitialPosition  = 1;
01934    fMCMCRValueCriterion      = 0.1;
01935    fMCMCRValueParametersCriterion = 0.1;
01936    fMCMCNIterationsConvergenceGlobal = -1;
01937    fMCMCFlagConvergenceGlobal = false;
01938    fMCMCRValue               = 100;
01939    fMCMCFlagPCA              = false;
01940    fMCMCFlagPCATruncate      = false;
01941    fMCMCPCA                  = 0;
01942    fMCMCPCAMinimumRatio      = 1e-7;
01943    fMCMCNIterationsUpdate    = 1000;
01944    fMCMCFlagOrderParameters  = true;
01945 }

void BCEngineMCMC::MCMCSetValuesQuick ( )

Set the values for a quick MCMC run.

Definition at line 1897 of file BCEngineMCMC.cxx.

01898 {
01899    fMCMCNChains              = 1;
01900    fMCMCNIterationsMax       = 1000;
01901    fMCMCNIterationsRun       = 10000;
01902    fMCMCNIterationsBurnIn    = 0;
01903    fMCMCNIterationsPCA       = 0;
01904    fMCMCNIterationsPreRunMin = 0;
01905    fMCMCFlagIterationsAuto   = false;
01906    fMCMCTrialFunctionScale   = 1.0;
01907    fMCMCFlagInitialPosition  = 1;
01908    fMCMCRValueCriterion      = 0.1;
01909    fMCMCRValueParametersCriterion = 0.1;
01910    fMCMCNIterationsConvergenceGlobal = -1;
01911    fMCMCFlagConvergenceGlobal = false;
01912    fMCMCRValue               = 100;
01913    fMCMCFlagPCA              = false;
01914    fMCMCFlagPCATruncate      = false;
01915    fMCMCPCA                  = 0;
01916    fMCMCPCAMinimumRatio      = 1e-7;
01917    fMCMCNIterationsUpdate    = 1000;
01918    fMCMCFlagOrderParameters  = true;
01919 }

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

Definition at line 354 of file BCEngineMCMC.h.

00355          { fMCMCFlagWriteChainToFile = flag; };

void BCEngineMCMC::MCMCTrialFunction ( std::vector< double > & x )

Definition at line 250 of file BCEngineMCMC.cxx.

00251 {
00252    // use uniform distribution for now
00253    for (int i = 0; i < fMCMCNParameters; ++i)
00254       x[i] = fMCMCTrialFunctionScale * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00255 }

void BCEngineMCMC::MCMCTrialFunctionAuto ( std::vector< double > & x )

Definition at line 285 of file BCEngineMCMC.cxx.

00286 {
00287    // use uniform distribution for now
00288    for (int i = 0; i < fMCMCNParameters; ++i)
00289       x[i] = fMCMCRandom -> Gaus(fMCMCPCAMean[i], sqrt(fMCMCPCAVariance[i]));
00290 }

double BCEngineMCMC::MCMCTrialFunctionIndependent	(	std::vector< double > &	xnew,
		std::vector< double > &	xold,
		bool	newpoint
	)

Definition at line 269 of file BCEngineMCMC.cxx.

00270 {
00271    // use uniform distribution for now
00272    if (newpoint)
00273       for (int i = 0; i < fMCMCNParameters; ++i)
00274          xnew[i] = fMCMCRandom -> Rndm() * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));
00275 
00276    double prob = 1.0;
00277    for (int i = 0; i < fMCMCNParameters; ++i)
00278       prob /= fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i);
00279 
00280    return prob;
00281 }

void BCEngineMCMC::MCMCTrialFunctionSingle	(	int	ichain,
		int	iparameter,
		std::vector< double > &	x
	)

Definition at line 259 of file BCEngineMCMC.cxx.

00260 {
00261    // no check of range for performance reasons
00262 
00263    // use uniform distribution
00264    x[iparameter] = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00265 }

void BCEngineMCMC::MCMCUpdateStatistics ( )

Definition at line 904 of file BCEngineMCMC.cxx.

00905 {
00906    // check if maximum is reached
00907    this -> MCMCUpdateStatisticsCheckMaximum();
00908 
00909    // fill histograms of marginalized distributions
00910    this -> MCMCUpdateStatisticsFillHistograms();
00911 
00912    // test convergence of single Markov chains
00913    //  --> not implemented yet
00914 
00915    // test convergence of all Markov chains
00916    // debugKK -> this will be called explicitly
00917    // this -> MCMCUpdateStatisticsTestConvergenceAllChains();
00918 
00919    // fill Markov chains
00920    if (fMCMCFlagWriteChainToFile)
00921       this -> MCMCUpdateStatisticsWriteChains();
00922 }

void BCEngineMCMC::MCMCUpdateStatisticsCheckMaximum ( )

Definition at line 743 of file BCEngineMCMC.cxx.

00744 {
00745    // loop over all chains
00746    for (int i = 0; i < fMCMCNChains; ++i)
00747    {
00748       if (fMCMCLogProbx[i] > fMCMCMaximumLogProb[i] || fMCMCNIterations[i] == 1)
00749       {
00750          fMCMCMaximumLogProb[i] = fMCMCLogProbx[i];
00751          for (int j = 0; j < fMCMCNParameters; ++j)
00752             fMCMCMaximumPoints[i * fMCMCNParameters + j] = fMCMCx[i * fMCMCNParameters + j];
00753       }
00754    }
00755 }

void BCEngineMCMC::MCMCUpdateStatisticsFillHistograms ( )

Definition at line 759 of file BCEngineMCMC.cxx.

00760 {
00761    // loop over chains
00762    for (int i = 0; i < fMCMCNChains; ++i)
00763    {
00764       // fill each 1-dimensional histogram
00765       for (int j = 0; j < fMCMCNParameters; ++j)
00766          fMCMCH1Marginalized[j] -> Fill(fMCMCx[i * fMCMCNParameters + j]);
00767 
00768       // fill each 2-dimensional histogram
00769       int counter = 0;
00770 
00771       for (int j = 0; j < fMCMCNParameters; ++j)
00772          for (int k = 0; k < j; ++k)
00773          {
00774             fMCMCH2Marginalized[counter] -> Fill(
00775                   fMCMCx[i * fMCMCNParameters + k],
00776                   fMCMCx[i * fMCMCNParameters + j]);
00777             counter ++;
00778          }
00779    }
00780 }

void BCEngineMCMC::MCMCUpdateStatisticsModeConvergence ( )

Definition at line 705 of file BCEngineMCMC.cxx.

00706 {
00707    double * mode_minimum = new double[fMCMCNParameters];
00708    double * mode_maximum = new double[fMCMCNParameters];
00709    double * mode_average = new double[fMCMCNParameters];
00710 
00711    // set initial values
00712    for (int j = 0; j < fMCMCNParameters; ++j)
00713    {
00714       mode_minimum[j] = fMCMCMaximumPoints[j];
00715       mode_maximum[j] = fMCMCMaximumPoints[j];
00716       mode_average[j] = 0;
00717    }
00718 
00719    // calculate the maximum difference in each dimension
00720    for (int i = 0; i < fMCMCNChains; ++i)
00721       for (int j = 0; j < fMCMCNParameters; ++j)
00722       {
00723          if (fMCMCMaximumPoints[i * fMCMCNParameters + j] < mode_minimum[j])
00724             mode_minimum[j] = fMCMCMaximumPoints[i * fMCMCNParameters + j];
00725 
00726          if (fMCMCMaximumPoints[i * fMCMCNParameters + j] > mode_maximum[j])
00727             mode_maximum[j] = fMCMCMaximumPoints[i * fMCMCNParameters + j];
00728 
00729          mode_average[j] += fMCMCMaximumPoints[i * fMCMCNParameters + j] / double(fMCMCNChains);
00730       }
00731 
00732    for (int j = 0; j < fMCMCNParameters; ++j)
00733       fMCMCNumericalPrecisionModeValues[j] = (mode_maximum[j] - mode_minimum[j]);
00734 //    fMCMCRelativePrecisionModeValues[j] = (mode_maximum[j] - mode_minimum[j]) / mode_average[j];
00735 
00736    delete [] mode_minimum;
00737    delete [] mode_maximum;
00738    delete [] mode_average;
00739 }

void BCEngineMCMC::MCMCUpdateStatisticsTestConvergenceAllChains ( )

Definition at line 784 of file BCEngineMCMC.cxx.

00785 {
00786   // calculate number of entries in this part of the chain   
00787   int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0]; 
00788 
00789   // do not evaluate the convergence criterion if the numbers of
00790   // elements is too small.
00791    //  if (npoints < fMCMCNIterationsPreRunMin) 
00792    //    {
00793    //      fMCMCNIterationsConvergenceGlobal = -1;      
00794    //      return; 
00795    //    }
00796   
00797   if (fMCMCNChains > 1 && npoints > 1)
00798     {
00799       // define flag for convergence       
00800       bool flag_convergence = true; 
00801       
00802       // loop over parameters       
00803       for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00804             {
00805                double sum = 0; 
00806                double sum2 = 0; 
00807                double sumv = 0; 
00808      
00809                // loop over chains              
00810                for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00811                   {
00812                      // get parameter index
00813                      int index = ichains * fMCMCNParameters + iparameters; 
00814                      
00815                      // calculate mean value of each parameter in the chain for this part                      
00816                      fMCMCMeanx[index] += (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints);   
00817                      
00818                      // calculate variance of each chain for this part                    
00819                      fMCMCVariancex[index] = (1.0 - 1/double(npoints)) * fMCMCVariancex[index] 
00820                         + (fMCMCx[index] - fMCMCMeanx[index]) * (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints - 1); 
00821                      
00822                      sum  += fMCMCMeanx[index]; 
00823                      sum2 += fMCMCMeanx[index] * fMCMCMeanx[index]; 
00824                      sumv += fMCMCVariancex[index]; 
00825                   }
00826                
00827                // calculate r-value for each parameter               
00828                double mean = sum / double(fMCMCNChains); 
00829                double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints); 
00830                double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains); 
00831                
00832                double r = 100.0; 
00833                
00834                if (W > 0)
00835                   {
00836                      r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W); 
00837                      
00838                      fMCMCRValueParameters[iparameters] = r; 
00839                   }
00840                // debugKK
00841                //             std::cout << " here : " << W << " " << fMCMCRValueParameters[iparameters] << std::endl; 
00842                
00843                // set flag to false if convergence criterion is not fulfilled for the parameter                
00844                if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
00845                   flag_convergence = false; 
00846             }
00847       
00848       // convergence criterion applied on the log-likelihood         
00849       double sum = 0; 
00850       double sum2 = 0; 
00851       double sumv = 0; 
00852       
00853       // loop over chains       
00854       for (int i = 0; i < fMCMCNChains; ++i)
00855             {
00856                // calculate mean value of each chain for this part               
00857                fMCMCMean[i] += (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints); 
00858                
00859                // calculate variance of each chain for this part              
00860                if (npoints > 1) 
00861                   fMCMCVariance[i] = (1.0 - 1/double(npoints)) * fMCMCVariance[i] 
00862                      + (fMCMCLogProbx[i] - fMCMCMean[i]) * (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints - 1); 
00863                
00864                sum  += fMCMCMean[i]; 
00865                sum2 += fMCMCMean[i] * fMCMCMean[i]; ; 
00866                sumv += fMCMCVariance[i]; 
00867             }
00868       
00869       // calculate r-value for log-likelihood   
00870       double mean = sum / double(fMCMCNChains); 
00871       double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints); 
00872       double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);       
00873       double r = 100.0; 
00874       
00875       if (W > 0)
00876             {
00877                r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W); 
00878                
00879                fMCMCRValue = r; 
00880             }
00881          
00882       // set flag to false if convergence criterion is not fulfilled for the log-likelihood        
00883       if (! ((fMCMCRValue-1.0) < fMCMCRValueCriterion))
00884             flag_convergence = false; 
00885       
00886       // remember number of iterations needed to converge       
00887       if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
00888             fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters; 
00889     }
00890   
00891 }

void BCEngineMCMC::MCMCUpdateStatisticsWriteChains ( )

Definition at line 895 of file BCEngineMCMC.cxx.

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

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

Defaut assignment operator

Definition at line 108 of file BCEngineMCMC.cxx.

00109 {
00110    if (this != &enginemcmc)
00111       enginemcmc.Copy(* this);
00112 
00113    return * this;
00114 }

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

Definition at line 1860 of file BCEngineMCMC.cxx.

01861 {
01862    int counter = 0;
01863    int index = 0;
01864 
01865    // search for correct combination
01866    for(int i = 0; i < fMCMCNParameters; i++)
01867       for (int j = 0; j < i; ++j)
01868       {
01869          if(j == index1 && i == index2)
01870             index = counter;
01871          counter++;
01872       }
01873 
01874    if((int)fMCMCH2Marginalized.size()<=index)
01875       return 0;
01876 
01877    if(h==0)
01878       return 0;
01879 
01880    if((int)fMCMCH2Marginalized.size()==index)
01881       fMCMCH2Marginalized.push_back(h);
01882    else
01883       fMCMCH2Marginalized[index]=h;
01884 
01885    return index;
01886 }

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

Definition at line 1842 of file BCEngineMCMC.cxx.

01843 {
01844    if((int)fMCMCH1Marginalized.size()<=index)
01845       return 0;
01846 
01847    if(h==0)
01848       return 0;
01849 
01850    if((int)fMCMCH1Marginalized.size()==index)
01851       fMCMCH1Marginalized.push_back(h);
01852    else
01853       fMCMCH1Marginalized[index]=h;
01854 
01855    return index;
01856 }