BCEngineMCMC Class Reference

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

#include <BCEngineMCMC.h>

List of all members.

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 MCMCSetInitialPositions (std::vector< std::vector< double > > x0s)

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 MCMCSetNIterationsPreRunMin (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

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.

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    fMCMCFlagInitialPosition  = 1;
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 206 of file BCEngineMCMC.cxx.

00207 {}

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

Reimplemented in BCIntegrate, and BCModel.

Definition at line 881 of file BCEngineMCMC.cxx.

00882 {
00883    // test function for now
00884    // this will be overloaded by the user
00885    return 0.0;
00886 }

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

Definition at line 1627 of file BCEngineMCMC.cxx.

01628 {
01629    // add the boundaries to the corresponding vectors
01630    fMCMCBoundaryMin.push_back(min);
01631    fMCMCBoundaryMax.push_back(max);
01632 
01633    // increase the number of parameters by one
01634    fMCMCNParameters++;
01635 
01636    // return the number of parameters
01637    return fMCMCNParameters;
01638 }

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

00323 {
00324    // check if ichain is in range
00325    if (ichain < 0 || ichain >= fMCMCNChains)
00326       return -1;
00327 
00328    // return log of the probability at the current point in the ith chain
00329    return fMCMCLogProbx.at(ichain);
00330 }

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

00479 {
00480    // calculate index
00481    int index = chain * fMCMCNParameters;
00482 
00483    // increase counter
00484    fMCMCNIterations[chain]++;
00485 
00486    // get proposal point
00487    if (!this -> MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal, pca))
00488       return false;
00489 
00490    // calculate probabilities of the old and new points
00491    double p0 = fMCMCLogProbx[chain];
00492    double p1 = this -> LogEval(fMCMCxLocal);
00493 
00494    // flag for accept
00495    bool accept = false;
00496 
00497    // if the new point is more probable, keep it ...
00498    if (p1 >= p0)
00499       accept = true;
00500    // ... or else throw dice.
00501    else
00502    {
00503       double r = log(fMCMCRandom -> Rndm());
00504 
00505       if(r < p1 - p0)
00506          accept = true;
00507    }
00508 
00509    // fill the new point
00510    if(accept)
00511    {
00512       // increase counter
00513       fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++;
00514 
00515       // copy the point
00516       for(int i = 0; i < fMCMCNParameters; ++i)
00517       {
00518          // save the point
00519          fMCMCx[index + i] = fMCMCxLocal[i];
00520 
00521          // save the probability of the point
00522          fMCMCLogProbx[chain] = p1;
00523       }
00524    }
00525    else
00526    {
00527       // increase counter
00528       fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00529    }
00530 
00531    return accept;
00532 }

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

Definition at line 536 of file BCEngineMCMC.cxx.

00537 {
00538    // calculate index
00539    int index = chain * fMCMCNParameters;
00540 
00541    // increase counter
00542    fMCMCNIterations[chain]++;
00543 
00544    // get proposal point
00545    if (!this -> MCMCGetProposalPointMetropolis(chain, fMCMCxLocal, pca))
00546       return false;
00547 
00548    // calculate probabilities of the old and new points
00549    double p0 = fMCMCLogProbx[chain];
00550    double p1 = this -> LogEval(fMCMCxLocal);
00551 
00552    // flag for accept
00553    bool accept = false;
00554 
00555    // if the new point is more probable, keep it ...
00556    if (p1 >= p0)
00557       accept = true;
00558 
00559    // ... or else throw dice.
00560    else
00561    {
00562       double r = log(fMCMCRandom -> Rndm());
00563 
00564       if(r < p1 - p0)
00565          accept = true;
00566    }
00567 
00568    // fill the new point
00569    if(accept)
00570    {
00571       // increase counter
00572       for (int i = 0; i < fMCMCNParameters; ++i)
00573          fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00574 
00575       // copy the point
00576       for(int i = 0; i < fMCMCNParameters; ++i)
00577       {
00578          // save the point
00579          fMCMCx[index + i] = fMCMCxLocal[i];
00580 
00581          // save the probability of the point
00582          fMCMCLogProbx[chain] = p1;
00583       }
00584    }
00585    else
00586    {
00587       // increase counter
00588       for (int i = 0; i < fMCMCNParameters; ++i)
00589          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00590    }
00591 
00592    return accept;
00593 }

bool BCEngineMCMC::MCMCGetNewPointMetropolisHastings ( int chain = 0 )

Definition at line 597 of file BCEngineMCMC.cxx.

00598 {
00599    // calculate index
00600    int index = chain * fMCMCNParameters;
00601 
00602    // save old point
00603    std::vector <double> xold;
00604    for (int i = 0; i < fMCMCNParameters; ++i)
00605       xold.push_back(fMCMCxLocal.at(i));
00606 
00607    // increase counter
00608    fMCMCNIterations[chain]++;
00609 
00610    // get proposal point
00611    if (!this -> MCMCGetProposalPointMetropolisHastings(chain, fMCMCxLocal, xold))
00612       return false;
00613 
00614    // calculate probabilities of the old and new points
00615    double p0 = fMCMCLogProbx[chain] + log(this -> MCMCTrialFunctionIndependent(xold, fMCMCxLocal, false));
00616    double p1 = this -> LogEval(fMCMCxLocal) + log(this -> MCMCTrialFunctionIndependent(xold, fMCMCxLocal, false));
00617 
00618    // flag for accept
00619    bool accept = false;
00620 
00621    // if the new point is more probable, keep it ...
00622    if (p1 >= p0)
00623       accept = true;
00624    // ... or else throw dice.
00625    else
00626    {
00627       if( log(fMCMCRandom -> Rndm()) < p1 - p0)
00628          accept = true;
00629    }
00630 
00631    // fill the new point
00632    if(accept)
00633    {
00634       // increase counter
00635       for (int i = 0; i < fMCMCNParameters; ++i)
00636          fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00637 
00638       // copy the point
00639       for(int i = 0; i < fMCMCNParameters; ++i)
00640       {
00641          // save the point
00642          fMCMCx[index + i] = fMCMCxLocal[i];
00643 
00644          // save the probability of the point
00645          fMCMCLogProbx[chain] = p1;
00646       }
00647    }
00648    else
00649    {
00650       // increase counter
00651       for (int i = 0; i < fMCMCNParameters; ++i)
00652          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00653    }
00654 
00655    return accept;
00656 }

void BCEngineMCMC::MCMCGetNewPointPCA ( )

Definition at line 469 of file BCEngineMCMC.cxx.

00470 {
00471    // get random point in allowed parameter space
00472    for (int i = 0; i < fMCMCNParameters; ++i)
00473       fMCMCx[i] = fMCMCBoundaryMin.at(i) + (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i)) * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00474 }

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

00394 {
00395    // get unscaled random point in the dimension of the chosen
00396    // parameter. this point might not be in the correct volume.
00397    this -> MCMCTrialFunctionSingle(chain, parameter, x);
00398 
00399    // shift the point to the old point (x0) and scale it.
00400    if (pca == false)
00401    {
00402       // copy the old point into the new
00403       for (int i = 0; i < fMCMCNParameters; ++i)
00404          if (i != parameter)
00405             x[i] = fMCMCx[chain * fMCMCNParameters + i];
00406 
00407       // modify the parameter under study
00408       x[parameter] = fMCMCx[chain * fMCMCNParameters + parameter] + x[parameter] * (fMCMCBoundaryMax.at(parameter) - fMCMCBoundaryMin.at(parameter));
00409    }
00410    else
00411    {
00412       // create temporary points in x and p space
00413       double * newp = new double[fMCMCNParameters];
00414       double * newx = new double[fMCMCNParameters];
00415 
00416       for (int i = 0; i < fMCMCNParameters; i++)
00417       {
00418          newp[i] = 0.;
00419          newx[i] = 0.;
00420       }
00421 
00422       // get the old point in x space
00423       for (int i = 0; i < fMCMCNParameters; i++)
00424          newx[i] = fMCMCx[chain * fMCMCNParameters + i];
00425 
00426       // transform the old point into p space
00427       fMCMCPCA -> X2P(newx, newp);
00428 
00429       // add new trial point to old point
00430       newp[parameter] += x[parameter] * sqrt(fMCMCPCAVariance[parameter]);
00431 
00432       // transform new point back to x space
00433 //    fMCMCPCA -> P2X(newp, newx, fMCMCNParameters);
00434       fMCMCPCA -> P2X(newp, newx, fMCMCNDimensionsPCA);
00435 
00436       // copy point into vector
00437       for (int i = 0; i < fMCMCNParameters; ++i)
00438          x[i] = newx[i];
00439 
00440       delete [] newp;
00441       delete [] newx;
00442    }
00443 
00444    // check if the point is in the correct volume.
00445    for (int i = 0; i < fMCMCNParameters; ++i)
00446       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00447          return false;
00448 
00449    return true;
00450 }

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

Definition at line 334 of file BCEngineMCMC.cxx.

00335 {
00336    // get unscaled random point. this point might not be in the correct volume.
00337    this -> MCMCTrialFunction(x);
00338 
00339    // shift the point to the old point (x0) and scale it.
00340    if (pca == false)
00341    {
00342       // get a proposal point from the trial function and scale it
00343       for (int i = 0; i < fMCMCNParameters; ++i)
00344          x[i] = fMCMCx[chain * fMCMCNParameters + i] + x[i] * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));
00345    }
00346    else
00347    {
00348       // create temporary points in x and p space
00349       double * newp = new double[fMCMCNParameters];
00350       double * newx = new double[fMCMCNParameters];
00351 
00352       for (int i = 0; i < fMCMCNParameters; i++)
00353       {
00354          newp[i] = 0.;
00355          newx[i] = 0.;
00356       }
00357 
00358       // get a new trial point
00359       this -> MCMCTrialFunctionAuto(x);
00360 
00361       // get the old point in x space
00362       for (int i = 0; i < fMCMCNParameters; i++)
00363          newx[i] = fMCMCx[chain * fMCMCNParameters + i];
00364 
00365       // transform the old point into p space
00366       fMCMCPCA -> X2P(newx, newp);
00367 
00368       // add new trial point to old point
00369       for (int i = 0; i < fMCMCNParameters; i++)
00370          newp[i] += fMCMCTrialFunctionScale * x[i];
00371 
00372       // transform new point back to x space
00373 //    fMCMCPCA -> P2X(newp, newx, fMCMCNParameters);
00374       fMCMCPCA -> P2X(newp, newx, fMCMCNDimensionsPCA);
00375 
00376       // copy point into vector
00377       for (int i = 0; i < fMCMCNParameters; ++i)
00378          x[i] = newx[i];
00379 
00380       delete [] newp;
00381       delete [] newx;
00382    }
00383 
00384    // check if the point is in the correct volume.
00385    for (int i = 0; i < fMCMCNParameters; ++i)
00386       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00387          return false;
00388 
00389    return true;
00390 }

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

Definition at line 454 of file BCEngineMCMC.cxx.

00455 {
00456    // get a scaled random point.
00457    this -> MCMCTrialFunctionIndependent(xnew, xold, true);
00458 
00459    // check if the point is in the correct volume.
00460    for (int i = 0; i < fMCMCNParameters; ++i)
00461       if ((xnew[i] < fMCMCBoundaryMin[i]) || (xnew[i] > fMCMCBoundaryMax[i]))
00462          return false;
00463 
00464    return true;
00465 }

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

00274 {
00275    // check if ichain is in range
00276    if (ichain < 0 || ichain >= fMCMCNChains)
00277       return 0;
00278 
00279    // check if ipar is in range
00280    if (ipar < 0 || ipar >= fMCMCNParameters)
00281       return 0;
00282 
00283    // return component of ipar point in the ichain chain
00284    return fMCMCTrialFunctionScaleFactor.at(ichain *  fMCMCNChains + ipar);
00285 }

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

Definition at line 255 of file BCEngineMCMC.cxx.

00256 {
00257    // create a new vector with the length of fMCMCNParameters
00258    std::vector <double> x;
00259 
00260    // check if ichain is in range
00261    if (ichain < 0 || ichain >= fMCMCNChains)
00262       return x;
00263 
00264    // copy the scale factors into the temporary vector
00265    for (int j = 0; j < fMCMCNParameters; ++j)
00266       x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNParameters + j));
00267 
00268    return x;
00269 }

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

00308 {
00309    // check if ichain is in range
00310    if (ichain < 0 || ichain >= fMCMCNChains)
00311       return 0;
00312 
00313    // check if ipar is in range
00314    if (ipar < 0 || ipar >= fMCMCNParameters)
00315       return 0;
00316 
00317    // return component of jth point in the ith chain
00318    return fMCMCx.at(ichain *  fMCMCNParameters + ipar);
00319 }

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

Definition at line 289 of file BCEngineMCMC.cxx.

00290 {
00291    // create a new vector with the length of fMCMCNParameters
00292    std::vector <double> x;
00293 
00294    // check if ichain is in range
00295    if (ichain < 0 || ichain >= fMCMCNChains)
00296       return x;
00297 
00298    // copy the point in the ichain chain into the temporary vector
00299    for (int j = 0; j < fMCMCNParameters; ++j)
00300       x.push_back(fMCMCx.at(ichain * fMCMCNParameters + j));
00301 
00302    return x;
00303 }

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

Definition at line 179 of file BCEngineMCMC.h.

00180          { return fMCMCx; };

int BCEngineMCMC::MCMCInitialize ( )

Definition at line 1660 of file BCEngineMCMC.cxx.

01661 {
01662    // reset variables
01663    fMCMCNIterations.clear();
01664    fMCMCNIterationsConvergenceLocal.clear();
01665    fMCMCNTrialsTrue.clear();
01666    fMCMCNTrialsFalse.clear();
01667    fMCMCTrialFunctionScaleFactor.clear();
01668    fMCMCMean.clear();
01669    fMCMCVariance.clear();
01670    fMCMCMeanx.clear();
01671    fMCMCVariancex.clear();
01672    fMCMCx.clear();
01673    fMCMCLogProbx.clear();
01674    fMCMCMaximumPoints.clear();
01675    fMCMCMaximumLogProb.clear();
01676 // fMCMCRelativePrecisionModeValues.clear();
01677    fMCMCNumericalPrecisionModeValues.clear();
01678    fMCMCNIterationsConvergenceGlobal = -1;
01679    fMCMCFlagConvergenceGlobal = false;
01680    fMCMCRValueParameters.clear();
01681 
01682    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01683       if (fMCMCH1Marginalized[i])
01684          delete fMCMCH1Marginalized[i];
01685 
01686    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01687       if (fMCMCH2Marginalized[i])
01688          delete fMCMCH2Marginalized[i];
01689 
01690    // clear plots
01691    fMCMCH1Marginalized.clear();
01692    fMCMCH2Marginalized.clear();
01693 
01694 // if (fMCMCH1RValue)
01695 //    delete fMCMCH1RValue;
01696 
01697 // if (fMCMCH1Efficiency)
01698 //    delete fMCMCH1Efficiency;
01699 
01700 // free memory for vectors
01701    fMCMCNIterationsConvergenceLocal.assign(fMCMCNChains, -1);
01702    fMCMCNIterations.assign(fMCMCNChains, 0);
01703    fMCMCMean.assign(fMCMCNChains, 0);
01704    fMCMCVariance.assign(fMCMCNChains, 0);
01705    fMCMCLogProbx.assign(fMCMCNChains, -1.0);
01706    fMCMCMaximumLogProb.assign(fMCMCNChains, -1.0);
01707 
01708    fMCMCNumericalPrecisionModeValues.assign(fMCMCNParameters, 0.0);
01709 
01710    fMCMCNTrialsTrue.assign(fMCMCNChains * fMCMCNParameters, 0);
01711    fMCMCNTrialsFalse.assign(fMCMCNChains * fMCMCNParameters, 0);
01712    fMCMCMaximumPoints.assign(fMCMCNChains * fMCMCNParameters, 0.);
01713    fMCMCMeanx.assign(fMCMCNChains * fMCMCNParameters, 0);
01714    fMCMCVariancex.assign(fMCMCNChains * fMCMCNParameters, 0);
01715 
01716    fMCMCRValueParameters.assign(fMCMCNParameters, 0.);
01717 
01718    if (fMCMCTrialFunctionScaleFactorStart.size() == 0)
01719       fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fMCMCNParameters, 1.0);
01720    else
01721       for (int i = 0; i < fMCMCNChains; ++i)
01722          for (int j = 0; j < fMCMCNParameters; ++j)
01723             fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j));
01724 
01725    // set initial position
01726    if (fMCMCFlagInitialPosition == 2) // user defined points
01727    {
01728       // define flag
01729       bool flag = true;
01730 
01731       // check the length of the array of initial positions
01732       if (int(fMCMCInitialPosition.size()) != (fMCMCNChains * fMCMCNParameters))
01733       {
01734          BCLog::Out(BCLog::error, BCLog::error, "BCEngine::MCMCInitialize : Length of vector containing initial positions does not have required length.");
01735          flag = false;
01736       }
01737 
01738       // check the boundaries
01739       if (flag)
01740       {
01741          for (int j = 0; j < fMCMCNChains; ++j)
01742             for (int i = 0; i < fMCMCNParameters; ++i)
01743             {
01744                if (fMCMCInitialPosition[j * fMCMCNParameters + i] < fMCMCBoundaryMin[i] ||
01745                      fMCMCInitialPosition[j * fMCMCNParameters + i] > fMCMCBoundaryMax[i])
01746                {
01747                   BCLog::OutError("BCEngine::MCMCInitialize : Initial position out of boundaries.");
01748                   flag = false;
01749                }
01750             }
01751       }
01752 
01753       // check flag
01754       if (!flag)
01755          fMCMCFlagInitialPosition = 1;
01756    }
01757 
01758    if (fMCMCFlagInitialPosition == 0) // center of the interval
01759       for (int j = 0; j < fMCMCNChains; ++j)
01760          for (int i = 0; i < fMCMCNParameters; ++i)
01761             fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01762 
01763    else if (fMCMCFlagInitialPosition == 2) // user defined
01764    {
01765       for (int j = 0; j < fMCMCNChains; ++j)
01766          for (int i = 0; i < fMCMCNParameters; ++i)
01767             fMCMCx.push_back(fMCMCInitialPosition.at(j * fMCMCNParameters + i));
01768    }
01769 
01770    else
01771       for (int j = 0; j < fMCMCNChains; ++j) // random number (default)
01772          for (int i = 0; i < fMCMCNParameters; ++i)
01773             fMCMCx.push_back(fMCMCBoundaryMin[i] + fMCMCRandom -> Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01774 
01775    // copy the point of the first chain
01776    fMCMCxLocal.clear();
01777    for (int i = 0; i < fMCMCNParameters; ++i)
01778       fMCMCxLocal.push_back(fMCMCx[i]);
01779 
01780    // define 1-dimensional histograms for marginalization
01781    for(int i = 0; i < fMCMCNParameters; ++i)
01782    {
01783       double hmin1 = fMCMCBoundaryMin.at(i);
01784       double hmax1 = fMCMCBoundaryMax.at(i);
01785 
01786       TH1D * h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i), "", fMCMCH1NBins[i], hmin1, hmax1);
01787       fMCMCH1Marginalized.push_back(h1);
01788    }
01789 
01790    for(int i = 0; i < fMCMCNParameters; ++i)
01791       for (int k = 0; k < i; ++k)
01792       {
01793          double hmin1 = fMCMCBoundaryMin.at(k);
01794          double hmax1 = fMCMCBoundaryMax.at(k);
01795          double hmin2 = fMCMCBoundaryMin.at(i);
01796          double hmax2 = fMCMCBoundaryMax.at(i);
01797 
01798          TH2D * h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, k), "",
01799                fMCMCH1NBins[i], hmin1, hmax1,
01800                fMCMCH1NBins[k], hmin2, hmax2);
01801          fMCMCH2Marginalized.push_back(h2);
01802       }
01803 
01804    // define plot for R-value
01805 // if (fMCMCNChains > 1)
01806 // {
01807 //    fMCMCH1RValue = new TH1D("h1_rvalue", ";Iteration;R-value",
01808 //           100, 0, double(fMCMCNIterationsMax));
01809 //    fMCMCH1RValue -> SetStats(false);
01810 // }
01811 
01812 // fMCMCH1Efficiency = new TH1D("h1_efficiency", ";Chain;Efficiency",
01813 //       fMCMCNParameters, -0.5, fMCMCNParameters - 0.5);
01814 // fMCMCH1Efficiency -> SetStats(false);
01815 
01816    return 1;
01817 }

void BCEngineMCMC::MCMCInitializeMarkovChains ( )

Definition at line 1642 of file BCEngineMCMC.cxx.

01643 {
01644    // evaluate function at the starting point
01645    std::vector <double> x0;
01646 
01647    for (int j = 0; j < fMCMCNChains; ++j)
01648    {
01649       x0.clear();
01650       for (int i = 0; i < fMCMCNParameters; ++i)
01651          x0.push_back(fMCMCx[j * fMCMCNParameters + i]);
01652       fMCMCLogProbx[j] = this -> LogEval(x0);
01653    }
01654 
01655    x0.clear();
01656 }

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

Reimplemented in BCIntegrate.

Definition at line 549 of file BCEngineMCMC.h.

00550          {};

int BCEngineMCMC::MCMCMetropolis ( )

Definition at line 1422 of file BCEngineMCMC.cxx.

01423 {
01424    // check if prerun has been performed
01425    if (!fMCMCFlagPreRun)
01426       this -> MCMCMetropolisPreRun();
01427 
01428    // print to screen
01429    BCLog::Out(BCLog::summary, BCLog::summary, "Run Metropolis MCMC...");
01430 
01431    if (fMCMCFlagPCA)
01432       BCLog::Out(BCLog::detail, BCLog::detail, " --> PCA switched on");
01433 
01434    // reset run statistics
01435    this -> MCMCResetRunStatistics();
01436 
01437    // perform run
01438    BCLog::Out(BCLog::summary, BCLog::summary,
01439          Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));
01440 
01441 // int counterupdate = 0;
01442 // bool convergence = false;
01443 // bool flagefficiency = false;
01444 
01445 // std::vector <double> efficiency;
01446 
01447 // for (int i = 0; i < fMCMCNParameters; ++i)
01448 //    for (int j = 0; j < fMCMCNChains; ++j)
01449 //       efficiency.push_back(0.0);
01450 
01451    int nwrite = fMCMCNIterationsRun/10;
01452    if(nwrite < 100)
01453       nwrite=100;
01454    else if(nwrite < 500)
01455       nwrite=1000;
01456    else if(nwrite < 10000)
01457       nwrite=1000;
01458    else
01459       nwrite=10000;
01460 
01461    // start the run
01462    for (int iiterations = 0; iiterations < fMCMCNIterationsRun; ++iiterations)
01463    {
01464       if ( (iiterations+1)%nwrite == 0 )
01465          BCLog::Out(BCLog::detail, BCLog::detail,
01466                Form(" --> iteration number %i (%.2f%%)", iiterations+1, (double)(iiterations+1)/(double)fMCMCNIterationsRun*100.));
01467 
01468       // if the flag is set then run over the parameters one after the other.
01469       if (fMCMCFlagOrderParameters)
01470       {
01471          // loop over parameters
01472          for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
01473             {
01474                // loop over chains
01475                for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01476                   this -> MCMCGetNewPointMetropolis(ichains, iparameters, fMCMCFlagPCA);
01477 
01478                // update statistics
01479                this -> MCMCUpdateStatistics();
01480 
01481                // do anything interface
01482                this -> MCMCIterationInterface();
01483             } // end loop over all parameters
01484       }
01485       // if the flag is not set then run over the parameters at the same time.
01486       else
01487       {
01488          // loop over chains
01489          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01490             // get new point
01491             this -> MCMCGetNewPointMetropolis(ichains, false);
01492 
01493          // update statistics
01494          this -> MCMCUpdateStatistics();
01495 
01496          // do anything interface
01497          this -> MCMCIterationInterface();
01498       }
01499    } // end run
01500 
01501    // print convergence status
01502    BCLog::Out(BCLog::summary, BCLog::summary,
01503          Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));
01504 
01505    // print modes
01506 
01507    // find global maximum
01508    double probmax = fMCMCMaximumLogProb.at(0);
01509    int probmaxindex = 0;
01510 
01511    // loop over all chains and find the maximum point
01512    for (int i = 1; i < fMCMCNChains; ++i)
01513    if (fMCMCMaximumLogProb.at(i) > probmax)
01514    {
01515       probmax = fMCMCMaximumLogProb.at(i);
01516       probmaxindex = i;
01517    }
01518 
01519    BCLog::Out(BCLog::detail, BCLog::detail, " --> Global mode from MCMC:");
01520    int ndigits = (int) log10(fMCMCNParameters);
01521    for (int i = 0; i < fMCMCNParameters; ++i)
01522       BCLog::Out(BCLog::detail, BCLog::detail,
01523             Form( TString::Format(" -->      parameter %%%di:   %%.4g", ndigits+1),
01524                   i, fMCMCMaximumPoints[probmaxindex * fMCMCNParameters + i]));
01525 
01526    // set flag
01527    fMCMCFlagPreRun = false;
01528 
01529    return 1;
01530 }

int BCEngineMCMC::MCMCMetropolisHastings ( )

Definition at line 1534 of file BCEngineMCMC.cxx.

01535 {
01536    BCLog::Out(BCLog::summary, BCLog::summary, "Run Metropolis-Hastings MCMC.");
01537 
01538    // initialize Markov chain
01539    this -> MCMCInitializeMarkovChains();
01540 
01541    // perform burn-in run
01542    BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Start burn-in run with %i iterations.", fMCMCNIterationsBurnIn));
01543    for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
01544       for (int j = 0; j < fMCMCNChains; ++j)
01545          this -> MCMCGetNewPointMetropolisHastings(j);
01546 
01547    // reset run statistics
01548    this -> MCMCResetRunStatistics();
01549 
01550    // perform run
01551    BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Perform MCMC run with %i iterations.", fMCMCNIterationsMax));
01552    for (int i = 0; i < fMCMCNIterationsMax; ++i)
01553    {
01554       for (int j = 0; j < fMCMCNChains; ++j)
01555          // get new point and increase counters
01556          this -> MCMCGetNewPointMetropolisHastings(j);
01557 
01558       // update statistics
01559       this -> MCMCUpdateStatistics();
01560 
01561       // call user interface
01562       this -> MCMCIterationInterface();
01563    }
01564 
01565    if (fMCMCNIterationsConvergenceGlobal > 0)
01566       BCLog::Out(BCLog::detail, BCLog::detail,
01567             Form(" --> Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
01568    else
01569       BCLog::Out(BCLog::detail, BCLog::detail,
01570             Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));
01571 
01572    // debug
01573    if (DEBUG)
01574    {
01575       for (int i = 0; i < fMCMCNChains; ++i)
01576       {
01577          std::cout << i << " " << fMCMCMaximumLogProb[i] << std::endl;
01578          for (int j = 0; j < fMCMCNParameters; ++j)
01579             std::cout << fMCMCMaximumPoints[i * fMCMCNParameters + j] << " ";
01580          std::cout << std::endl;
01581       }
01582    }
01583 
01584    return 1;
01585 }

int BCEngineMCMC::MCMCMetropolisPreRun ( )

Definition at line 890 of file BCEngineMCMC.cxx.

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

void BCEngineMCMC::MCMCResetRunStatistics ( )

Definition at line 1589 of file BCEngineMCMC.cxx.

01590 {
01591    // debugKK -> this is set explicitly
01592    // fMCMCNIterationsConvergenceGlobal = -1;
01593 
01594    for (int j = 0; j < fMCMCNChains; ++j)
01595    {
01596       fMCMCNIterations[j]  = 0;
01597       fMCMCNTrialsTrue[j]  = 0;
01598       fMCMCNTrialsFalse[j] = 0;
01599       fMCMCMean[j]         = 0;
01600       fMCMCVariance[j]     = 0;
01601 
01602       for (int k = 0; k < fMCMCNParameters; ++k)
01603       {
01604          fMCMCNTrialsTrue[j * fMCMCNParameters + k]  = 0;
01605          fMCMCNTrialsFalse[j * fMCMCNParameters + k] = 0;
01606       }
01607    }
01608 
01609    // reset marginalized distributions
01610    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01611       fMCMCH1Marginalized[i] -> Reset();
01612 
01613    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01614       fMCMCH2Marginalized[i] -> Reset();
01615 
01616 // if (fMCMCH1RValue)
01617 //   fMCMCH1RValue -> Reset();
01618 
01619 // if (fMCMCH1Efficiency)
01620 //   fMCMCH1Efficiency -> Reset();
01621 
01622    fMCMCRValue = 100;
01623 }

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

Definition at line 371 of file BCEngineMCMC.h.

00372          { fMCMCFlagInitialPosition = flag; };

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

Definition at line 377 of file BCEngineMCMC.h.

00378       { fMCMCFlagOrderParameters = flag; };

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

Definition at line 397 of file BCEngineMCMC.h.

00398          { fMCMCFlagPCA = flag; };

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

Definition at line 407 of file BCEngineMCMC.h.

00408          { fMCMCFlagPCATruncate = flag; };

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

Definition at line 180 of file BCEngineMCMC.cxx.

00181 {
00182    // create new vector 
00183    std::vector <double> y0s; 
00184    
00185    // loop over vector elements
00186    for (int i = 0; i < int(x0s.size()); ++i)
00187       for (int j = 0; j < int((x0s.at(i)).size()); ++j)
00188          y0s.push_back((x0s.at(i)).at(j)); 
00189 
00190    this -> MCMCSetInitialPositions(y0s); 
00191 }

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

Definition at line 164 of file BCEngineMCMC.cxx.

00165 {
00166    // clear the existing initial position vector
00167    fMCMCInitialPosition.clear(); 
00168 
00169    // copy the initial positions
00170    int n = int(x0s.size()); 
00171 
00172    for (int i = 0; i < n; ++i)
00173       fMCMCInitialPosition.push_back(x0s.at(i)); 
00174 
00175    // use these intial positions for the Markov chain
00176    this -> MCMCSetFlagInitialPosition(2);
00177 }

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

Definition at line 344 of file BCEngineMCMC.h.

00345          { fMCMCFlagIterationsAuto = flag; };

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

Definition at line 194 of file BCEngineMCMC.cxx.

00195 {
00196 // clear vector
00197    fMCMCTrees.clear();
00198 
00199    // copy tree
00200    for (int i = 0; i < int(trees.size()); ++i)
00201       fMCMCTrees.push_back(trees[i]);
00202 }

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

Definition at line 354 of file BCEngineMCMC.h.

00355          { fMCMCEfficiencyMax = efficiency; };

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

Definition at line 349 of file BCEngineMCMC.h.

00350          { 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 335 of file BCEngineMCMC.h.

00336          { fMCMCNIterationsPCA = n; };

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

Definition at line 330 of file BCEngineMCMC.h.

00331          { fMCMCNIterationsPreRunMin = 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 338 of file BCEngineMCMC.h.

00339       { fMCMCNIterationsUpdate = n; };

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

Definition at line 413 of file BCEngineMCMC.h.

00414          { fMCMCPCAMinimumRatio = ratio; };

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

Definition at line 382 of file BCEngineMCMC.h.

00383          { fMCMCRValueCriterion = r; };

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

Definition at line 387 of file BCEngineMCMC.h.

00388          { 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 1869 of file BCEngineMCMC.cxx.

01870 {
01871    this -> MCMCSetValuesDetail();
01872 }

void BCEngineMCMC::MCMCSetValuesDetail ( )

Set the values for a detailed MCMC run.

Definition at line 1902 of file BCEngineMCMC.cxx.

01903 {
01904    fMCMCNChains              = 5;
01905    fMCMCNIterationsMax       = 1000000;
01906    fMCMCNIterationsRun       = 100000;
01907    fMCMCNIterationsBurnIn    = 0;
01908    fMCMCNIterationsPCA       = 10000;
01909    fMCMCNIterationsPreRunMin = 500;
01910    fMCMCFlagIterationsAuto   = false;
01911    fMCMCTrialFunctionScale   = 1.0;
01912    fMCMCFlagInitialPosition  = 1;
01913    fMCMCRValueCriterion      = 0.1;
01914    fMCMCRValueParametersCriterion = 0.1;
01915    fMCMCNIterationsConvergenceGlobal = -1;
01916    fMCMCFlagConvergenceGlobal = false;
01917    fMCMCRValue               = 100;
01918    fMCMCFlagPCA              = false;
01919    fMCMCFlagPCATruncate      = false;
01920    fMCMCPCA                  = 0;
01921    fMCMCPCAMinimumRatio      = 1e-7;
01922    fMCMCNIterationsUpdate    = 1000;
01923    fMCMCFlagOrderParameters  = true;
01924 }

void BCEngineMCMC::MCMCSetValuesQuick ( )

Set the values for a quick MCMC run.

Definition at line 1876 of file BCEngineMCMC.cxx.

01877 {
01878    fMCMCNChains              = 1;
01879    fMCMCNIterationsMax       = 1000;
01880    fMCMCNIterationsRun       = 10000;
01881    fMCMCNIterationsBurnIn    = 0;
01882    fMCMCNIterationsPCA       = 0;
01883    fMCMCNIterationsPreRunMin = 0;
01884    fMCMCFlagIterationsAuto   = false;
01885    fMCMCTrialFunctionScale   = 1.0;
01886    fMCMCFlagInitialPosition  = 1;
01887    fMCMCRValueCriterion      = 0.1;
01888    fMCMCRValueParametersCriterion = 0.1;
01889    fMCMCNIterationsConvergenceGlobal = -1;
01890    fMCMCFlagConvergenceGlobal = false;
01891    fMCMCRValue               = 100;
01892    fMCMCFlagPCA              = false;
01893    fMCMCFlagPCATruncate      = false;
01894    fMCMCPCA                  = 0;
01895    fMCMCPCAMinimumRatio      = 1e-7;
01896    fMCMCNIterationsUpdate    = 1000;
01897    fMCMCFlagOrderParameters  = true;
01898 }

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

Definition at line 359 of file BCEngineMCMC.h.

00360          { fMCMCFlagWriteChainToFile = flag; };

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

Definition at line 211 of file BCEngineMCMC.cxx.

00212 {
00213    // use uniform distribution for now
00214    for (int i = 0; i < fMCMCNParameters; ++i)
00215       x[i] = fMCMCTrialFunctionScale * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00216 }

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

Definition at line 246 of file BCEngineMCMC.cxx.

00247 {
00248    // use uniform distribution for now
00249    for (int i = 0; i < fMCMCNParameters; ++i)
00250       x[i] = fMCMCRandom -> Gaus(fMCMCPCAMean[i], sqrt(fMCMCPCAVariance[i]));
00251 }

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

Definition at line 230 of file BCEngineMCMC.cxx.

00231 {
00232    // use uniform distribution for now
00233    if (newpoint)
00234       for (int i = 0; i < fMCMCNParameters; ++i)
00235          xnew[i] = fMCMCRandom -> Rndm() * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));
00236 
00237    double prob = 1.0;
00238    for (int i = 0; i < fMCMCNParameters; ++i)
00239       prob /= fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i);
00240 
00241    return prob;
00242 }

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

Definition at line 220 of file BCEngineMCMC.cxx.

00221 {
00222    // no check of range for performance reasons
00223 
00224    // use uniform distribution
00225    x[iparameter] = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00226 }

void BCEngineMCMC::MCMCUpdateStatistics ( )

Definition at line 859 of file BCEngineMCMC.cxx.

00860 {
00861    // check if maximum is reached
00862    this -> MCMCUpdateStatisticsCheckMaximum();
00863 
00864    // fill histograms of marginalized distributions
00865    this -> MCMCUpdateStatisticsFillHistograms();
00866 
00867    // test convergence of single Markov chains
00868    //  --> not implemented yet
00869 
00870    // test convergence of all Markov chains
00871    // debugKK -> this will be called explicitly
00872    // this -> MCMCUpdateStatisticsTestConvergenceAllChains();
00873 
00874    // fill Markov chains
00875    if (fMCMCFlagWriteChainToFile)
00876       this -> MCMCUpdateStatisticsWriteChains();
00877 }

void BCEngineMCMC::MCMCUpdateStatisticsCheckMaximum ( )

Definition at line 698 of file BCEngineMCMC.cxx.

00699 {
00700    // loop over all chains
00701    for (int i = 0; i < fMCMCNChains; ++i)
00702    {
00703       if (fMCMCLogProbx[i] > fMCMCMaximumLogProb[i] || fMCMCNIterations[i] == 1)
00704       {
00705          fMCMCMaximumLogProb[i] = fMCMCLogProbx[i];
00706          for (int j = 0; j < fMCMCNParameters; ++j)
00707             fMCMCMaximumPoints[i * fMCMCNParameters + j] = fMCMCx[i * fMCMCNParameters + j];
00708       }
00709    }
00710 }

void BCEngineMCMC::MCMCUpdateStatisticsFillHistograms ( )

Definition at line 714 of file BCEngineMCMC.cxx.

00715 {
00716    // loop over chains
00717    for (int i = 0; i < fMCMCNChains; ++i)
00718    {
00719       // fill each 1-dimensional histogram
00720       for (int j = 0; j < fMCMCNParameters; ++j)
00721          fMCMCH1Marginalized[j] -> Fill(fMCMCx[i * fMCMCNParameters + j]);
00722 
00723       // fill each 2-dimensional histogram
00724       int counter = 0;
00725 
00726       for (int j = 0; j < fMCMCNParameters; ++j)
00727          for (int k = 0; k < j; ++k)
00728          {
00729             fMCMCH2Marginalized[counter] -> Fill(
00730                   fMCMCx[i * fMCMCNParameters + k],
00731                   fMCMCx[i * fMCMCNParameters + j]);
00732             counter ++;
00733          }
00734    }
00735 }

void BCEngineMCMC::MCMCUpdateStatisticsModeConvergence ( )

Definition at line 660 of file BCEngineMCMC.cxx.

00661 {
00662    double * mode_minimum = new double[fMCMCNParameters];
00663    double * mode_maximum = new double[fMCMCNParameters];
00664    double * mode_average = new double[fMCMCNParameters];
00665 
00666    // set initial values
00667    for (int j = 0; j < fMCMCNParameters; ++j)
00668    {
00669       mode_minimum[j] = fMCMCMaximumPoints[j];
00670       mode_maximum[j] = fMCMCMaximumPoints[j];
00671       mode_average[j] = 0;
00672    }
00673 
00674    // calculate the maximum difference in each dimension
00675    for (int i = 0; i < fMCMCNChains; ++i)
00676       for (int j = 0; j < fMCMCNParameters; ++j)
00677       {
00678          if (fMCMCMaximumPoints[i * fMCMCNParameters + j] < mode_minimum[j])
00679             mode_minimum[j] = fMCMCMaximumPoints[i * fMCMCNParameters + j];
00680 
00681          if (fMCMCMaximumPoints[i * fMCMCNParameters + j] > mode_maximum[j])
00682             mode_maximum[j] = fMCMCMaximumPoints[i * fMCMCNParameters + j];
00683 
00684          mode_average[j] += fMCMCMaximumPoints[i * fMCMCNParameters + j] / double(fMCMCNChains);
00685       }
00686 
00687    for (int j = 0; j < fMCMCNParameters; ++j)
00688       fMCMCNumericalPrecisionModeValues[j] = (mode_maximum[j] - mode_minimum[j]);
00689 //    fMCMCRelativePrecisionModeValues[j] = (mode_maximum[j] - mode_minimum[j]) / mode_average[j];
00690 
00691    delete [] mode_minimum;
00692    delete [] mode_maximum;
00693    delete [] mode_average;
00694 }

void BCEngineMCMC::MCMCUpdateStatisticsTestConvergenceAllChains ( )

Definition at line 739 of file BCEngineMCMC.cxx.

00740 {
00741    // calculate number of entries in this part of the chain
00742    int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0];
00743 
00744    // do not evaluate the convergence criterion if the numbers of
00745    // elements is too small.
00746    //  if (npoints < fMCMCNIterationsPreRunMin)
00747    //    {
00748    //      fMCMCNIterationsConvergenceGlobal = -1;
00749    //      return;
00750    //    }
00751 
00752    if (fMCMCNChains > 1 && npoints > 1)
00753    {
00754       // define flag for convergence
00755       bool flag_convergence = true;
00756 
00757       // loop over parameters
00758       for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00759             {
00760                double sum = 0; 
00761                double sum2 = 0; 
00762                double sumv = 0; 
00763      
00764                // loop over chains              
00765                for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00766                   {
00767                      // get parameter index
00768                      int index = ichains * fMCMCNParameters + iparameters; 
00769                      
00770                      // calculate mean value of each parameter in the chain for this part                      
00771                      fMCMCMeanx[index] += (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints);   
00772                      
00773                      // calculate variance of each chain for this part                    
00774                      fMCMCVariancex[index] = (1.0 - 1/double(npoints)) * fMCMCVariancex[index] 
00775                         + (fMCMCx[index] - fMCMCMeanx[index]) * (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints - 1); 
00776                      
00777                      sum  += fMCMCMeanx[index]; 
00778                      sum2 += fMCMCMeanx[index] * fMCMCMeanx[index]; 
00779                      sumv += fMCMCVariancex[index]; 
00780                   }
00781                
00782                // calculate r-value for each parameter               
00783                double mean = sum / double(fMCMCNChains); 
00784                double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints); 
00785                double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains); 
00786                
00787                double r = 100.0; 
00788                
00789                if (W > 0)
00790                   {
00791                      r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W); 
00792                      
00793                      fMCMCRValueParameters[iparameters] = r; 
00794                   }
00795                // debugKK
00796                //             std::cout << " here : " << W << " " << fMCMCRValueParameters[iparameters] << std::endl; 
00797                
00798                // set flag to false if convergence criterion is not fulfilled for the parameter                
00799                if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
00800                   flag_convergence = false; 
00801             }
00802       
00803       // convergence criterion applied on the log-likelihood         
00804       double sum = 0; 
00805       double sum2 = 0; 
00806       double sumv = 0; 
00807       
00808       // loop over chains       
00809       for (int i = 0; i < fMCMCNChains; ++i)
00810             {
00811                // calculate mean value of each chain for this part               
00812                fMCMCMean[i] += (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints); 
00813                
00814                // calculate variance of each chain for this part              
00815                if (npoints > 1) 
00816                   fMCMCVariance[i] = (1.0 - 1/double(npoints)) * fMCMCVariance[i] 
00817                      + (fMCMCLogProbx[i] - fMCMCMean[i]) * (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints - 1); 
00818                
00819                sum  += fMCMCMean[i]; 
00820                sum2 += fMCMCMean[i] * fMCMCMean[i]; ; 
00821                sumv += fMCMCVariance[i]; 
00822             }
00823       
00824       // calculate r-value for log-likelihood   
00825       double mean = sum / double(fMCMCNChains); 
00826       double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints); 
00827       double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);       
00828       double r = 100.0; 
00829       
00830       if (W > 0)
00831             {
00832                r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W); 
00833                
00834                fMCMCRValue = r; 
00835             }
00836          
00837       // set flag to false if convergence criterion is not fulfilled for the log-likelihood        
00838       if (! ((fMCMCRValue-1.0) < fMCMCRValueCriterion))
00839             flag_convergence = false; 
00840       
00841       // remember number of iterations needed to converge       
00842       if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
00843             fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters; 
00844     }
00845   
00846 }

void BCEngineMCMC::MCMCUpdateStatisticsWriteChains ( )

Definition at line 850 of file BCEngineMCMC.cxx.

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

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

01840 {
01841    int counter = 0;
01842    int index = 0;
01843 
01844    // search for correct combination
01845    for(int i = 0; i < fMCMCNParameters; i++)
01846       for (int j = 0; j < i; ++j)
01847       {
01848          if(j == index1 && i == index2)
01849             index = counter;
01850          counter++;
01851       }
01852 
01853    if((int)fMCMCH2Marginalized.size()<=index)
01854       return 0;
01855 
01856    if(h==0)
01857       return 0;
01858 
01859    if((int)fMCMCH2Marginalized.size()==index)
01860       fMCMCH2Marginalized.push_back(h);
01861    else
01862       fMCMCH2Marginalized[index]=h;
01863 
01864    return index;
01865 }

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

Definition at line 1821 of file BCEngineMCMC.cxx.

01822 {
01823    if((int)fMCMCH1Marginalized.size()<=index)
01824       return 0;
01825 
01826    if(h==0)
01827       return 0;
01828 
01829    if((int)fMCMCH1Marginalized.size()==index)
01830       fMCMCH1Marginalized.push_back(h);
01831    else
01832       fMCMCH1Marginalized[index]=h;
01833 
01834    return index;
01835 }

Member Data Documentation

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

Definition at line 602 of file BCEngineMCMC.h.

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

Definition at line 601 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCEfficiencyMax [protected]

Definition at line 738 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCEfficiencyMin [protected]

Definition at line 734 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagConvergenceGlobal [protected]

Definition at line 628 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCFlagInitialPosition [protected]

Definition at line 744 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagIterationsAuto [protected]

Definition at line 705 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagOrderParameters [protected]

Definition at line 749 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagPCA [protected]

Definition at line 810 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagPCATruncate [protected]

Definition at line 663 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagPreRun [protected]

Definition at line 723 of file BCEngineMCMC.h.

bool BCEngineMCMC::fMCMCFlagWriteChainToFile [protected]

Definition at line 709 of file BCEngineMCMC.h.

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

Definition at line 818 of file BCEngineMCMC.h.

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

Definition at line 814 of file BCEngineMCMC.h.

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

Definition at line 819 of file BCEngineMCMC.h.

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

Definition at line 730 of file BCEngineMCMC.h.

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

Definition at line 765 of file BCEngineMCMC.h.

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

Definition at line 776 of file BCEngineMCMC.h.

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

Definition at line 771 of file BCEngineMCMC.h.

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

Definition at line 686 of file BCEngineMCMC.h.

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

Definition at line 699 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNChains [protected]

Definition at line 606 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNDimensionsPCA [protected]

Definition at line 673 of file BCEngineMCMC.h.

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

Definition at line 611 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsBurnIn [protected]

Definition at line 645 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsConvergenceGlobal [protected]

Definition at line 624 of file BCEngineMCMC.h.

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

Definition at line 620 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsMax [protected]

Definition at line 632 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsPCA [protected]

Definition at line 650 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsPreRunMin [protected]

Definition at line 640 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsRun [protected]

Definition at line 636 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNIterationsUpdate [protected]

Definition at line 615 of file BCEngineMCMC.h.

int BCEngineMCMC::fMCMCNParameters [protected]

Definition at line 597 of file BCEngineMCMC.h.

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

Definition at line 681 of file BCEngineMCMC.h.

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

Definition at line 680 of file BCEngineMCMC.h.

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

Definition at line 798 of file BCEngineMCMC.h.

TPrincipal* BCEngineMCMC::fMCMCPCA [protected]

Definition at line 806 of file BCEngineMCMC.h.

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

Definition at line 654 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCPCAMinimumRatio [protected]

Definition at line 668 of file BCEngineMCMC.h.

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

Definition at line 658 of file BCEngineMCMC.h.

MCMCPointerToGetProposalPoint BCEngineMCMC::fMCMCPointerToGetProposalPoint [private]

Definition at line 591 of file BCEngineMCMC.h.

TRandom3* BCEngineMCMC::fMCMCRandom [protected]

Definition at line 802 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValue [protected]

Definition at line 788 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValueCriterion [protected]

Definition at line 780 of file BCEngineMCMC.h.

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

Definition at line 790 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCRValueParametersCriterion [protected]

Definition at line 784 of file BCEngineMCMC.h.

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

Definition at line 693 of file BCEngineMCMC.h.

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

Definition at line 694 of file BCEngineMCMC.h.

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

Definition at line 829 of file BCEngineMCMC.h.

double BCEngineMCMC::fMCMCTrialFunctionScale [protected]

Definition at line 713 of file BCEngineMCMC.h.

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

Definition at line 718 of file BCEngineMCMC.h.

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

Definition at line 719 of file BCEngineMCMC.h.

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

Definition at line 687 of file BCEngineMCMC.h.

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

Definition at line 700 of file BCEngineMCMC.h.

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

Definition at line 756 of file BCEngineMCMC.h.

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

Definition at line 760 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	MCMCSetInitialPositions (std::vector< std::vector< double > > x0s)
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	MCMCSetNIterationsPreRunMin (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