An engine class for Markov Chain Monte Carlo. More...
#include <BCEngineMCMC.h>
Inherited by BCIntegrate.
An engine class for Markov Chain Monte Carlo.
Definition at line 35 of file BCEngineMCMC.h.
typedef bool(BCEngineMCMC::* BCEngineMCMC::MCMCPointerToGetProposalPoint)(int chain, std::vector< double > xnew, std::vector< double > xold) const [private] |
Defines a type of a pointer to a member function.
Definition at line 579 of file BCEngineMCMC.h.
An enumerator for the status of a test.
Definition at line 44 of file BCEngineMCMC.h.
{ kLow, kMedium, kHigh, kVeryHigh };
BCEngineMCMC::BCEngineMCMC | ( | ) |
Default constructor.
Definition at line 24 of file BCEngineMCMC.cxx.
{ // set default parameters for the mcmc MCMCSetValuesDefault(); // initialize random number generator fRandom = new TRandom3(0); }
BCEngineMCMC::BCEngineMCMC | ( | int | n | ) |
Constructor.
n | number of chains |
Definition at line 34 of file BCEngineMCMC.cxx.
{ // set number of chains to n fMCMCNChains = n; // call default constructor BCEngineMCMC(); }
BCEngineMCMC::BCEngineMCMC | ( | const BCEngineMCMC & | enginemcmc | ) |
Default copy constructor.
Definition at line 183 of file BCEngineMCMC.cxx.
{ fMCMCPointerToGetProposalPoint = enginemcmc.fMCMCPointerToGetProposalPoint; fMCMCNParameters = enginemcmc.fMCMCNParameters; fMCMCBoundaryMin = enginemcmc.fMCMCBoundaryMin; fMCMCBoundaryMax = enginemcmc.fMCMCBoundaryMax; fMCMCFlagsFillHistograms = enginemcmc.fMCMCFlagsFillHistograms; fMCMCNChains = enginemcmc.fMCMCNChains; fMCMCNLag = enginemcmc.fMCMCNLag; fMCMCNIterations = enginemcmc.fMCMCNIterations; fMCMCCurrentIteration = enginemcmc.fMCMCCurrentIteration; fMCMCCurrentChain = enginemcmc.fMCMCCurrentChain; fMCMCNIterationsUpdate = enginemcmc.fMCMCNIterationsUpdate; fMCMCNIterationsUpdateMax = enginemcmc.fMCMCNIterationsUpdateMax; fMCMCNIterationsConvergenceGlobal = enginemcmc.fMCMCNIterationsConvergenceGlobal; fMCMCFlagConvergenceGlobal = enginemcmc.fMCMCFlagConvergenceGlobal; fMCMCNIterationsMax = enginemcmc.fMCMCNIterationsMax; fMCMCNIterationsRun = enginemcmc.fMCMCNIterationsRun; fMCMCNIterationsPreRunMin = enginemcmc.fMCMCNIterationsPreRunMin; fMCMCNTrialsTrue = enginemcmc.fMCMCNTrialsTrue; fMCMCNTrialsFalse = enginemcmc.fMCMCNTrialsFalse; fMCMCFlagWriteChainToFile = enginemcmc.fMCMCFlagWriteChainToFile; fMCMCFlagWritePreRunToFile = enginemcmc.fMCMCFlagWritePreRunToFile; fMCMCTrialFunctionScaleFactor = enginemcmc.fMCMCTrialFunctionScaleFactor; fMCMCTrialFunctionScaleFactorStart = enginemcmc.fMCMCTrialFunctionScaleFactorStart; fMCMCFlagPreRun = enginemcmc.fMCMCFlagPreRun; fMCMCFlagRun = enginemcmc.fMCMCFlagRun; fMCMCInitialPosition = enginemcmc.fMCMCInitialPosition; fMCMCEfficiencyMin = enginemcmc.fMCMCEfficiencyMin; fMCMCEfficiencyMax = enginemcmc.fMCMCEfficiencyMax; fMCMCFlagInitialPosition = enginemcmc.fMCMCFlagInitialPosition; fMCMCFlagOrderParameters = enginemcmc.fMCMCFlagOrderParameters; fMCMCFlagFillHistograms = enginemcmc.fMCMCFlagFillHistograms; fMCMCPhase = enginemcmc.fMCMCPhase; fMCMCCycle = enginemcmc.fMCMCCycle; fMCMCx = enginemcmc.fMCMCx; fMCMCxMax = enginemcmc.fMCMCxMax; fMCMCxMean = enginemcmc.fMCMCxMean; fMCMCxVar = enginemcmc.fMCMCxVar; fMCMCxLocal = enginemcmc.fMCMCxLocal; fMCMCprob = enginemcmc.fMCMCprob; fMCMCprobMax = enginemcmc.fMCMCprobMax; fMCMCprobMean = enginemcmc.fMCMCprobMean; fMCMCprobVar = enginemcmc.fMCMCprobVar; fMCMCRValueUseStrict = enginemcmc.fMCMCRValueUseStrict; fMCMCRValueCriterion = enginemcmc.fMCMCRValueCriterion ; fMCMCRValueParametersCriterion = enginemcmc.fMCMCRValueParametersCriterion; fMCMCRValue = enginemcmc.fMCMCRValue; fMCMCRValueParameters = enginemcmc.fMCMCRValueParameters; if (enginemcmc.fRandom) fRandom = new TRandom3(*enginemcmc.fRandom); else fRandom = 0; fMCMCH1NBins = enginemcmc.fMCMCH1NBins; for (int i = 0; i < int(enginemcmc.fMCMCH1Marginalized.size()); ++i) { if (enginemcmc.fMCMCH1Marginalized.at(i)) fMCMCH1Marginalized.push_back(new TH1D(*(enginemcmc.fMCMCH1Marginalized.at(i)))); else fMCMCH1Marginalized.push_back(0); } for (int i = 0; i < int(enginemcmc.fMCMCH2Marginalized.size()); ++i) { if (enginemcmc.fMCMCH2Marginalized.at(i)) fMCMCH2Marginalized.push_back(new TH2D(*(enginemcmc.fMCMCH2Marginalized.at(i)))); else fMCMCH2Marginalized.push_back(0); } for (int i = 0; i < int(enginemcmc.fMCMCTrees.size()); ++i) { // debugKK fMCMCTrees.push_back(0); } }
BCEngineMCMC::~BCEngineMCMC | ( | ) | [virtual] |
Default destructor.
Definition at line 163 of file BCEngineMCMC.cxx.
{ // delete random number generator if (fRandom) delete fRandom; // delete 1-d marginalized distributions for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i) if (fMCMCH1Marginalized[i]) delete fMCMCH1Marginalized[i]; fMCMCH1Marginalized.clear(); // delete 2-d marginalized distributions for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i) if (fMCMCH2Marginalized[i]) delete fMCMCH2Marginalized[i]; fMCMCH2Marginalized.clear(); }
double BCEngineMCMC::LogEval | ( | const std::vector< double > & | parameters | ) | [virtual] |
Needs to be overloaded in the derived class.
Reimplemented in BCIntegrate, and BCModel.
Definition at line 900 of file BCEngineMCMC.cxx.
{ // test function for now // this will be overloaded by the user return 0.0; }
int BCEngineMCMC::MCMCAddParameter | ( | double | min, | |
double | max | |||
) |
Adds a parameter.
min | minimum value of the parameter | |
max | maximum value of the parameter |
Definition at line 1474 of file BCEngineMCMC.cxx.
{ // add the boundaries to the corresponding vectors fMCMCBoundaryMin.push_back(min); fMCMCBoundaryMax.push_back(max); // set flag for individual parameters fMCMCFlagsFillHistograms.push_back(true); // increase the number of parameters by one fMCMCNParameters++; // return the number of parameters return fMCMCNParameters; }
virtual void BCEngineMCMC::MCMCCurrentPointInterface | ( | std::vector< double > & | point, | |
int | ichain, | |||
bool | accepted | |||
) | [inline, virtual] |
Interface allowing to execute arbitrary code for each new point of the MCMC. This method needs to be overloaded in the derived class
point | point that was generated and checked | |
ichain | index of the chain | |
accepted | flag whether or not the point was accepted for the chain |
Definition at line 564 of file BCEngineMCMC.h.
{ // suppress warnings for unused parameters // with optimization, no code should be generated (void)point; (void)ichain; (void)accepted; }
int BCEngineMCMC::MCMCGetCurrentChain | ( | ) | [inline] |
Definition at line 106 of file BCEngineMCMC.h.
{ return fMCMCCurrentChain; }
int BCEngineMCMC::MCMCGetCurrentIteration | ( | ) | [inline] |
Definition at line 101 of file BCEngineMCMC.h.
{ return fMCMCCurrentIteration; }
int BCEngineMCMC::MCMCGetCycle | ( | ) | [inline] |
Definition at line 216 of file BCEngineMCMC.h.
{ return fMCMCCycle; }
bool BCEngineMCMC::MCMCGetFlagConvergenceGlobal | ( | ) | [inline] |
Definition at line 117 of file BCEngineMCMC.h.
{ return fMCMCFlagConvergenceGlobal; }
int BCEngineMCMC::MCMCGetFlagInitialPosition | ( | ) | [inline] |
Definition at line 236 of file BCEngineMCMC.h.
{ return fMCMCFlagInitialPosition; }
bool BCEngineMCMC::MCMCGetFlagRun | ( | ) | [inline] |
Definition at line 266 of file BCEngineMCMC.h.
{ return fMCMCFlagRun; }
TH1D * BCEngineMCMC::MCMCGetH1Marginalized | ( | int | i | ) |
Retrieve a histogram of the 1D marginalized distribution of a single parameter.
i | index of the parameter |
Definition at line 339 of file BCEngineMCMC.cxx.
{ // check index if (index < 0 || index >= int(fMCMCH1Marginalized.size())) { BCLog::OutWarning("BCEngineMCMC::MCMCGetH1Marginalized. Index out of range."); return 0; } return fMCMCH1Marginalized[index]; }
TH2D * BCEngineMCMC::MCMCGetH2Marginalized | ( | int | i, | |
int | j | |||
) |
Retrieve a histogram of the 2D marginalized distribution for two parameters.
i | index of the first parameter | |
j | index of the second parameter |
Definition at line 352 of file BCEngineMCMC.cxx.
{ int counter = 0; int index = 0; // search for correct combination for(int i = 0; i < fMCMCNParameters; i++) for (int j = 0; j < i; ++j) { if(j == index1 && i == index2) index = counter; counter++; } // check index if (index < 0 || index >= int(fMCMCH2Marginalized.size())) { BCLog::OutWarning("BCEngineMCMC::MCMCGetH2Marginalized. Index out of range."); return 0; } return fMCMCH2Marginalized[index]; }
std::vector<double> BCEngineMCMC::MCMCGetLogProbx | ( | ) | [inline] |
Definition at line 201 of file BCEngineMCMC.h.
{ return fMCMCprob; }
double BCEngineMCMC::MCMCGetLogProbx | ( | int | ichain | ) |
ichain | chain index |
Definition at line 573 of file BCEngineMCMC.cxx.
{ // check if ichain is in range if (ichain < 0 || ichain >= fMCMCNChains) return -1; // return log of the probability at the current point in the ith chain return fMCMCprob.at(ichain); }
TTree* BCEngineMCMC::MCMCGetMarkovChainTree | ( | int | i | ) | [inline] |
Rtrieve the tree containing the Markov chain.
i | index of the Markov chain |
Definition at line 273 of file BCEngineMCMC.h.
{ return fMCMCTrees.at(i); }
std::vector<double> BCEngineMCMC::MCMCGetMaximumLogProb | ( | ) | [inline] |
Definition at line 231 of file BCEngineMCMC.h.
{ return fMCMCprobMax; }
std::vector< double > BCEngineMCMC::MCMCGetMaximumPoint | ( | int | i | ) |
i | The index of the Markov chain |
Definition at line 377 of file BCEngineMCMC.cxx.
{ // create a new vector with the lenght of fMCMCNParameters std::vector<double> x; // check if i is in range if (i < 0 || i >= fMCMCNChains) return x; // copy the point in the ith chain into the temporary vector for (int j = 0; j < fMCMCNParameters; ++j) x.push_back(fMCMCxMax.at(i * fMCMCNParameters + j)); return x; }
std::vector<double> BCEngineMCMC::MCMCGetMaximumPoints | ( | ) | [inline] |
Definition at line 221 of file BCEngineMCMC.h.
{ return fMCMCxMax; }
int BCEngineMCMC::MCMCGetNChains | ( | ) | [inline] |
Definition at line 86 of file BCEngineMCMC.h.
{ return fMCMCNChains; }
bool BCEngineMCMC::MCMCGetNewPointMetropolis | ( | int | chain = 0 |
) |
Generates a new point using the Metropolis algorithm.
chain | chain index |
Definition at line 693 of file BCEngineMCMC.cxx.
{ // calculate index int index = chain * fMCMCNParameters; fMCMCCurrentChain = chain; // increase counter fMCMCNIterations[chain]++; // get proposal point if (!MCMCGetProposalPointMetropolis(chain, fMCMCxLocal)) { // increase counter for (int i = 0; i < fMCMCNParameters; ++i) fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++; // execute user code for every point MCMCCurrentPointInterface(fMCMCxLocal, chain, false); return false; } // calculate probabilities of the old and new points double p0 = fMCMCprob[chain]; double p1 = LogEval(fMCMCxLocal); // flag for accept bool accept = false; // if the new point is more probable, keep it ... if (p1 >= p0) accept = true; // ... or else throw dice. else { double r = log(fRandom->Rndm()); if(r < p1 - p0) accept = true; } // fill the new point if(accept) { // increase counter for (int i = 0; i < fMCMCNParameters; ++i) fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++; // copy the point for(int i = 0; i < fMCMCNParameters; ++i) { // save the point fMCMCx[index + i] = fMCMCxLocal[i]; // save the probability of the point fMCMCprob[chain] = p1; } } else { // increase counter for (int i = 0; i < fMCMCNParameters; ++i) fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++; } // execute user code for every point MCMCCurrentPointInterface(fMCMCxLocal, chain, accept); return accept; }
bool BCEngineMCMC::MCMCGetNewPointMetropolis | ( | int | chain, | |
int | parameter | |||
) |
Definition at line 623 of file BCEngineMCMC.cxx.
{ // calculate index int index = chain * fMCMCNParameters; fMCMCCurrentChain = chain; // increase counter fMCMCNIterations[chain]++; // get proposal point if (!MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal)) { // increase counter fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++; // execute user code for every point MCMCCurrentPointInterface(fMCMCxLocal, chain, false); return false; } // calculate probabilities of the old and new points double p0 = fMCMCprob[chain]; double p1 = LogEval(fMCMCxLocal); // flag for accept bool accept = false; // if the new point is more probable, keep it ... if (p1 >= p0) accept = true; // ... or else throw dice. else { double r = log(fRandom->Rndm()); if(r < p1 - p0) accept = true; } // fill the new point if(accept) { // increase counter fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++; // copy the point for(int i = 0; i < fMCMCNParameters; ++i) { // save the point fMCMCx[index + i] = fMCMCxLocal[i]; // save the probability of the point fMCMCprob[chain] = p1; } } else { // increase counter fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++; } // execute user code for every point MCMCCurrentPointInterface(fMCMCxLocal, chain, accept); return accept; }
std::vector<int> BCEngineMCMC::MCMCGetNIterations | ( | ) | [inline] |
Definition at line 96 of file BCEngineMCMC.h.
{ return fMCMCNIterations; }
int BCEngineMCMC::MCMCGetNIterationsConvergenceGlobal | ( | ) | [inline] |
Definition at line 112 of file BCEngineMCMC.h.
{ return fMCMCNIterationsConvergenceGlobal; }
int BCEngineMCMC::MCMCGetNIterationsMax | ( | ) | [inline] |
Definition at line 122 of file BCEngineMCMC.h.
{ return fMCMCNIterationsMax; }
int BCEngineMCMC::MCMCGetNIterationsPreRunMin | ( | ) | [inline] |
Definition at line 132 of file BCEngineMCMC.h.
{ return fMCMCNIterationsPreRunMin; }
int BCEngineMCMC::MCMCGetNIterationsRun | ( | ) | [inline] |
Definition at line 127 of file BCEngineMCMC.h.
{ return fMCMCNIterationsRun; }
int BCEngineMCMC::MCMCGetNIterationsUpdate | ( | ) | [inline] |
Definition at line 137 of file BCEngineMCMC.h.
{ return fMCMCNIterationsUpdate; }
int BCEngineMCMC::MCMCGetNIterationsUpdateMax | ( | ) | [inline] |
Definition at line 142 of file BCEngineMCMC.h.
{ return fMCMCNIterationsUpdateMax; }
int BCEngineMCMC::MCMCGetNLag | ( | ) | [inline] |
Definition at line 91 of file BCEngineMCMC.h.
{ return fMCMCNLag; }
int BCEngineMCMC::MCMCGetNParameters | ( | ) | [inline] |
Definition at line 81 of file BCEngineMCMC.h.
{ return fMCMCNParameters; }
std::vector<int> BCEngineMCMC::MCMCGetNTrialsFalse | ( | ) | [inline] |
Definition at line 152 of file BCEngineMCMC.h.
{ return fMCMCNTrialsFalse; }
std::vector<int> BCEngineMCMC::MCMCGetNTrialsTrue | ( | ) | [inline] |
Definition at line 147 of file BCEngineMCMC.h.
{ return fMCMCNTrialsTrue; }
int BCEngineMCMC::MCMCGetPhase | ( | ) | [inline] |
Definition at line 211 of file BCEngineMCMC.h.
{ return fMCMCPhase; }
std::vector<double> BCEngineMCMC::MCMCGetprobMean | ( | ) | [inline] |
Definition at line 158 of file BCEngineMCMC.h.
{ return fMCMCprobMean; }
bool BCEngineMCMC::MCMCGetProposalPointMetropolis | ( | int | chain, | |
std::vector< double > & | x | |||
) |
Returns a trial point for the Metropolis algorithm.
chain | chain index | |
x | proposal point |
Definition at line 584 of file BCEngineMCMC.cxx.
{ // get unscaled random point. this point might not be in the correct volume. MCMCTrialFunction(chain, x); // get a proposal point from the trial function and scale it for (int i = 0; i < fMCMCNParameters; ++i) x[i] = fMCMCx[chain * fMCMCNParameters + i] + x[i] * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i)); // check if the point is in the correct volume. for (int i = 0; i < fMCMCNParameters; ++i) if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i])) return false; return true; }
bool BCEngineMCMC::MCMCGetProposalPointMetropolis | ( | int | chain, | |
int | parameter, | |||
std::vector< double > & | x | |||
) |
Returns a trial point for the Metropolis algorithm.
chain | chain index | |
x | proposal point |
Definition at line 602 of file BCEngineMCMC.cxx.
{ // get unscaled random point in the dimension of the chosen // parameter. this point might not be in the correct volume. double proposal = MCMCTrialFunctionSingle(ichain, ipar); // copy the old point into the new for (int i = 0; i < fMCMCNParameters; ++i) x[i] = fMCMCx[ichain * fMCMCNParameters + i]; // modify the parameter under study x[ipar] += proposal * (fMCMCBoundaryMax[ipar] - fMCMCBoundaryMin[ipar]); // check if the point is in the correct volume. if ((x[ipar] < fMCMCBoundaryMin[ipar]) || (x[ipar] > fMCMCBoundaryMax[ipar])) return false; return true; }
double BCEngineMCMC::MCMCGetRValue | ( | ) | [inline] |
double BCEngineMCMC::MCMCGetRValueCriterion | ( | ) | [inline] |
Definition at line 241 of file BCEngineMCMC.h.
{ return fMCMCRValueCriterion; }
double BCEngineMCMC::MCMCGetRValueParameters | ( | int | i | ) | [inline] |
i | parameter index |
Definition at line 257 of file BCEngineMCMC.h.
{ return fMCMCRValueParameters.at(i); }
double BCEngineMCMC::MCMCGetRValueParametersCriterion | ( | ) | [inline] |
Definition at line 246 of file BCEngineMCMC.h.
{ return fMCMCRValueParametersCriterion; }
bool BCEngineMCMC::MCMCGetRValueStrict | ( | ) | [inline] |
Use strict or relaxed rule for Gelman/Rubin R-value
Definition at line 261 of file BCEngineMCMC.h.
{ return fMCMCRValueUseStrict; }
TRandom3* BCEngineMCMC::MCMCGetTRandom3 | ( | ) | [inline] |
Return the random number generator. Any non-zero seed gives reproducible behavior,e.g. m->MCMCGetTRandom3()->SetSeed(21340)
Definition at line 293 of file BCEngineMCMC.h.
{ return fRandom; }
std::vector<double> BCEngineMCMC::MCMCGetTrialFunctionScaleFactor | ( | ) | [inline] |
Definition at line 169 of file BCEngineMCMC.h.
{ return fMCMCTrialFunctionScaleFactor; }
std::vector< double > BCEngineMCMC::MCMCGetTrialFunctionScaleFactor | ( | int | ichain | ) |
ichain | chain index |
Definition at line 509 of file BCEngineMCMC.cxx.
{ // create a new vector with the length of fMCMCNParameters std::vector<double> x; // check if ichain is in range if (ichain < 0 || ichain >= fMCMCNChains) return x; // copy the scale factors into the temporary vector for (int j = 0; j < fMCMCNParameters; ++j) x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNParameters + j)); return x; }
double BCEngineMCMC::MCMCGetTrialFunctionScaleFactor | ( | int | ichain, | |
int | ipar | |||
) |
ichain | chain index | |
ipar | parameter index |
Definition at line 526 of file BCEngineMCMC.cxx.
{ // check if ichain is in range if (ichain < 0 || ichain >= fMCMCNChains) return 0; // check if ipar is in range if (ipar < 0 || ipar >= fMCMCNParameters) return 0; // return component of ipar point in the ichain chain return fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNChains + ipar); }
std::vector<double> BCEngineMCMC::MCMCGetVariance | ( | ) | [inline] |
Definition at line 164 of file BCEngineMCMC.h.
{ return fMCMCprobVar; }
std::vector<double> BCEngineMCMC::MCMCGetx | ( | ) | [inline] |
Definition at line 185 of file BCEngineMCMC.h.
{ return fMCMCx; }
std::vector< double > BCEngineMCMC::MCMCGetx | ( | int | ichain | ) |
ichain | index of the Markov chain |
Definition at line 541 of file BCEngineMCMC.cxx.
{ // create a new vector with the length of fMCMCNParameters std::vector<double> x; // check if ichain is in range if (ichain < 0 || ichain >= fMCMCNChains) return x; // copy the point in the ichain chain into the temporary vector for (int j = 0; j < fMCMCNParameters; ++j) x.push_back(fMCMCx.at(ichain * fMCMCNParameters + j)); return x; }
double BCEngineMCMC::MCMCGetx | ( | int | ichain, | |
int | ipar | |||
) |
ichain | chain index | |
ipar | parameter index |
Definition at line 558 of file BCEngineMCMC.cxx.
{ // check if ichain is in range if (ichain < 0 || ichain >= fMCMCNChains) return 0; // check if ipar is in range if (ipar < 0 || ipar >= fMCMCNParameters) return 0; // return component of jth point in the ith chain return fMCMCx.at(ichain * fMCMCNParameters + ipar); }
void BCEngineMCMC::MCMCInChainCheckMaximum | ( | ) |
Updates statistics: find new maximum
Definition at line 767 of file BCEngineMCMC.cxx.
{ // loop over all chains for (int i = 0; i < fMCMCNChains; ++i) { // check if new maximum is found or chain is at the beginning if (fMCMCprob[i] > fMCMCprobMax[i] || fMCMCNIterations[i] == 1) { // copy maximum value fMCMCprobMax[i] = fMCMCprob[i]; // copy mode of chain for (int j = 0; j < fMCMCNParameters; ++j) fMCMCxMax[i * fMCMCNParameters + j] = fMCMCx[i * fMCMCNParameters + j]; } } }
void BCEngineMCMC::MCMCInChainFillHistograms | ( | ) |
Updates statistics: fill marginalized distributions
Definition at line 819 of file BCEngineMCMC.cxx.
{ // check if histograms are supposed to be filled if (!fMCMCFlagFillHistograms) return; // loop over chains for (int i = 0; i < fMCMCNChains; ++i) { // fill each 1-dimensional histogram (if supposed to be filled) for (int j = 0; j < fMCMCNParameters; ++j) if (fMCMCFlagsFillHistograms.at(j)) fMCMCH1Marginalized[j]->Fill(fMCMCx[i * fMCMCNParameters + j]); // fill each 2-dimensional histogram (if supposed to be filled) int counter = 0; for (int j = 0; j < fMCMCNParameters; ++j) for (int k = 0; k < j; ++k) { if (fMCMCFlagsFillHistograms.at(j) && fMCMCFlagsFillHistograms.at(k)) fMCMCH2Marginalized[counter]->Fill(fMCMCx[i*fMCMCNParameters+k],fMCMCx[i* fMCMCNParameters+j]); counter ++; } } }
void BCEngineMCMC::MCMCInChainTestConvergenceAllChains | ( | ) |
Updates statistics: check convergence
Definition at line 847 of file BCEngineMCMC.cxx.
{ // calculate number of entries in this part of the chain int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0]; if (fMCMCNChains > 1 && npoints > 1) { // define flag for convergence bool flag_convergence = true; // loop over parameters for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters){ // extract means and variances std::vector<double> means(fMCMCNChains); std::vector<double> variances(fMCMCNChains); // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) { // get parameter index int index = ichains * fMCMCNParameters + iparameters; means[ichains] = fMCMCxMean[index]; variances[ichains] = fMCMCxVar[index]; } fMCMCRValueParameters[iparameters] = BCMath::Rvalue(means, variances, npoints, fMCMCRValueUseStrict); // set flag to false if convergence criterion is not fulfilled for the parameter if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion)) flag_convergence = false; // else: leave convergence flag true for that parameter } fMCMCRValue = BCMath::Rvalue(fMCMCprobMean, fMCMCprobVar, npoints, fMCMCRValueUseStrict); // set flag to false if convergence criterion is not fulfilled for the log-likelihood if (!((fMCMCRValue - 1.0) < fMCMCRValueCriterion)) flag_convergence = false; // remember number of iterations needed to converge if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true) fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters; } }
void BCEngineMCMC::MCMCInChainUpdateStatistics | ( | ) |
Updates statistics:
Definition at line 786 of file BCEngineMCMC.cxx.
{ // calculate number of entries in this part of the chain int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0]; // length of vectors int nentries = fMCMCNParameters * fMCMCNChains; // loop over all parameters of all chains for (int i = 0; i < nentries; ++i) { // calculate mean value of each parameter in the chain for this part fMCMCxMean[i] += (fMCMCx[i] - fMCMCxMean[i]) / double(npoints); // calculate variance of each chain for this part if (npoints > 1) fMCMCxVar[i] = (1.0 - 1./double(npoints)) * fMCMCxVar[i] + (fMCMCx[i] - fMCMCxMean[i]) * (fMCMCx[i] - fMCMCxMean[i]) / double(npoints - 1); } // loop over chains for (int i = 0; i < fMCMCNChains; ++i) { // calculate mean value of each chain for this part fMCMCprobMean[i] += (fMCMCprob[i] - fMCMCprobMean[i]) / double(npoints); // calculate variance of each chain for this part if (npoints > 1) fMCMCprobVar[i] = (1.0 - 1/double(npoints)) * fMCMCprobVar[i] + (fMCMCprob[i] - fMCMCprobMean[i]) * (fMCMCprob[i] - fMCMCprobMean[i]) / double(npoints - 1); } }
void BCEngineMCMC::MCMCInChainWriteChains | ( | ) |
Updates statistics: write chains to file
Definition at line 892 of file BCEngineMCMC.cxx.
{ // loop over all chains for (int i = 0; i < fMCMCNChains; ++i) fMCMCTrees[i]->Fill(); }
int BCEngineMCMC::MCMCInitialize | ( | ) |
Initializes the engine.
Definition at line 1548 of file BCEngineMCMC.cxx.
{ // reset values MCMCResetResults(); // free memory for vectors fMCMCNIterations.assign(fMCMCNChains, 0); fMCMCprobMean.assign(fMCMCNChains, 0); fMCMCprobVar.assign(fMCMCNChains, 0); fMCMCprob.assign(fMCMCNChains, -1.0); fMCMCprobMax.assign(fMCMCNChains, -1.0); fMCMCNTrialsTrue.assign(fMCMCNChains * fMCMCNParameters, 0); fMCMCNTrialsFalse.assign(fMCMCNChains * fMCMCNParameters, 0); fMCMCxMax.assign(fMCMCNChains * fMCMCNParameters, 0.); fMCMCxMean.assign(fMCMCNChains * fMCMCNParameters, 0); fMCMCxVar.assign(fMCMCNChains * fMCMCNParameters, 0); fMCMCRValueParameters.assign(fMCMCNParameters, 0.); if (fMCMCTrialFunctionScaleFactorStart.size() == 0) fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fMCMCNParameters, 1.0); else for (int i = 0; i < fMCMCNChains; ++i) for (int j = 0; j < fMCMCNParameters; ++j) fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j)); // set initial position if (fMCMCFlagInitialPosition == 2) // user defined points { // define flag bool flag = true; // check the length of the array of initial positions if (int(fMCMCInitialPosition.size()) != (fMCMCNChains * fMCMCNParameters)) { BCLog::OutError("BCEngine::MCMCInitialize : Length of vector containing initial positions does not have required length."); flag = false; } // check the boundaries if (flag) { for (int j = 0; j < fMCMCNChains; ++j) for (int i = 0; i < fMCMCNParameters; ++i) if (fMCMCInitialPosition[j * fMCMCNParameters + i] < fMCMCBoundaryMin[i] || fMCMCInitialPosition[j * fMCMCNParameters + i] > fMCMCBoundaryMax[i]) { BCLog::OutError("BCEngine::MCMCInitialize : Initial position out of boundaries."); flag = false; } } // check flag if (!flag) fMCMCFlagInitialPosition = 1; } if (fMCMCFlagInitialPosition == 0) // center of the interval for (int j = 0; j < fMCMCNChains; ++j) for (int i = 0; i < fMCMCNParameters; ++i) fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i])); else if (fMCMCFlagInitialPosition == 2) // user defined { for (int j = 0; j < fMCMCNChains; ++j) for (int i = 0; i < fMCMCNParameters; ++i) fMCMCx.push_back(fMCMCInitialPosition.at(j * fMCMCNParameters + i)); } else for (int j = 0; j < fMCMCNChains; ++j) // random number (default) for (int i = 0; i < fMCMCNParameters; ++i) fMCMCx.push_back(fMCMCBoundaryMin[i] + fRandom->Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i])); // copy the point of the first chain fMCMCxLocal.clear(); for (int i = 0; i < fMCMCNParameters; ++i) fMCMCxLocal.push_back(fMCMCx[i]); // define 1-dimensional histograms for marginalization for(int i = 0; i < fMCMCNParameters; ++i) { TH1D * h1 = 0; if (fMCMCFlagsFillHistograms.at(i)) h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i), "", fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i]); fMCMCH1Marginalized.push_back(h1); } for(int i = 0; i < fMCMCNParameters; ++i) for (int k = 0; k < i; ++k) { TH2D * h2 = 0; if (fMCMCFlagsFillHistograms.at(i) && fMCMCFlagsFillHistograms.at(k)) h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, k), "", fMCMCH1NBins[k], fMCMCBoundaryMin[k], fMCMCBoundaryMax[k], fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i] ); fMCMCH2Marginalized.push_back(h2); } return 1; }
void BCEngineMCMC::MCMCInitializeMarkovChains | ( | ) |
Initializes Markov chains.
Definition at line 1491 of file BCEngineMCMC.cxx.
{ // evaluate function at the starting point std::vector<double> x0; for (int j = 0; j < fMCMCNChains; ++j) { x0.clear(); for (int i = 0; i < fMCMCNParameters; ++i) x0.push_back(fMCMCx[j * fMCMCNParameters + i]); fMCMCprob[j] = LogEval(x0); } x0.clear(); }
void BCEngineMCMC::MCMCInitializeMarkovChainTrees | ( | ) |
Initialize trees containing the Markov chains.
Definition at line 467 of file BCEngineMCMC.cxx.
{ // clear vector fMCMCTrees.clear(); // create new trees for (int i = 0; i < fMCMCNChains; ++i) { fMCMCTrees.push_back(new TTree(TString::Format("MarkovChainTree_%i", i), "MarkovChainTree")); fMCMCTrees[i]->Branch("Iteration", &fMCMCNIterations[i], "iteration/I"); fMCMCTrees[i]->Branch("NParameters", &fMCMCNParameters, "parameters/I"); fMCMCTrees[i]->Branch("LogProbability", &fMCMCprob[i], "log(probability)/D"); fMCMCTrees[i]->Branch("Phase", &fMCMCPhase, "phase/I"); fMCMCTrees[i]->Branch("Cycle", &fMCMCCycle, "cycle/I"); for (int j = 0; j < fMCMCNParameters; ++j) fMCMCTrees[i]->Branch(TString::Format("Parameter%i", j), &fMCMCx[i * fMCMCNParameters + j], TString::Format("parameter %i/D", j)); } }
virtual void BCEngineMCMC::MCMCIterationInterface | ( | ) | [inline, virtual] |
Interface allowing to execute arbitrary code for each iteration of the MCMC. The frequency of calling this method is influenced by the setup of the Lag and whether or not the MCMC is run with ordered parameters. This method needs to be overloaded in the derived class.
Reimplemented in BCIntegrate, and BCRooInterface.
Definition at line 553 of file BCEngineMCMC.h.
{}
int BCEngineMCMC::MCMCMetropolis | ( | ) |
Runs Metropolis algorithm.
Definition at line 1296 of file BCEngineMCMC.cxx.
{ // check if prerun should be performed if (fMCMCFlagPreRun) MCMCMetropolisPreRun(); else BCLog::OutWarning("BCEngineMCMC::MCMCMetropolis. Not running prerun. This can cause trouble if the data have changed."); // print to screen BCLog::OutSummary( "Run Metropolis MCMC..."); // reset run statistics MCMCResetRunStatistics(); // set phase and cycle number fMCMCPhase = 2; fMCMCCycle = 0; // perform run BCLog::OutSummary(Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun)); // int counterupdate = 0; // bool convergence = false; // bool flagefficiency = false; // std::vector<double> efficiency; // for (int i = 0; i < fMCMCNParameters; ++i) // for (int j = 0; j < fMCMCNChains; ++j) // efficiency.push_back(0.0); int nwrite = fMCMCNIterationsRun/10; if(nwrite < 100) nwrite=100; else if(nwrite < 500) nwrite=1000; else if(nwrite < 10000) nwrite=1000; else nwrite=10000; // start the run for (fMCMCCurrentIteration = 1; fMCMCCurrentIteration <= fMCMCNIterationsRun; ++fMCMCCurrentIteration) { if ( (fMCMCCurrentIteration)%nwrite == 0 ) BCLog::OutDetail(Form(" --> iteration number %i (%.2f%%)", fMCMCCurrentIteration, (double)(fMCMCCurrentIteration)/(double)fMCMCNIterationsRun*100.)); // if the flag is set then run over the parameters one after the other. if (fMCMCFlagOrderParameters) { // loop over parameters for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters) { // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) MCMCGetNewPointMetropolis(ichains, iparameters); // reset current chain fMCMCCurrentChain = -1; // update search for maximum MCMCInChainCheckMaximum(); } // end loop over all parameters // check if the current iteration is consistent with the lag if ( fMCMCCurrentIteration % fMCMCNLag == 0) { // fill histograms MCMCInChainFillHistograms(); // write chain to file if (fMCMCFlagWriteChainToFile) MCMCInChainWriteChains(); // do anything interface MCMCIterationInterface(); } } // if the flag is not set then run over the parameters at the same time. else { // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) // get new point MCMCGetNewPointMetropolis(ichains); // reset current chain fMCMCCurrentChain = -1; // update search for maximum MCMCInChainCheckMaximum(); // check if the current iteration is consistent with the lag if (fMCMCCurrentIteration % fMCMCNLag == 0) { // fill histograms MCMCInChainFillHistograms(); // write chain to file if (fMCMCFlagWriteChainToFile) MCMCInChainWriteChains(); // do anything interface MCMCIterationInterface(); } } } // end run // print convergence status BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun)); // print modes // find global maximum double probmax = fMCMCprobMax.at(0); int probmaxindex = 0; // loop over all chains and find the maximum point for (int i = 1; i < fMCMCNChains; ++i) { if (fMCMCprobMax.at(i) > probmax) { probmax = fMCMCprobMax.at(i); probmaxindex = i; } } BCLog::OutDetail(" --> Global mode from MCMC:"); BCLog::OutDebug(Form(" --> Posterior value: %g", probmax)); int ndigits = (int) log10(fMCMCNParameters); for (int i = 0; i < fMCMCNParameters; ++i) BCLog::OutDetail(Form( TString::Format(" --> parameter %%%di: %%.4g", ndigits+1), i, fMCMCxMax[probmaxindex * fMCMCNParameters + i])); // reset counter fMCMCCurrentIteration = -1; // reset current chain fMCMCCurrentChain = -1; // set flags fMCMCFlagRun = true; return 1; }
int BCEngineMCMC::MCMCMetropolisPreRun | ( | ) |
Runs a pre run for the Metropolis algorithm.
Definition at line 908 of file BCEngineMCMC.cxx.
{ // print on screen BCLog::OutSummary("Pre-run Metropolis MCMC..."); // initialize Markov chain MCMCInitialize(); MCMCInitializeMarkovChains(); // helper variable containing number of digits in the number of parameters int ndigits = (int)log10(fMCMCNParameters) +1; if(ndigits<4) ndigits=4; // reset run statistics MCMCResetRunStatistics(); fMCMCNIterationsConvergenceGlobal = -1; // perform run BCLog::OutSummary(Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax)); // don't write to file during pre run bool tempflag_writetofile = fMCMCFlagWriteChainToFile; fMCMCFlagWriteChainToFile = false; // initialize counter variables and flags fMCMCCurrentIteration = 1; // counts the number of iterations int counterupdate = 1; // after how many iterations is an update needed? bool convergence = false; // convergence reached? bool flagefficiency = false; // efficiency reached? // array of efficiencies std::vector<double> efficiency; efficiency.assign(fMCMCNParameters * fMCMCNChains, 0.0); // how often to check convergence and efficiencies? // it's either every fMCMCNParameters*nMCMCNIterationsUpdate (for 5 parameters the default would be 5000) // or it's fMCMCNIterationsUpdateMax (10000 by default) // whichever of the two is smaller int updateLimit = ( fMCMCNIterationsUpdateMax<fMCMCNIterationsUpdate*(fMCMCNParameters) && fMCMCNIterationsUpdateMax>0 ) ? fMCMCNIterationsUpdateMax : fMCMCNIterationsUpdate*(fMCMCNParameters); // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) { // loop over parameters for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter){ // global index of the parameter (throughout all the chains) int index = ichains * fMCMCNParameters + iparameter; // reset counters fMCMCNTrialsTrue[index] = 0; fMCMCNTrialsFalse[index] = 0; fMCMCxMean[index] = fMCMCx[index]; fMCMCxVar[index] = 0; } fMCMCprobMean[ichains] = fMCMCprob[ichains]; fMCMCprobVar[ichains] = 0; } // set phase and cycle number fMCMCPhase = 1; fMCMCCycle = 0; // run chain ... // (a) for at least a minimum number of iterations, // (b) until a maximum number of iterations is reached, // (c) or until convergence is reached and the efficiency is in the // specified region while (fMCMCCurrentIteration < fMCMCNIterationsPreRunMin || (fMCMCCurrentIteration < fMCMCNIterationsMax && !(convergence && flagefficiency))) { //------------------------------------------- // reset flags and counters //------------------------------------------- // set convergence to false by default convergence = false; // set number of iterations needed to converge to negative fMCMCNIterationsConvergenceGlobal = -1; //------------------------------------------- // get new point in n-dim space //------------------------------------------- // if the flag is set then run over the parameters one after the other. if (fMCMCFlagOrderParameters) { // loop over parameters for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters) { // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) MCMCGetNewPointMetropolis(ichains, iparameters); // search for global maximum MCMCInChainCheckMaximum(); } } // if the flag is not set then run over the parameters at the same time. else { // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) // get new point MCMCGetNewPointMetropolis(ichains); // search for global maximum MCMCInChainCheckMaximum(); } //------------------------------------------- // print out message to log //------------------------------------------- // progress printout if ( fMCMCCurrentIteration > 0 && fMCMCCurrentIteration % fMCMCNIterationsUpdate == 0 ) BCLog::OutDetail(Form(" --> Iteration %i", fMCMCNIterations[0]/fMCMCNParameters)); //------------------------------------------- // update statistics //------------------------------------------- if (counterupdate > 1) MCMCInChainUpdateStatistics(); //------------------------------------------- // update scale factors and check convergence //------------------------------------------- // debugKK // check if this line makes sense if ( counterupdate % updateLimit == 0 && counterupdate > 0 && fMCMCCurrentIteration >= fMCMCNIterationsPreRunMin) // if ( fMCMCCurrentIteration % fMCMCNIterationsUpdate == 0 && counterupdate > 1 && fMCMCCurrentIteration >= fMCMCNIterationsPreRunMin) { // ----------------------------- // reset flags and counters // ----------------------------- bool rvalues_ok = true; static bool has_converged = false; // reset the number of iterations needed for convergence to // negative fMCMCNIterationsConvergenceGlobal = -1; // ----------------------------- // check convergence status // ----------------------------- // test convergence MCMCInChainTestConvergenceAllChains(); // set convergence flag if (fMCMCNIterationsConvergenceGlobal > 0) convergence = true; // print convergence status: if (convergence && fMCMCNChains > 1) BCLog::OutDetail(Form(" * Convergence status: Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal)); else if (!convergence && fMCMCNChains > 1) { BCLog::OutDetail(Form(" * Convergence status: Set of %i Markov chains did not converge after %i iterations.", fMCMCNChains, fMCMCCurrentIteration)); BCLog::OutDetail(" - R-Values:"); for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter) { if(fabs(fMCMCRValueParameters[iparameter]-1.) < fMCMCRValueParametersCriterion) BCLog::OutDetail(Form( TString::Format(" parameter %%%di : %%.06f",ndigits), iparameter, fMCMCRValueParameters.at(iparameter))); else { BCLog::OutDetail(Form( TString::Format(" parameter %%%di : %%.06f <--",ndigits), iparameter, fMCMCRValueParameters.at(iparameter))); rvalues_ok = false; } } if(fabs(fMCMCRValue-1.) < fMCMCRValueCriterion) BCLog::OutDetail(Form(" log-likelihood : %.06f", fMCMCRValue)); else { BCLog::OutDetail(Form(" log-likelihood : %.06f <--", fMCMCRValue)); rvalues_ok = false; } } // set convergence flag if(!has_converged) if(rvalues_ok) has_converged = true; // ----------------------------- // check efficiency status // ----------------------------- // set flag flagefficiency = true; bool flagprintefficiency = true; // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) { // loop over parameters for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter) { // global index of the parameter (throughout all the chains) int index = ichains * fMCMCNParameters + iparameter; // calculate efficiency efficiency[index] = double(fMCMCNTrialsTrue[index]) / double(fMCMCNTrialsTrue[index] + fMCMCNTrialsFalse[index]); // adjust scale factors if efficiency is too low if (efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01) { if (flagprintefficiency) { BCLog::OutDetail(Form(" * Efficiency status: Efficiencies not within pre-defined range.")); BCLog::OutDetail(Form(" - Efficiencies:")); flagprintefficiency = false; } double fscale=2.; if(has_converged && fMCMCEfficiencyMin/efficiency[index] > 2.) fscale = 4.; fMCMCTrialFunctionScaleFactor[index] /= fscale; BCLog::OutDetail(Form(" Efficiency of parameter %i dropped below %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g", iparameter, 100. * fMCMCEfficiencyMin, 100. * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index])); } // adjust scale factors if efficiency is too high else if (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.0) { if (flagprintefficiency) { BCLog::OutDetail(Form(" * Efficiency status: Efficiencies not within pre-defined ranges.")); BCLog::OutDetail(Form(" - Efficiencies:")); flagprintefficiency = false; } fMCMCTrialFunctionScaleFactor[index] *= 2.; BCLog::OutDetail(Form(" Efficiency of parameter %i above %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g", iparameter, 100.0 * fMCMCEfficiencyMax, 100.0 * efficiency[index], ichains, fMCMCTrialFunctionScaleFactor[index])); } // check flag if ((efficiency[index] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[index] > .01) || (efficiency[index] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[index] < 1.)) flagefficiency = false; } // end of running over all parameters } // end of running over all chains // print to screen if (flagefficiency) BCLog::OutDetail(Form(" * Efficiency status: Efficiencies within pre-defined ranges.")); // ----------------------------- // reset counters // ----------------------------- counterupdate = 0; // loop over chains for (int ichains = 0; ichains < fMCMCNChains; ++ichains) { // loop over parameters for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter){ // global index of the parameter (throughout all the chains) int index = ichains * fMCMCNParameters + iparameter; // reset counters fMCMCNTrialsTrue[index] = 0; fMCMCNTrialsFalse[index] = 0; fMCMCxMean[index] = fMCMCx[index]; fMCMCxVar[index] = 0; } fMCMCprobMean[ichains] = fMCMCprob[ichains]; fMCMCprobVar[ichains] = 0; } } // end if update scale factors and check convergence //------------------------------------------- // write chain to file //------------------------------------------- // write chain to file if (fMCMCFlagWritePreRunToFile) MCMCInChainWriteChains(); //------------------------------------------- // increase counters //------------------------------------------- if (counterupdate == 0 && fMCMCCurrentIteration > 1) fMCMCCycle++; fMCMCCurrentIteration++; counterupdate++; } // end of running // decrease counter by one since it didn't really run that long // fMCMCCurrentIteration--; // counterupdate--; // did we check convergence at least once ? if (fMCMCCurrentIteration<updateLimit) { BCLog::OutWarning(" Convergence never checked !"); BCLog::OutWarning(" Increase maximum number of iterations in the pre-run /MCMCSetNIterationsMax()/"); BCLog::OutWarning(" or decrease maximum number of iterations for update /MCMCSetNIterationsUpdateMax()/"); } // --------------- // after chain run // --------------- // define convergence status if (fMCMCNIterationsConvergenceGlobal > 0) fMCMCFlagConvergenceGlobal = true; else fMCMCFlagConvergenceGlobal = false; // print convergence status if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && !flagefficiency) BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations but could not adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal)); else if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && flagefficiency) BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations and all scales are adjusted.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal)); else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && flagefficiency) BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax)); else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && !flagefficiency) BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations and could not adjust scales.", fMCMCNChains, fMCMCNIterationsMax)); else if(fMCMCNChains == 1) BCLog::OutSummary(" --> No convergence criterion for a single chain defined."); else BCLog::OutSummary(" --> Only one Markov chain. No global convergence criterion defined."); BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCCurrentIteration)); // print efficiencies std::vector<double> efficiencies; for (int i = 0; i < fMCMCNParameters; ++i) efficiencies.push_back(0.); BCLog::OutDetail(" --> Average efficiencies:"); for (int i = 0; i < fMCMCNParameters; ++i) { for (int j = 0; j < fMCMCNChains; ++j) efficiencies[i] += efficiency[j * fMCMCNParameters + i] / double(fMCMCNChains); BCLog::OutDetail(Form(TString::Format(" --> parameter %%%di : %%.02f%%%%",ndigits), i, 100. * efficiencies[i])); } // print scale factors std::vector<double> scalefactors; for (int i = 0; i < fMCMCNParameters; ++i) scalefactors.push_back(0.0); BCLog::OutDetail(" --> Average scale factors:"); for (int i = 0; i < fMCMCNParameters; ++i) { for (int j = 0; j < fMCMCNChains; ++j) scalefactors[i] += fMCMCTrialFunctionScaleFactor[j * fMCMCNParameters + i] / double(fMCMCNChains); BCLog::OutDetail(Form( TString::Format(" --> parameter %%%di : %%.02f%%%%",ndigits), i, 100. * scalefactors[i])); } // reset flag fMCMCFlagWriteChainToFile = tempflag_writetofile; // reset current iteration fMCMCCurrentIteration = -1; // reset current chain fMCMCCurrentChain = -1; // no error return 1; }
int BCEngineMCMC::MCMCResetResults | ( | ) |
Reset the MCMC variables.
Definition at line 1508 of file BCEngineMCMC.cxx.
{ // reset variables fMCMCNIterations.clear(); fMCMCNTrialsTrue.clear(); fMCMCNTrialsFalse.clear(); fMCMCTrialFunctionScaleFactor.clear(); fMCMCprobMean.clear(); fMCMCprobVar.clear(); fMCMCxMean.clear(); fMCMCxVar.clear(); fMCMCx.clear(); fMCMCprob.clear(); fMCMCxMax.clear(); fMCMCprobMax.clear(); fMCMCNIterationsConvergenceGlobal = -1; fMCMCRValueParameters.clear(); for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i) if (fMCMCH1Marginalized[i]) delete fMCMCH1Marginalized[i]; for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i) if (fMCMCH2Marginalized[i]) delete fMCMCH2Marginalized[i]; // clear plots fMCMCH1Marginalized.clear(); fMCMCH2Marginalized.clear(); // reset flags fMCMCFlagPreRun = true; fMCMCFlagRun = false; fMCMCFlagConvergenceGlobal = false; // no errors return 1; }
void BCEngineMCMC::MCMCResetRunStatistics | ( | ) |
Resets the run statistics.
Definition at line 1444 of file BCEngineMCMC.cxx.
{ for (int j = 0; j < fMCMCNChains; ++j) { fMCMCNIterations[j] = 0; fMCMCNTrialsTrue[j] = 0; fMCMCNTrialsFalse[j] = 0; fMCMCprobMean[j] = 0; fMCMCprobVar[j] = 0; for (int k = 0; k < fMCMCNParameters; ++k) { fMCMCNTrialsTrue[j * fMCMCNParameters + k] = 0; fMCMCNTrialsFalse[j * fMCMCNParameters + k] = 0; } } // reset marginalized distributions for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i) if (fMCMCH1Marginalized[i]) fMCMCH1Marginalized[i]->Reset(); for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i) if (fMCMCH2Marginalized[i]) fMCMCH2Marginalized[i]->Reset(); fMCMCRValue = 100; }
void BCEngineMCMC::MCMCSetFlagFillHistograms | ( | bool | flag | ) |
Sets the flag for all parameters to either fill histograms or not.
Definition at line 433 of file BCEngineMCMC.cxx.
{ fMCMCFlagFillHistograms = flag; for (int i = 0; i < fMCMCNParameters; ++i) fMCMCFlagsFillHistograms[i] = flag; }
void BCEngineMCMC::MCMCSetFlagFillHistograms | ( | int | index, | |
bool | flag | |||
) |
Sets the flag for a single parameter to either fill histograms or not.
Definition at line 442 of file BCEngineMCMC.cxx.
{ // check if index is within range if (index < 0 || index > fMCMCNParameters) { BCLog::OutWarning("BCEngineMCMC :MCMCSetFlagFillHistograms. Index out of range."); return; } // set flag fMCMCFlagsFillHistograms[index] = flag; }
void BCEngineMCMC::MCMCSetFlagInitialPosition | ( | int | flag | ) | [inline] |
Sets flag which defines initial position.
Definition at line 377 of file BCEngineMCMC.h.
{ fMCMCFlagInitialPosition = flag; }
void BCEngineMCMC::MCMCSetFlagOrderParameters | ( | bool | flag | ) | [inline] |
Sets the flag which controls the sequence parameters during the running of the MCMC.
Definition at line 383 of file BCEngineMCMC.h.
{ fMCMCFlagOrderParameters = flag; }
void BCEngineMCMC::MCMCSetFlagPreRun | ( | bool | flag | ) | [inline] |
Sets the flag if a prerun should be performed or not.
Definition at line 393 of file BCEngineMCMC.h.
{ fMCMCFlagPreRun = flag; }
void BCEngineMCMC::MCMCSetInitialPositions | ( | std::vector< std::vector< double > > | x0s | ) |
Sets the initial positions for all chains.
x0s | initial positions for all chains. |
Definition at line 419 of file BCEngineMCMC.cxx.
{ // create new vector std::vector<double> y0s; // loop over vector elements for (int i = 0; i < int(x0s.size()); ++i) for (int j = 0; j < int((x0s.at(i)).size()); ++j) y0s.push_back((x0s.at(i)).at(j)); MCMCSetInitialPositions(y0s); }
void BCEngineMCMC::MCMCSetInitialPositions | ( | const std::vector< double > & | x0s | ) |
Sets the initial positions for all chains.
x0s | initial positions for all chains. |
Definition at line 403 of file BCEngineMCMC.cxx.
{ // clear the existing initial position vector fMCMCInitialPosition.clear(); // copy the initial positions int n = int(x0s.size()); for (int i = 0; i < n; ++i) fMCMCInitialPosition.push_back(x0s.at(i)); // use these intial positions for the Markov chain MCMCSetFlagInitialPosition(2); }
void BCEngineMCMC::MCMCSetMarkovChainTrees | ( | std::vector< TTree * > | trees | ) |
Sets the tree containing the Markov chains.
Definition at line 456 of file BCEngineMCMC.cxx.
{ // clear vector fMCMCTrees.clear(); // copy tree for (int i = 0; i < int(trees.size()); ++i) fMCMCTrees.push_back(trees[i]); }
void BCEngineMCMC::MCMCSetMaximumEfficiency | ( | double | efficiency | ) | [inline] |
Sets the maximum efficiency required for a chain.
Definition at line 352 of file BCEngineMCMC.h.
{ fMCMCEfficiencyMax = efficiency; }
void BCEngineMCMC::MCMCSetMinimumEfficiency | ( | double | efficiency | ) | [inline] |
Sets the minimum efficiency required for a chain.
Definition at line 347 of file BCEngineMCMC.h.
{ fMCMCEfficiencyMin = efficiency; }
void BCEngineMCMC::MCMCSetNChains | ( | int | n | ) |
Sets the number of Markov chains which are run in parallel.
Definition at line 394 of file BCEngineMCMC.cxx.
{ fMCMCNChains = n; // re-initialize MCMCInitialize(); }
void BCEngineMCMC::MCMCSetNIterationsMax | ( | int | n | ) | [inline] |
Sets the maximum number of iterations in the pre-run.
Definition at line 317 of file BCEngineMCMC.h.
{ fMCMCNIterationsMax = n; }
void BCEngineMCMC::MCMCSetNIterationsPreRunMin | ( | int | n | ) | [inline] |
Sets the minimum number of iterations in the pre-run
Definition at line 327 of file BCEngineMCMC.h.
{ fMCMCNIterationsPreRunMin = n; }
void BCEngineMCMC::MCMCSetNIterationsRun | ( | int | n | ) | [inline] |
Sets the number of iterations.
Definition at line 322 of file BCEngineMCMC.h.
{ fMCMCNIterationsRun = n; }
void BCEngineMCMC::MCMCSetNIterationsUpdate | ( | int | n | ) | [inline] |
Sets the number of iterations in the pre-run after which an update on the statistics (convergence, efficiency, etc.) is done.
n | The number of iterations. |
Definition at line 334 of file BCEngineMCMC.h.
{ fMCMCNIterationsUpdate = n; }
void BCEngineMCMC::MCMCSetNIterationsUpdateMax | ( | int | n | ) | [inline] |
Sets the maximum number of iterations in the pre-run after which an update on the statistics (convergence, efficiency, etc.) is done. If set to 0 no maximum is set.
n | maximum number of iterations. |
Definition at line 342 of file BCEngineMCMC.h.
{ fMCMCNIterationsUpdateMax = n; }
void BCEngineMCMC::MCMCSetNLag | ( | int | n | ) | [inline] |
void BCEngineMCMC::MCMCSetPrecision | ( | BCEngineMCMC::Precision | precision | ) |
Set the precision for the MCMC run.
Definition at line 105 of file BCEngineMCMC.cxx.
{ switch(precision) { case BCEngineMCMC::kLow: fMCMCNChains = 1; fMCMCNLag = 1; fMCMCNIterationsMax = 10000; fMCMCNIterationsRun = 10000; fMCMCNIterationsPreRunMin = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCRValue = 100; break; case BCEngineMCMC::kMedium: fMCMCNChains = 5; fMCMCNLag = 1; fMCMCNIterationsMax = 100000; fMCMCNIterationsRun = 100000; fMCMCNIterationsPreRunMin = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCRValue = 100; break; case BCEngineMCMC::kHigh: fMCMCNChains = 10; fMCMCNLag = 10; fMCMCNIterationsMax = 1000000; fMCMCNIterationsRun = 1000000; fMCMCNIterationsPreRunMin = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCRValue = 100; break; case BCEngineMCMC::kVeryHigh: fMCMCNChains = 10; fMCMCNLag = 10; fMCMCNIterationsMax = 10000000; fMCMCNIterationsRun = 10000000; fMCMCNIterationsPreRunMin = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCRValue = 100; break; } // re-initialize MCMCInitialize(); }
void BCEngineMCMC::MCMCSetRValueCriterion | ( | double | r | ) | [inline] |
Sets the R-value criterion for convergence of all chains.
Definition at line 398 of file BCEngineMCMC.h.
{ fMCMCRValueCriterion = r; }
void BCEngineMCMC::MCMCSetRValueParametersCriterion | ( | double | r | ) | [inline] |
Sets the parameter R-value criterion for convergence of all chains
Definition at line 403 of file BCEngineMCMC.h.
{ fMCMCRValueParametersCriterion = r; }
void BCEngineMCMC::MCMCSetRValueStrict | ( | bool | strict = true |
) | [inline] |
Use strict or relaxed rule for Gelman/Rubin R-value
Definition at line 407 of file BCEngineMCMC.h.
{ fMCMCRValueUseStrict = strict; }
void BCEngineMCMC::MCMCSetTrialFunctionScaleFactor | ( | std::vector< double > | scale | ) | [inline] |
Set the scale factors for the trial functions
scale | a vector of doubles containing the scale factors |
Definition at line 303 of file BCEngineMCMC.h.
{ fMCMCTrialFunctionScaleFactorStart = scale; }
void BCEngineMCMC::MCMCSetValuesDefault | ( | ) |
Set the default values for the MCMC chain.
Definition at line 44 of file BCEngineMCMC.cxx.
{ fMCMCNParameters = 0; fMCMCFlagWriteChainToFile = false; fMCMCFlagWritePreRunToFile = false; fMCMCFlagPreRun = true; fMCMCFlagRun = false; fMCMCFlagFillHistograms = true; fMCMCEfficiencyMin = 0.15; fMCMCEfficiencyMax = 0.50; fMCMCFlagInitialPosition = 1; fMCMCNLag = 1; fMCMCCurrentIteration = -1; fMCMCCurrentChain = -1; MCMCSetValuesDetail(); }
void BCEngineMCMC::MCMCSetValuesDetail | ( | ) |
Set the values for a detailed MCMC run.
Definition at line 84 of file BCEngineMCMC.cxx.
{ fMCMCNChains = 5; fMCMCNIterationsMax = 1000000; fMCMCNIterationsRun = 100000; fMCMCNIterationsPreRunMin = 500; fMCMCFlagInitialPosition = 1; fMCMCRValueUseStrict = false; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCNIterationsConvergenceGlobal = -1; fMCMCFlagConvergenceGlobal = false; fMCMCRValue = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCFlagOrderParameters = true; fMCMCCurrentIteration = -1; fMCMCCurrentChain = -1; }
void BCEngineMCMC::MCMCSetValuesQuick | ( | ) |
Set the values for a quick MCMC run.
Definition at line 63 of file BCEngineMCMC.cxx.
{ fMCMCNChains = 2; fMCMCNIterationsMax = 10000; fMCMCNIterationsRun = 10000; fMCMCNIterationsPreRunMin = 500; fMCMCFlagInitialPosition = 1; fMCMCRValueUseStrict = false; fMCMCRValueCriterion = 0.1; fMCMCRValueParametersCriterion = 0.1; fMCMCNIterationsConvergenceGlobal = -1; fMCMCFlagConvergenceGlobal = false; fMCMCRValue = 100; fMCMCNIterationsUpdate = 1000; fMCMCNIterationsUpdateMax = 10000; fMCMCFlagOrderParameters = true; fMCMCCurrentIteration = -1; fMCMCCurrentChain = -1; }
void BCEngineMCMC::MCMCSetWriteChainToFile | ( | bool | flag | ) | [inline] |
Sets flag to write Markov chains to file.
Definition at line 357 of file BCEngineMCMC.h.
{ fMCMCFlagWriteChainToFile = flag; }
void BCEngineMCMC::MCMCSetWritePreRunToFile | ( | bool | flag | ) | [inline] |
Sets flag to write pre run to file.
Definition at line 362 of file BCEngineMCMC.h.
{ fMCMCFlagWritePreRunToFile = flag; }
void BCEngineMCMC::MCMCTrialFunction | ( | int | ichain, | |
std::vector< double > & | x | |||
) | [virtual] |
Random walk trial function. The default trial function is a Breit-Wigner. It can be overloaded by the user to set the trial function.
ichain | the chain index | |
x | point with the dimension fMCMCNParameters |
Definition at line 489 of file BCEngineMCMC.cxx.
{ // call MCMCTrialFunctionSingle() for all parameters by default for (int i = 0; i < fMCMCNParameters; ++i) x[i] = MCMCTrialFunctionSingle(ichain, i); }
double BCEngineMCMC::MCMCTrialFunctionSingle | ( | int | ichain, | |
int | ipar | |||
) | [virtual] |
Random walk trial function. The default trial function is a Breit-Wigner. It can be overloaded by the user to set the trial function.
ichain | the chain index | |
ipar | the parameter index |
Definition at line 497 of file BCEngineMCMC.cxx.
{ // no check of range for performance reasons // use uniform distribution // return = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fRandom->Rndm()); // Breit-Wigner width adjustable width return fRandom->BreitWigner(0.0, fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter]); }
BCEngineMCMC & BCEngineMCMC::operator= | ( | const BCEngineMCMC & | engineMCMC | ) |
Defaut assignment operator
Definition at line 259 of file BCEngineMCMC.cxx.
{ fMCMCPointerToGetProposalPoint = enginemcmc.fMCMCPointerToGetProposalPoint; fMCMCNParameters = enginemcmc.fMCMCNParameters; fMCMCBoundaryMin = enginemcmc.fMCMCBoundaryMin; fMCMCBoundaryMax = enginemcmc.fMCMCBoundaryMax; fMCMCFlagsFillHistograms = enginemcmc.fMCMCFlagsFillHistograms; fMCMCNChains = enginemcmc.fMCMCNChains; fMCMCNLag = enginemcmc.fMCMCNLag; fMCMCNIterations = enginemcmc.fMCMCNIterations; fMCMCCurrentIteration = enginemcmc.fMCMCCurrentIteration; fMCMCCurrentChain = enginemcmc.fMCMCCurrentChain; fMCMCNIterationsUpdate = enginemcmc.fMCMCNIterationsUpdate; fMCMCNIterationsUpdateMax = enginemcmc.fMCMCNIterationsUpdateMax; fMCMCNIterationsConvergenceGlobal = enginemcmc.fMCMCNIterationsConvergenceGlobal; fMCMCFlagConvergenceGlobal = enginemcmc.fMCMCFlagConvergenceGlobal; fMCMCNIterationsMax = enginemcmc.fMCMCNIterationsMax; fMCMCNIterationsRun = enginemcmc.fMCMCNIterationsRun; fMCMCNIterationsPreRunMin = enginemcmc.fMCMCNIterationsPreRunMin; fMCMCNTrialsTrue = enginemcmc.fMCMCNTrialsTrue; fMCMCNTrialsFalse = enginemcmc.fMCMCNTrialsFalse; fMCMCFlagWriteChainToFile = enginemcmc.fMCMCFlagWriteChainToFile; fMCMCFlagWritePreRunToFile = enginemcmc.fMCMCFlagWritePreRunToFile; fMCMCTrialFunctionScaleFactor = enginemcmc.fMCMCTrialFunctionScaleFactor; fMCMCTrialFunctionScaleFactorStart = enginemcmc.fMCMCTrialFunctionScaleFactorStart; fMCMCFlagPreRun = enginemcmc.fMCMCFlagPreRun; fMCMCFlagRun = enginemcmc.fMCMCFlagRun; fMCMCInitialPosition = enginemcmc.fMCMCInitialPosition; fMCMCEfficiencyMin = enginemcmc.fMCMCEfficiencyMin; fMCMCEfficiencyMax = enginemcmc.fMCMCEfficiencyMax; fMCMCFlagInitialPosition = enginemcmc.fMCMCFlagInitialPosition; fMCMCFlagOrderParameters = enginemcmc.fMCMCFlagOrderParameters; fMCMCFlagFillHistograms = enginemcmc.fMCMCFlagFillHistograms; fMCMCPhase = enginemcmc.fMCMCPhase; fMCMCCycle = enginemcmc.fMCMCCycle; fMCMCx = enginemcmc.fMCMCx; fMCMCxMax = enginemcmc.fMCMCxMax; fMCMCxMean = enginemcmc.fMCMCxMean; fMCMCxVar = enginemcmc.fMCMCxVar; fMCMCxLocal = enginemcmc.fMCMCxLocal; fMCMCprob = enginemcmc.fMCMCprob; fMCMCprobMax = enginemcmc.fMCMCprobMax; fMCMCprobMean = enginemcmc.fMCMCprobMean; fMCMCprobVar = enginemcmc.fMCMCprobVar; fMCMCRValueUseStrict = enginemcmc.fMCMCRValueUseStrict; fMCMCRValueCriterion = enginemcmc.fMCMCRValueCriterion ; fMCMCRValueParametersCriterion = enginemcmc.fMCMCRValueParametersCriterion; fMCMCRValue = enginemcmc.fMCMCRValue; fMCMCRValueParameters = enginemcmc.fMCMCRValueParameters; if (enginemcmc.fRandom) fRandom = new TRandom3(*enginemcmc.fRandom); else fRandom = 0; fMCMCH1NBins = enginemcmc.fMCMCH1NBins; for (int i = 0; i < int(enginemcmc.fMCMCH1Marginalized.size()); ++i) { if (enginemcmc.fMCMCH1Marginalized.at(i)) fMCMCH1Marginalized.push_back(new TH1D(*(enginemcmc.fMCMCH1Marginalized.at(i)))); else fMCMCH1Marginalized.push_back(0); } for (int i = 0; i < int(enginemcmc.fMCMCH2Marginalized.size()); ++i) { if (enginemcmc.fMCMCH2Marginalized.at(i)) fMCMCH2Marginalized.push_back(new TH2D(*(enginemcmc.fMCMCH2Marginalized.at(i)))); else fMCMCH2Marginalized.push_back(0); } for (int i = 0; i < int(enginemcmc.fMCMCTrees.size()); ++i) { // debugKK fMCMCTrees.push_back(0); } // return this return *this; }
int BCEngineMCMC::SetMarginalized | ( | int | index, | |
TH1D * | h | |||
) |
Sets the histogram with 1D marginalized distributions for parameter.
i | index of the parameter | |
h | pointer to an existing histogram |
Definition at line 1653 of file BCEngineMCMC.cxx.
{ if((int)fMCMCH1Marginalized.size()<=index) return 0; if(h==0) return 0; if((int)fMCMCH1Marginalized.size()==index) fMCMCH1Marginalized.push_back(h); else fMCMCH1Marginalized[index]=h; return index; }
int BCEngineMCMC::SetMarginalized | ( | int | index1, | |
int | index2, | |||
TH2D * | h | |||
) |
Sets the histogram with 2D marginalized distributions for two parameters.
index1 | index of the first parameter | |
index2 | index of the second parameter | |
h | pointer to an existing histogram |
Definition at line 1670 of file BCEngineMCMC.cxx.
{ int counter = 0; int index = 0; // search for correct combination for(int i = 0; i < fMCMCNParameters; i++) for (int j = 0; j < i; ++j) { if(j == index1 && i == index2) index = counter; counter++; } if((int)fMCMCH2Marginalized.size()<=index) return 0; if(h==0) return 0; if((int)fMCMCH2Marginalized.size()==index) fMCMCH2Marginalized.push_back(h); else fMCMCH2Marginalized[index]=h; return index; }
std::vector<double> BCEngineMCMC::fMCMCBoundaryMax [protected] |
Definition at line 594 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCBoundaryMin [protected] |
Parameter boundaries
Definition at line 593 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCCurrentChain [protected] |
The current chain index. If not called within the running of the algorithm, return -1.
Definition at line 621 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCCurrentIteration [protected] |
The current iteration number. If not called within the running of the algorithm, return -1.
Definition at line 616 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCCycle [protected] |
The cycle of the pre-run
Definition at line 727 of file BCEngineMCMC.h.
double BCEngineMCMC::fMCMCEfficiencyMax [protected] |
The maximum allowed efficiency for MCMC
Definition at line 701 of file BCEngineMCMC.h.
double BCEngineMCMC::fMCMCEfficiencyMin [protected] |
The minimum required efficiency for MCMC
Definition at line 697 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagConvergenceGlobal [protected] |
Flag for convergence
Definition at line 638 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagFillHistograms [protected] |
Flag which controls fill histograms during main run.
Definition at line 716 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCFlagInitialPosition [protected] |
Variable which defines the initial position. 0 (default) center of the allowed region, (1) random initial position (2) pre-defined intial position.
Definition at line 707 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagOrderParameters [protected] |
Flag which controls the sequence parameters during the running of the MCMC.
Definition at line 712 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagPreRun [protected] |
Defines if a prerun should be performed or not
Definition at line 682 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagRun [protected] |
Defines if MCMC has been performed or not
Definition at line 686 of file BCEngineMCMC.h.
std::vector<bool> BCEngineMCMC::fMCMCFlagsFillHistograms [protected] |
Parameter flags for marginalization
Definition at line 598 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagWriteChainToFile [protected] |
Flag to write Markov chains to file
Definition at line 664 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCFlagWritePreRunToFile [protected] |
Flag to write pre run to file
Definition at line 668 of file BCEngineMCMC.h.
std::vector<TH1D *> BCEngineMCMC::fMCMCH1Marginalized [protected] |
An array of marginalized distributions
Definition at line 809 of file BCEngineMCMC.h.
std::vector<int> BCEngineMCMC::fMCMCH1NBins [protected] |
Number of bins per dimension for the marginalized distributions.
Definition at line 805 of file BCEngineMCMC.h.
std::vector<TH2D *> BCEngineMCMC::fMCMCH2Marginalized [protected] |
Definition at line 810 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCInitialPosition [protected] |
The intial position of each Markov chain. The length of the vectors is equal to fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on
Definition at line 693 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNChains [protected] |
Number of Markov chains ran in parallel
Definition at line 602 of file BCEngineMCMC.h.
std::vector<int> BCEngineMCMC::fMCMCNIterations [protected] |
Number of total iterations of the Markov chains. The length of the vector is equal to fMCMCNChains.
Definition at line 611 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsConvergenceGlobal [protected] |
Number of iterations needed for all chains to convergence simulaneously
Definition at line 634 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsMax [protected] |
Maximum number of iterations for a Markov chain prerun
Definition at line 642 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsPreRunMin [protected] |
Minimum number of iterations for the pre-run
Definition at line 650 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsRun [protected] |
Number of iterations for a Markov chain run
Definition at line 646 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsUpdate [protected] |
Number of iterations for updating scale factors
Definition at line 625 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNIterationsUpdateMax [protected] |
Maximum number of iterations for updating scale factors
Definition at line 629 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNLag [protected] |
The lag for the Markov Chain
Definition at line 606 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCNParameters [protected] |
Number of parameters
Definition at line 589 of file BCEngineMCMC.h.
std::vector<int> BCEngineMCMC::fMCMCNTrialsFalse [protected] |
Number of not accepted trials for each chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.
Definition at line 660 of file BCEngineMCMC.h.
std::vector<int> BCEngineMCMC::fMCMCNTrialsTrue [protected] |
Number of accepted trials for each chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.
Definition at line 655 of file BCEngineMCMC.h.
int BCEngineMCMC::fMCMCPhase [protected] |
The phase of the run. 1: pre-run, 2: main run.
Definition at line 722 of file BCEngineMCMC.h.
Pointer to a member function
Definition at line 583 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCprob [protected] |
The log of the probability of the current points of each Markov chain. The length of the vectors is fMCMCNChains.
Definition at line 759 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCprobMax [protected] |
The maximum (log) probability of each Markov chain. The length of the vector is fMCMCNChains.
Definition at line 764 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCprobMean [protected] |
The mean of all log prob values of each Markov chain. The length of the vector is equal to fMCMCNChains.
Definition at line 769 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCprobVar [protected] |
The variance of all log prob values of each Markov chain. The length of the vector is equal to fMCMCNChains.
Definition at line 774 of file BCEngineMCMC.h.
double BCEngineMCMC::fMCMCRValue [protected] |
The R-value at which the chains did converge
Definition at line 794 of file BCEngineMCMC.h.
double BCEngineMCMC::fMCMCRValueCriterion [protected] |
The R-value criterion for convergence of log-likelihood
Definition at line 786 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCRValueParameters [protected] |
The R-values for each parameter
Definition at line 797 of file BCEngineMCMC.h.
double BCEngineMCMC::fMCMCRValueParametersCriterion [protected] |
The R-value criterion for convergence of parameters
Definition at line 790 of file BCEngineMCMC.h.
bool BCEngineMCMC::fMCMCRValueUseStrict [protected] |
flag: use exactly the R-value definition of Gelman/Rubin (R_strict) or a relaxed, simplified version (R_relax) [default]. Note that R_relax <= R_strict, and in some cases even R_relax < 1, but we always have R_strict >= 1.
Definition at line 782 of file BCEngineMCMC.h.
std::vector<TTree *> BCEngineMCMC::fMCMCTrees [protected] |
The trees containing the Markov chains. The length of the vector is fMCMCNChains.
Definition at line 815 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCTrialFunctionScaleFactor [protected] |
Scales the width of the trial functions by a scale factor for each parameter and chain
Definition at line 673 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCTrialFunctionScaleFactorStart [protected] |
Start values of the scale factors for the trial functions.
Definition at line 678 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCx [protected] |
The current points of each Markov chain. The length of the vectors is equal to fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on.
Definition at line 734 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCxLocal [protected] |
A temporary vector for a single Markov chain
Definition at line 754 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCxMax [protected] |
The maximum points of each Markov chain. The length of the vector is fMCMCNChains * fMCMCNParameters. First, the values of the first Markov chain are saved, then those of the second and so on.
Definition at line 740 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCxMean [protected] |
The mean of all parameters of each Markov chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.
Definition at line 745 of file BCEngineMCMC.h.
std::vector<double> BCEngineMCMC::fMCMCxVar [protected] |
The variance of all parameters of each Markov chain. The length of the vector is equal to fMCMCNChains * fMCMCNParameters.
Definition at line 750 of file BCEngineMCMC.h.
TRandom3* BCEngineMCMC::fRandom [protected] |
Random number generator
Definition at line 801 of file BCEngineMCMC.h.