BCEngineMCMC.cxx

Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 2008, Daniel Kollar and Kevin Kroeninger.
00003  * All rights reserved.
00004  *
00005  * For the licensing terms see doc/COPYING.
00006  */
00007 
00008 // ---------------------------------------------------------
00009 
00010 #include "BAT/BCEngineMCMC.h"
00011 
00012 #include "BAT/BCParameter.h"
00013 #include "BAT/BCDataPoint.h"
00014 #include "BAT/BCLog.h"
00015 
00016 #include <TH1D.h>
00017 #include <TH2D.h>
00018 #include <TTree.h>
00019 #include <TRandom3.h>
00020 
00021 #include <iostream>
00022 #include <math.h>
00023 
00024 // ---------------------------------------------------------
00025 BCEngineMCMC::BCEngineMCMC()
00026 {
00027    // set default parameters for the mcmc
00028    this -> MCMCSetValuesDefault();
00029 
00030    // initialize random number generator
00031    fMCMCRandom = new TRandom3(0);
00032 }
00033 
00034 // ---------------------------------------------------------
00035 BCEngineMCMC::BCEngineMCMC(int n)
00036 {
00037    // set number of chains to n
00038    fMCMCNChains = n;
00039 
00040    // call default constructor
00041    BCEngineMCMC();
00042 }
00043 
00044 // ---------------------------------------------------------
00045 void BCEngineMCMC::MCMCSetValuesDefault()
00046 {
00047    fMCMCNParameters          = 0;
00048    fMCMCFlagWriteChainToFile = false;
00049    fMCMCFlagPreRun           = false;
00050    fMCMCFlagRun              = false;
00051    fMCMCFlagFillHistograms   = true;
00052    fMCMCEfficiencyMin        = 0.15;
00053    fMCMCEfficiencyMax        = 0.50;
00054    fMCMCFlagInitialPosition  = 1;
00055    fMCMCNLag                 = 1;
00056    fMCMCNIterationsUpdateMax = 10000;
00057 
00058    this -> MCMCSetValuesDetail();
00059 }
00060 
00061 // ---------------------------------------------------------
00062 void BCEngineMCMC::MCMCSetValuesQuick()
00063 {
00064    fMCMCNChains              = 1;
00065    fMCMCNIterationsMax       = 1000;
00066    fMCMCNIterationsRun       = 10000;
00067    fMCMCNIterationsBurnIn    = 0;
00068    fMCMCNIterationsPreRunMin = 0;
00069    fMCMCFlagInitialPosition  = 1;
00070    fMCMCRValueCriterion      = 0.1;
00071    fMCMCRValueParametersCriterion = 0.1;
00072    fMCMCNIterationsConvergenceGlobal = -1;
00073    fMCMCFlagConvergenceGlobal = false;
00074    fMCMCRValue               = 100;
00075    fMCMCNIterationsUpdate    = 1000;
00076    fMCMCFlagOrderParameters  = true;
00077 }
00078 
00079 // ---------------------------------------------------------
00080 void BCEngineMCMC::MCMCSetValuesDetail()
00081 {
00082    fMCMCNChains              = 5;
00083    fMCMCNIterationsMax       = 1000000;
00084    fMCMCNIterationsRun       = 100000;
00085    fMCMCNIterationsBurnIn    = 0;
00086    fMCMCNIterationsPreRunMin = 500;
00087    fMCMCFlagInitialPosition  = 1;
00088    fMCMCRValueCriterion      = 0.1;
00089    fMCMCRValueParametersCriterion = 0.1;
00090    fMCMCNIterationsConvergenceGlobal = -1;
00091    fMCMCFlagConvergenceGlobal = false;
00092    fMCMCRValue               = 100;
00093    fMCMCNIterationsUpdate    = 1000;
00094    fMCMCFlagOrderParameters  = true;
00095 }
00096 
00097 // ---------------------------------------------------------
00098 BCEngineMCMC::~BCEngineMCMC()
00099 {
00100    // delete random number generator
00101    if (fMCMCRandom)
00102       delete fMCMCRandom;
00103 
00104    // delete 1-d marginalized distributions
00105    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
00106       if (fMCMCH1Marginalized[i])
00107          delete fMCMCH1Marginalized[i];
00108    fMCMCH1Marginalized.clear();
00109 
00110    // delete 2-d marginalized distributions
00111    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
00112       if (fMCMCH2Marginalized[i])
00113          delete fMCMCH2Marginalized[i];
00114    fMCMCH2Marginalized.clear();
00115 }
00116 
00117 // ---------------------------------------------------------
00118 BCEngineMCMC::BCEngineMCMC(const BCEngineMCMC & enginemcmc)
00119 {
00120    enginemcmc.Copy(*this);
00121 }
00122 
00123 // ---------------------------------------------------------
00124 BCEngineMCMC & BCEngineMCMC::operator = (const BCEngineMCMC & enginemcmc)
00125 {
00126    if (this != &enginemcmc)
00127       enginemcmc.Copy(* this);
00128 
00129    return * this;
00130 }
00131 
00132 // --------------------------------------------------------
00133 TH1D * BCEngineMCMC::MCMCGetH1Marginalized(int index)
00134 {
00135       // check index
00136    if (index < 0 || index >= int(fMCMCH1Marginalized.size()))
00137    {
00138       BCLog::OutWarning("BCEngineMCMC::MCMCGetH1Marginalized. Index out of range.");
00139       return 0;
00140    }
00141 
00142    return fMCMCH1Marginalized[index];
00143 }
00144 
00145 // --------------------------------------------------------
00146 TH2D * BCEngineMCMC::MCMCGetH2Marginalized(int index1, int index2)
00147 {
00148    int counter = 0;
00149    int index = 0;
00150 
00151    // search for correct combination
00152    for(int i = 0; i < fMCMCNParameters; i++)
00153       for (int j = 0; j < i; ++j)
00154       {
00155          if(j == index1 && i == index2)
00156             index = counter;
00157          counter++;
00158       }
00159 
00160    // check index
00161    if (index < 0 || index >= int(fMCMCH2Marginalized.size()))
00162    {
00163       BCLog::OutWarning("BCEngineMCMC::MCMCGetH2Marginalized. Index out of range.");
00164       return 0;
00165    }
00166 
00167    return fMCMCH2Marginalized[index];
00168 }
00169 
00170 // --------------------------------------------------------
00171 std::vector <double> BCEngineMCMC::MCMCGetMaximumPoint(int i)
00172 {
00173    // create a new vector with the lenght of fMCMCNParameters
00174    std::vector <double> x;
00175 
00176    // check if i is in range
00177    if (i < 0 || i >= fMCMCNChains)
00178       return x;
00179 
00180    // copy the point in the ith chain into the temporary vector
00181    for (int j = 0; j < fMCMCNParameters; ++j)
00182    x.push_back(fMCMCMaximumPoints.at(i * fMCMCNParameters + j));
00183 
00184    return x;
00185 }
00186 
00187 // --------------------------------------------------------
00188 void BCEngineMCMC::MCMCSetNChains(int n)
00189 {
00190    fMCMCNChains = n;
00191 
00192    // re-initialize
00193    this -> MCMCInitialize();
00194 }
00195 
00196 // --------------------------------------------------------
00197 void BCEngineMCMC::MCMCSetInitialPositions(std::vector<double> x0s)
00198 {
00199    // clear the existing initial position vector
00200    fMCMCInitialPosition.clear();
00201 
00202    // copy the initial positions
00203    int n = int(x0s.size());
00204 
00205    for (int i = 0; i < n; ++i)
00206       fMCMCInitialPosition.push_back(x0s.at(i));
00207 
00208    // use these intial positions for the Markov chain
00209    this -> MCMCSetFlagInitialPosition(2);
00210 }
00211 
00212 // --------------------------------------------------------
00213 void BCEngineMCMC::MCMCSetInitialPositions(std::vector< std::vector<double> > x0s)
00214 {
00215    // create new vector
00216    std::vector <double> y0s;
00217 
00218    // loop over vector elements
00219    for (int i = 0; i < int(x0s.size()); ++i)
00220       for (int j = 0; j < int((x0s.at(i)).size()); ++j)
00221          y0s.push_back((x0s.at(i)).at(j));
00222 
00223    this -> MCMCSetInitialPositions(y0s);
00224 }
00225 
00226 // --------------------------------------------------------
00227 void BCEngineMCMC::MCMCSetFlagFillHistograms(bool flag)
00228 {
00229    fMCMCFlagFillHistograms = flag;
00230 
00231    for (int i = 0; i < fMCMCNParameters; ++i)
00232       fMCMCFlagsFillHistograms[i] = flag;
00233 }
00234 
00235 // --------------------------------------------------------
00236 void BCEngineMCMC::MCMCSetFlagFillHistograms(int index, bool flag)
00237 {
00238    // check if index is within range
00239    if (index < 0 || index > fMCMCNParameters)
00240    {
00241       BCLog::OutWarning("BCEngineMCMC :MCMCSetFlagFillHistograms. Index out of range.");
00242       return;
00243    }
00244 
00245    // set flag
00246    fMCMCFlagsFillHistograms[index] = flag;
00247 }
00248 
00249 // --------------------------------------------------------
00250 void BCEngineMCMC::MCMCSetMarkovChainTrees(std::vector <TTree *> trees)
00251 {
00252 // clear vector
00253    fMCMCTrees.clear();
00254 
00255    // copy tree
00256    for (int i = 0; i < int(trees.size()); ++i)
00257       fMCMCTrees.push_back(trees[i]);
00258 }
00259 
00260 // --------------------------------------------------------
00261 void BCEngineMCMC::Copy(BCEngineMCMC & enginemcmc) const
00262 {}
00263 
00264 // --------------------------------------------------------
00265 void BCEngineMCMC::MCMCTrialFunction(int chain, std::vector <double> &x)
00266 {
00267    // use uniform distribution for now
00268 // for (int i = 0; i < fMCMCNParameters; ++i)
00269 //    x[i] = fMCMCTrialFunctionScaleFactor[i] * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00270 
00271    // Breit-Wigner with adjustable width
00272    for (int i = 0; i < fMCMCNParameters; ++i)
00273       x[i] = fMCMCRandom -> BreitWigner(0.0, fMCMCTrialFunctionScaleFactor[chain * fMCMCNParameters + i]);
00274 }
00275 
00276 // --------------------------------------------------------
00277 void BCEngineMCMC::MCMCTrialFunctionSingle(int ichain, int iparameter, std::vector <double> &x)
00278 {
00279    // no check of range for performance reasons
00280 
00281    // use uniform distribution
00282 // x[iparameter] = fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter] * 2.0 * (0.5 - fMCMCRandom -> Rndm());
00283 
00284    // Breit-Wigner width adjustable width
00285    x[iparameter] = fMCMCRandom -> BreitWigner(0.0, fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter]);
00286 }
00287 
00288 // --------------------------------------------------------
00289 std::vector <double> BCEngineMCMC::MCMCGetTrialFunctionScaleFactor(int ichain)
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 scale factors into the temporary vector
00299    for (int j = 0; j < fMCMCNParameters; ++j)
00300       x.push_back(fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNParameters + j));
00301 
00302    return x;
00303 }
00304 
00305 // --------------------------------------------------------
00306 double BCEngineMCMC::MCMCGetTrialFunctionScaleFactor(int ichain, int ipar)
00307 {
00308    // check if ichain is in range
00309    if (ichain < 0 || ichain >= fMCMCNChains)
00310       return 0;
00311 
00312    // check if ipar is in range
00313    if (ipar < 0 || ipar >= fMCMCNParameters)
00314       return 0;
00315 
00316    // return component of ipar point in the ichain chain
00317    return fMCMCTrialFunctionScaleFactor.at(ichain *  fMCMCNChains + ipar);
00318 }
00319 
00320 // --------------------------------------------------------
00321 std::vector <double> BCEngineMCMC::MCMCGetx(int ichain)
00322 {
00323    // create a new vector with the length of fMCMCNParameters
00324    std::vector <double> x;
00325 
00326    // check if ichain is in range
00327    if (ichain < 0 || ichain >= fMCMCNChains)
00328       return x;
00329 
00330    // copy the point in the ichain chain into the temporary vector
00331    for (int j = 0; j < fMCMCNParameters; ++j)
00332       x.push_back(fMCMCx.at(ichain * fMCMCNParameters + j));
00333 
00334    return x;
00335 }
00336 
00337 // --------------------------------------------------------
00338 double BCEngineMCMC::MCMCGetx(int ichain, int ipar)
00339 {
00340    // check if ichain is in range
00341    if (ichain < 0 || ichain >= fMCMCNChains)
00342       return 0;
00343 
00344    // check if ipar is in range
00345    if (ipar < 0 || ipar >= fMCMCNParameters)
00346       return 0;
00347 
00348    // return component of jth point in the ith chain
00349    return fMCMCx.at(ichain *  fMCMCNParameters + ipar);
00350 }
00351 
00352 // --------------------------------------------------------
00353 double BCEngineMCMC::MCMCGetLogProbx(int ichain)
00354 {
00355    // check if ichain is in range
00356    if (ichain < 0 || ichain >= fMCMCNChains)
00357       return -1;
00358 
00359    // return log of the probability at the current point in the ith chain
00360    return fMCMCLogProbx.at(ichain);
00361 }
00362 
00363 // --------------------------------------------------------
00364 bool BCEngineMCMC::MCMCGetProposalPointMetropolis(int chain, std::vector <double> &x)
00365 {
00366    // get unscaled random point. this point might not be in the correct volume.
00367    this -> MCMCTrialFunction(chain, x);
00368 
00369    // get a proposal point from the trial function and scale it
00370    for (int i = 0; i < fMCMCNParameters; ++i)
00371       x[i] = fMCMCx[chain * fMCMCNParameters + i] + x[i] * (fMCMCBoundaryMax.at(i) - fMCMCBoundaryMin.at(i));
00372 
00373    // check if the point is in the correct volume.
00374    for (int i = 0; i < fMCMCNParameters; ++i)
00375       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00376          return false;
00377 
00378    return true;
00379 }
00380 
00381 // --------------------------------------------------------
00382 bool BCEngineMCMC::MCMCGetProposalPointMetropolis(int chain, int parameter, std::vector <double> &x)
00383 {
00384    // get unscaled random point in the dimension of the chosen
00385    // parameter. this point might not be in the correct volume.
00386    this -> MCMCTrialFunctionSingle(chain, parameter, x);
00387 
00388    // copy the old point into the new
00389    for (int i = 0; i < fMCMCNParameters; ++i)
00390       if (i != parameter)
00391          x[i] = fMCMCx[chain * fMCMCNParameters + i];
00392 
00393    // modify the parameter under study
00394    x[parameter] = fMCMCx[chain * fMCMCNParameters + parameter] + x[parameter] * (fMCMCBoundaryMax.at(parameter) - fMCMCBoundaryMin.at(parameter));
00395 
00396    // check if the point is in the correct volume.
00397    for (int i = 0; i < fMCMCNParameters; ++i)
00398       if ((x[i] < fMCMCBoundaryMin[i]) || (x[i] > fMCMCBoundaryMax[i]))
00399          return false;
00400 
00401    return true;
00402 }
00403 
00404 // --------------------------------------------------------
00405 bool BCEngineMCMC::MCMCGetNewPointMetropolis(int chain, int parameter)
00406 {
00407    // calculate index
00408    int index = chain * fMCMCNParameters;
00409 
00410    // increase counter
00411    fMCMCNIterations[chain]++;
00412 
00413    // get proposal point
00414    if (!this -> MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal))
00415    {
00416       // increase counter
00417       fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00418 
00419       return false;
00420    }
00421 
00422    // calculate probabilities of the old and new points
00423    double p0 = fMCMCLogProbx[chain];
00424    double p1 = this -> LogEval(fMCMCxLocal);
00425 
00426    // flag for accept
00427    bool accept = false;
00428 
00429    // if the new point is more probable, keep it ...
00430    if (p1 >= p0)
00431       accept = true;
00432    // ... or else throw dice.
00433    else
00434    {
00435       double r = log(fMCMCRandom -> Rndm());
00436 
00437       if(r < p1 - p0)
00438          accept = true;
00439    }
00440 
00441    // fill the new point
00442    if(accept)
00443    {
00444       // increase counter
00445       fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++;
00446 
00447       // copy the point
00448       for(int i = 0; i < fMCMCNParameters; ++i)
00449       {
00450          // save the point
00451          fMCMCx[index + i] = fMCMCxLocal[i];
00452 
00453          // save the probability of the point
00454          fMCMCLogProbx[chain] = p1;
00455       }
00456    }
00457    else
00458    {
00459       // increase counter
00460       fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00461    }
00462 
00463    return accept;
00464 }
00465 
00466 // --------------------------------------------------------
00467 bool BCEngineMCMC::MCMCGetNewPointMetropolis(int chain)
00468 {
00469    // calculate index
00470    int index = chain * fMCMCNParameters;
00471 
00472    // increase counter
00473    fMCMCNIterations[chain]++;
00474 
00475    // get proposal point
00476    if (!this -> MCMCGetProposalPointMetropolis(chain, fMCMCxLocal))
00477    {
00478       // increase counter
00479       for (int i = 0; i < fMCMCNParameters; ++i)
00480          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00481 
00482       return false;
00483    }
00484 
00485    // calculate probabilities of the old and new points
00486    double p0 = fMCMCLogProbx[chain];
00487    double p1 = this -> LogEval(fMCMCxLocal);
00488 
00489    // flag for accept
00490    bool accept = false;
00491 
00492    // if the new point is more probable, keep it ...
00493    if (p1 >= p0)
00494       accept = true;
00495 
00496    // ... or else throw dice.
00497    else
00498    {
00499       double r = log(fMCMCRandom -> Rndm());
00500 
00501       if(r < p1 - p0)
00502          accept = true;
00503    }
00504 
00505    // fill the new point
00506    if(accept)
00507    {
00508       // increase counter
00509       for (int i = 0; i < fMCMCNParameters; ++i)
00510          fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00511 
00512       // copy the point
00513       for(int i = 0; i < fMCMCNParameters; ++i)
00514       {
00515          // save the point
00516          fMCMCx[index + i] = fMCMCxLocal[i];
00517 
00518          // save the probability of the point
00519          fMCMCLogProbx[chain] = p1;
00520       }
00521    }
00522    else
00523    {
00524       // increase counter
00525       for (int i = 0; i < fMCMCNParameters; ++i)
00526          fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00527    }
00528 
00529    return accept;
00530 }
00531 
00532 // --------------------------------------------------------
00533 void BCEngineMCMC::MCMCUpdateStatisticsCheckMaximum()
00534 {
00535    // loop over all chains
00536    for (int i = 0; i < fMCMCNChains; ++i)
00537    {
00538       // check if new maximum is found or chain is at the beginning
00539       if (fMCMCLogProbx[i] > fMCMCMaximumLogProb[i] || fMCMCNIterations[i] == 1)
00540       {
00541          // copy maximum value
00542          fMCMCMaximumLogProb[i] = fMCMCLogProbx[i];
00543 
00544          // copy mode of chain
00545          for (int j = 0; j < fMCMCNParameters; ++j)
00546             fMCMCMaximumPoints[i * fMCMCNParameters + j] = fMCMCx[i * fMCMCNParameters + j];
00547       }
00548    }
00549 }
00550 
00551 // --------------------------------------------------------
00552 void BCEngineMCMC::MCMCUpdateStatisticsFillHistograms()
00553 {
00554    // check if histograms are supposed to be filled
00555    if (!fMCMCFlagFillHistograms)
00556       return;
00557 
00558    // loop over chains
00559    for (int i = 0; i < fMCMCNChains; ++i)
00560    {
00561       // fill each 1-dimensional histogram (if supposed to be filled)
00562       for (int j = 0; j < fMCMCNParameters; ++j)
00563          if (fMCMCFlagsFillHistograms.at(j))
00564             fMCMCH1Marginalized[j] -> Fill(fMCMCx[i * fMCMCNParameters + j]);
00565 
00566       // fill each 2-dimensional histogram (if supposed to be filled)
00567       int counter = 0;
00568 
00569       for (int j = 0; j < fMCMCNParameters; ++j)
00570          for (int k = 0; k < j; ++k)
00571          {
00572             if (fMCMCFlagsFillHistograms.at(j) && fMCMCFlagsFillHistograms.at(k))
00573                fMCMCH2Marginalized[counter]->Fill(fMCMCx[i*fMCMCNParameters+k],fMCMCx[i* fMCMCNParameters+j]);
00574             counter ++;
00575          }
00576    }
00577 }
00578 
00579 // --------------------------------------------------------
00580 void BCEngineMCMC::MCMCUpdateStatisticsTestConvergenceAllChains()
00581 {
00582    // calculate number of entries in this part of the chain
00583    int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0];
00584 
00585    if (fMCMCNChains > 1 && npoints > 1)
00586    {
00587       // define flag for convergence
00588       bool flag_convergence = true;
00589 
00590       // loop over parameters
00591       for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00592       {
00593          double sum = 0;
00594          double sum2 = 0;
00595          double sumv = 0;
00596 
00597          // loop over chains
00598          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00599          {
00600             // get parameter index
00601             int index = ichains * fMCMCNParameters + iparameters;
00602 
00603             // calculate mean value of each parameter in the chain for this part
00604             fMCMCMeanx[index] += (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints);
00605 
00606             // calculate variance of each chain for this part
00607             fMCMCVariancex[index] = (1.0 - 1/double(npoints)) * fMCMCVariancex[index]
00608                   + (fMCMCx[index] - fMCMCMeanx[index]) * (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints - 1);
00609 
00610             sum  += fMCMCMeanx[index];
00611             sum2 += fMCMCMeanx[index] * fMCMCMeanx[index];
00612             sumv += fMCMCVariancex[index];
00613          }
00614 
00615          // calculate r-value for each parameter
00616          double mean = sum / double(fMCMCNChains);
00617          double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints);
00618          double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);
00619 
00620          double r = 100.0;
00621 
00622          if (W > 0)
00623          {
00624             r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W);
00625             fMCMCRValueParameters[iparameters] = r;
00626          }
00627          // set flag to false if convergence criterion is not fulfilled for the parameter
00628          if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
00629             flag_convergence = false;
00630       }
00631 
00632       // convergence criterion applied on the log-likelihood
00633       double sum = 0;
00634       double sum2 = 0;
00635       double sumv = 0;
00636 
00637       // loop over chains
00638       for (int i = 0; i < fMCMCNChains; ++i)
00639       {
00640          // calculate mean value of each chain for this part
00641          fMCMCMean[i] += (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints);
00642 
00643          // calculate variance of each chain for this part
00644          if (npoints > 1)
00645             fMCMCVariance[i] = (1.0 - 1/double(npoints)) * fMCMCVariance[i]
00646                   + (fMCMCLogProbx[i] - fMCMCMean[i]) * (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints - 1);
00647 
00648          sum  += fMCMCMean[i];
00649          sum2 += fMCMCMean[i] * fMCMCMean[i]; ;
00650          sumv += fMCMCVariance[i];
00651       }
00652 
00653       // calculate r-value for log-likelihood
00654       double mean = sum / double(fMCMCNChains);
00655       double B = (sum2 / double(fMCMCNChains) - mean * mean) * double(fMCMCNChains) / double(fMCMCNChains-1) * double(npoints);
00656       double W = sumv * double(npoints) / double(npoints - 1) / double(fMCMCNChains);
00657       double r = 100.0;
00658 
00659       if (W > 0)
00660       {
00661          r = sqrt( ( (1-1/double(npoints)) * W + 1/double(npoints) * B ) / W);
00662          fMCMCRValue = r;
00663       }
00664 
00665       // set flag to false if convergence criterion is not fulfilled for the log-likelihood
00666       if (! ((fMCMCRValue-1.0) < fMCMCRValueCriterion))
00667          flag_convergence = false;
00668 
00669       // remember number of iterations needed to converge
00670       if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
00671          fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters;
00672    }
00673 }
00674 
00675 // --------------------------------------------------------
00676 void BCEngineMCMC::MCMCUpdateStatisticsWriteChains()
00677 {
00678    // loop over all chains
00679    for (int i = 0; i < fMCMCNChains; ++i)
00680       fMCMCTrees[i] -> Fill();
00681 }
00682 
00683 // --------------------------------------------------------
00684 double BCEngineMCMC::LogEval(std::vector <double> parameters)
00685 {
00686    // test function for now
00687    // this will be overloaded by the user
00688    return 0.0;
00689 }
00690 
00691 // --------------------------------------------------------
00692 int BCEngineMCMC::MCMCMetropolisPreRun()
00693 {
00694    // print on screen
00695    BCLog::OutSummary("Pre-run Metropolis MCMC...");
00696 
00697    // initialize Markov chain
00698    this -> MCMCInitialize();
00699    this -> MCMCInitializeMarkovChains();
00700 
00701    // helper variable containing number of digits in the number of parameters
00702    int ndigits = (int)log10(fMCMCNParameters) +1;
00703    if(ndigits<4)
00704       ndigits=4;
00705 
00706    // perform burn-in run
00707    if (fMCMCNIterationsBurnIn > 0)
00708       BCLog::OutDetail(Form(" --> Start burn-in run with %i iterations", fMCMCNIterationsBurnIn));
00709 
00710    // if the flag is set then run over the parameters one after the other.
00711    if (fMCMCFlagOrderParameters)
00712    {
00713       for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
00714          for (int j = 0; j < fMCMCNChains; ++j)
00715             for (int k = 0; k < fMCMCNParameters; ++k)
00716                this -> MCMCGetNewPointMetropolis(j, k);
00717    }
00718    // if the flag is not set then run over the parameters at the same time.
00719    else
00720    {
00721       for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
00722          // loop over chains
00723          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00724             // get new point
00725             this -> MCMCGetNewPointMetropolis(ichains);
00726    }
00727 
00728    // reset run statistics
00729    this -> MCMCResetRunStatistics();
00730    fMCMCNIterationsConvergenceGlobal = -1;
00731 
00732    // perform run
00733    BCLog::OutSummary(Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax));
00734 
00735 
00736    // don't write to file during pre run
00737    bool tempflag_writetofile = fMCMCFlagWriteChainToFile;
00738    fMCMCFlagWriteChainToFile = false;
00739 
00740    // initialize counter variables and flags
00741    int counter = 0;
00742    int counterupdate = 0;
00743    bool convergence = false;     // convergence reached?
00744    bool flagefficiency = false;  // efficiency reached?
00745 
00746    // array of efficiencies
00747    std::vector <double> efficiency;
00748    efficiency.assign(fMCMCNParameters * fMCMCNChains, 0.0);
00749 
00750    // how often to check convergence and efficiencies?
00751    int updateLimit = ( fMCMCNIterationsUpdateMax<fMCMCNIterationsUpdate*(fMCMCNParameters+1)  && fMCMCNIterationsUpdateMax>0 ) ?
00752          fMCMCNIterationsUpdateMax : fMCMCNIterationsUpdate*(fMCMCNParameters+1);
00753 
00754    // run chain ...
00755    // (a) for at least a minimum number of iterations,
00756    // (b) until a maximum number of iterations is reached,
00757    // (c) or until convergence is reached and the efficiency is in the
00758    //     specified region
00759    while (counter < fMCMCNIterationsPreRunMin || (counter < fMCMCNIterationsMax && !(convergence && flagefficiency)))
00760    {
00761       // set convergence to false by default
00762       convergence = false;
00763 
00764       // set number of iterations needed to converge to negative
00765       fMCMCNIterationsConvergenceGlobal = -1;
00766 
00767       // if the flag is set then run over the parameters one after the other.
00768       if (fMCMCFlagOrderParameters)
00769       {
00770          // loop over parameters
00771          for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00772          {
00773             // loop over chains
00774             for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00775                this -> MCMCGetNewPointMetropolis(ichains, iparameters);
00776 
00777             // search for global maximum
00778             this -> MCMCUpdateStatisticsCheckMaximum();
00779          } // end loop over parameters
00780       } // end if fMCMCFlagOrderParameters
00781 
00782       // if the flag is not set then run over the parameters at the same time.
00783       else
00784       {
00785          // loop over chains
00786          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00787             // get new point
00788             this -> MCMCGetNewPointMetropolis(ichains);
00789 
00790          // search for global maximum
00791          this -> MCMCUpdateStatisticsCheckMaximum();
00792       }
00793 
00794       // progress printout
00795       if ( counter > 0 && counter % fMCMCNIterationsUpdate == 0 )
00796          BCLog::OutDetail(Form(" --> Iteration %i", fMCMCNIterations[0]/fMCMCNParameters - 1));
00797 
00798       //-------------------------------------------
00799       // update scale factors and check convergence
00800       //-------------------------------------------
00801       if (counterupdate % updateLimit == 0 && counterupdate > 0 && counter >= fMCMCNIterationsPreRunMin)
00802       {
00803          bool rvalues_ok = true;
00804          static bool has_converged = false;
00805 
00806          // -----------------------------
00807          // check convergence status
00808          // -----------------------------
00809 
00810          // reset the number of iterations needed for convergence to
00811          // negative
00812          fMCMCNIterationsConvergenceGlobal = -1;
00813 
00814          this -> MCMCUpdateStatisticsTestConvergenceAllChains();
00815 
00816          // check convergence
00817          if (fMCMCNIterationsConvergenceGlobal > 0)
00818             convergence = true;
00819 
00820          // print convergence status:
00821          if (convergence && fMCMCNChains > 1)
00822             BCLog::OutDetail(Form("     * Convergence status: Set of %i Markov chains converged within %i iterations.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
00823          else if (!convergence && fMCMCNChains > 1)
00824          {
00825             BCLog::OutDetail(Form("     * Convergence status: Set of %i Markov chains did not converge after %i iterations.", fMCMCNChains, counter));
00826 
00827             BCLog::OutDetail("       - R-Values:");
00828             for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
00829             {
00830                if(fabs(fMCMCRValueParameters[iparameter]-1.) < fMCMCRValueParametersCriterion)
00831                   BCLog::OutDetail(Form( TString::Format("         parameter %%%di :  %%.06f",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
00832                else
00833                {
00834                   BCLog::OutDetail(Form( TString::Format("         parameter %%%di :  %%.06f <--",ndigits), iparameter, fMCMCRValueParameters.at(iparameter)));
00835                   rvalues_ok = false;
00836                }
00837             }
00838             if(fabs(fMCMCRValue-1.) < fMCMCRValueCriterion)
00839                BCLog::OutDetail(Form("         log-likelihood :  %.06f", fMCMCRValue));
00840             else
00841             {
00842                BCLog::OutDetail(Form("         log-likelihood :  %.06f <--", fMCMCRValue));
00843                rvalues_ok = false;
00844             }
00845          }
00846 
00847          if(!has_converged)
00848             if(rvalues_ok)
00849                has_converged = true;
00850 
00851          // -----------------------------
00852          // check efficiency status
00853          // -----------------------------
00854 
00855          // set flag
00856          flagefficiency = true;
00857 
00858          bool flagprintefficiency = true;
00859 
00860          // loop over chains
00861          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00862          {
00863             // loop over parameters
00864             for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
00865             {
00866                // global index of the parameter (throughout all the chains)
00867                int ipc = ichains * fMCMCNParameters + iparameter;
00868 
00869                // calculate efficiency
00870                efficiency[ipc] = double(fMCMCNTrialsTrue[ipc]) / double(fMCMCNTrialsTrue[ipc] + fMCMCNTrialsFalse[ipc]);
00871 
00872                // adjust scale factors if efficiency is too low
00873                if (efficiency[ipc] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[ipc] > .01)
00874                {
00875                   if (flagprintefficiency)
00876                   {
00877                      BCLog::OutDetail(Form("     * Efficiency status: Efficiencies not within pre-defined range."));
00878                      BCLog::OutDetail(Form("       - Efficiencies:"));
00879                      flagprintefficiency = false;
00880                   }
00881 
00882                   double fscale=2.;
00883                   if(has_converged && fMCMCEfficiencyMin/efficiency[ipc] > 2.)
00884                      fscale = 4.;
00885                   fMCMCTrialFunctionScaleFactor[ipc] /= fscale;
00886 
00887                   BCLog::OutDetail(Form("         Efficiency of parameter %i dropped below %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
00888                         iparameter, 100. * fMCMCEfficiencyMin, 100. * efficiency[ipc], ichains, fMCMCTrialFunctionScaleFactor[ipc]));
00889                }
00890 
00891                // adjust scale factors if efficiency is too high
00892                else if (efficiency[ipc] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[ipc] < 1.0)
00893                {
00894                   if (flagprintefficiency)
00895                   {
00896                      BCLog::OutDetail(Form("   * Efficiency status: Efficiencies not within pre-defined ranges."));
00897                      BCLog::OutDetail(Form("     - Efficiencies:"));
00898                      flagprintefficiency = false;
00899                   }
00900 
00901                   fMCMCTrialFunctionScaleFactor[ipc] *= 2.;
00902 
00903                   BCLog::OutDetail(Form("         Efficiency of parameter %i above %.2lf%% (eps = %.2lf%%) in chain %i. Set scale to %.4g",
00904                         iparameter, 100.0 * fMCMCEfficiencyMax, 100.0 * efficiency[ipc], ichains, fMCMCTrialFunctionScaleFactor[ipc]));
00905                }
00906 
00907                // check flag
00908                if ((efficiency[ipc] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[ipc] > .01)
00909                      || (efficiency[ipc] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[ipc] < 1.))
00910                   flagefficiency = false;
00911 
00912                // reset counters
00913                counterupdate = 0;
00914                fMCMCNTrialsTrue[ipc] = 0;
00915                fMCMCNTrialsFalse[ipc] = 0;
00916             } // end of running over all parameters
00917 
00918             // reset counters
00919             fMCMCMean[ichains] = 0;
00920             fMCMCVariance[ichains] = 0;
00921          } // end of running over all chains
00922 
00923          // print to scrren
00924          if (flagefficiency)
00925             BCLog::OutDetail(Form("     * Efficiency status: Efficiencies within pre-defined ranges."));
00926 
00927       } // end if update scale factors and check convergence
00928 
00929       // increase counter
00930       counter++;
00931       counterupdate++;
00932    } // end of running
00933 
00934    // did we check convergence at least once ?
00935    if (counter-1<updateLimit)
00936    {
00937       BCLog::OutWarning(" Convergence never checked !");
00938       BCLog::OutWarning("   Increase maximum number of iterations in the pre-run /MCMCSetNIterationsMax()/");
00939       BCLog::OutWarning("   or decrease maximum number of iterations for update  /MCMCSetNIterationsUpdateMax()/");
00940    }
00941 
00942    // ---------------
00943    // after chain run
00944    // ---------------
00945 
00946    // define convergence status
00947    if (fMCMCNIterationsConvergenceGlobal > 0)
00948       fMCMCFlagConvergenceGlobal = true;
00949    else
00950       fMCMCFlagConvergenceGlobal = false;
00951 
00952    // print convergence status
00953    if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && !flagefficiency)
00954       BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations but could not adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
00955 
00956    else if (fMCMCFlagConvergenceGlobal && fMCMCNChains > 1 && flagefficiency)
00957       BCLog::OutSummary(Form(" --> Set of %i Markov chains converged within %i iterations and could adjust scales.", fMCMCNChains, fMCMCNIterationsConvergenceGlobal));
00958 
00959    else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && flagefficiency)
00960       BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations.", fMCMCNChains, fMCMCNIterationsMax));
00961 
00962    else if (!fMCMCFlagConvergenceGlobal && (fMCMCNChains > 1) && !flagefficiency)
00963       BCLog::OutSummary(Form(" --> Set of %i Markov chains did not converge within %i iterations and could not adjust scales.", fMCMCNChains, fMCMCNIterationsMax));
00964 
00965    else if(fMCMCNChains == 1)
00966       BCLog::OutSummary(" --> No convergence criterion for a single chain defined.");
00967 
00968    else
00969       BCLog::OutSummary(" --> Only one Markov chain. No global convergence criterion defined.");
00970 
00971    BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", counter));
00972 
00973 
00974    // print efficiencies
00975    std::vector <double> efficiencies;
00976 
00977    for (int i = 0; i < fMCMCNParameters; ++i)
00978       efficiencies.push_back(0.);
00979 
00980    BCLog::OutDetail(" --> Average efficiencies:");
00981    for (int i = 0; i < fMCMCNParameters; ++i)
00982    {
00983       for (int j = 0; j < fMCMCNChains; ++j)
00984          efficiencies[i] += efficiency[j * fMCMCNParameters + i] / double(fMCMCNChains);
00985 
00986       BCLog::OutDetail(Form(TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * efficiencies[i]));
00987    }
00988 
00989 
00990    // print scale factors
00991    std::vector <double> scalefactors;
00992 
00993    for (int i = 0; i < fMCMCNParameters; ++i)
00994       scalefactors.push_back(0.0);
00995 
00996    BCLog::OutDetail(" --> Average scale factors:");
00997    for (int i = 0; i < fMCMCNParameters; ++i)
00998    {
00999       for (int j = 0; j < fMCMCNChains; ++j)
01000          scalefactors[i] += fMCMCTrialFunctionScaleFactor[j * fMCMCNParameters + i] / double(fMCMCNChains);
01001 
01002       BCLog::OutDetail(Form( TString::Format(" -->      parameter %%%di :  %%.02f%%%%",ndigits), i, 100. * scalefactors[i]));
01003    }
01004 
01005    // reset flag
01006    fMCMCFlagWriteChainToFile = tempflag_writetofile;
01007 
01008    // set pre-run flag
01009    fMCMCFlagPreRun = true;
01010 
01011    // no error
01012    return 1;
01013 }
01014 
01015 // --------------------------------------------------------
01016 int BCEngineMCMC::MCMCMetropolis()
01017 {
01018    // check if prerun has been performed
01019    if (!fMCMCFlagPreRun)
01020       this -> MCMCMetropolisPreRun();
01021 
01022    // print to screen
01023    BCLog::OutSummary( "Run Metropolis MCMC...");
01024 
01025    // reset run statistics
01026    this -> MCMCResetRunStatistics();
01027 
01028    // perform run
01029    BCLog::OutSummary(Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));
01030 
01031 // int counterupdate = 0;
01032 // bool convergence = false;
01033 // bool flagefficiency = false;
01034 
01035 // std::vector <double> efficiency;
01036 
01037 // for (int i = 0; i < fMCMCNParameters; ++i)
01038 //    for (int j = 0; j < fMCMCNChains; ++j)
01039 //       efficiency.push_back(0.0);
01040 
01041    int nwrite = fMCMCNIterationsRun/10;
01042    if(nwrite < 100)
01043       nwrite=100;
01044    else if(nwrite < 500)
01045       nwrite=1000;
01046    else if(nwrite < 10000)
01047       nwrite=1000;
01048    else
01049       nwrite=10000;
01050 
01051    // start the run
01052    for (int iiterations = 0; iiterations < fMCMCNIterationsRun; ++iiterations)
01053    {
01054       if ( (iiterations+1)%nwrite == 0 )
01055          BCLog::OutDetail(Form(" --> iteration number %i (%.2f%%)", iiterations+1, (double)(iiterations+1)/(double)fMCMCNIterationsRun*100.));
01056 
01057       // if the flag is set then run over the parameters one after the other.
01058       if (fMCMCFlagOrderParameters)
01059       {
01060          // loop over parameters
01061          for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
01062          {
01063             // loop over chains
01064             for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01065                this -> MCMCGetNewPointMetropolis(ichains, iparameters);
01066 
01067             // update search for maximum
01068             this -> MCMCUpdateStatisticsCheckMaximum();
01069 
01070             // check if the current iteration is consistent with the lag
01071             if ( (fMCMCNParameters * iiterations + iparameters) % (fMCMCNLag * fMCMCNParameters) == 0)
01072             {
01073                // fill histograms
01074                this -> MCMCUpdateStatisticsFillHistograms();
01075 
01076                // write chain to file
01077                if (fMCMCFlagWriteChainToFile)
01078                   this -> MCMCUpdateStatisticsWriteChains();
01079 
01080                // do anything interface
01081                this -> MCMCIterationInterface();
01082             }
01083 
01084          } // end loop over all parameters
01085       }
01086       // if the flag is not set then run over the parameters at the same time.
01087       else
01088       {
01089          // loop over chains
01090          for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01091             // get new point
01092             this -> MCMCGetNewPointMetropolis(ichains);
01093 
01094          // update search for maximum
01095          this -> MCMCUpdateStatisticsCheckMaximum();
01096 
01097          // check if the current iteration is consistent with the lag
01098          if (iiterations % fMCMCNLag == 0)
01099          {
01100             // fill histograms
01101             this -> MCMCUpdateStatisticsFillHistograms();
01102 
01103             // write chain to file
01104             if (fMCMCFlagWriteChainToFile)
01105                this -> MCMCUpdateStatisticsWriteChains();
01106 
01107             // do anything interface
01108             this -> MCMCIterationInterface();
01109          }
01110       }
01111 
01112    } // end run
01113 
01114    // print convergence status
01115    BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));
01116 
01117    // print modes
01118 
01119    // find global maximum
01120    double probmax = fMCMCMaximumLogProb.at(0);
01121    int probmaxindex = 0;
01122 
01123    // loop over all chains and find the maximum point
01124    for (int i = 1; i < fMCMCNChains; ++i)
01125       if (fMCMCMaximumLogProb.at(i) > probmax)
01126       {
01127          probmax = fMCMCMaximumLogProb.at(i);
01128          probmaxindex = i;
01129       }
01130 
01131    BCLog::OutDetail(" --> Global mode from MCMC:");
01132    int ndigits = (int) log10(fMCMCNParameters);
01133    for (int i = 0; i < fMCMCNParameters; ++i)
01134       BCLog::OutDetail(Form( TString::Format(" -->      parameter %%%di:   %%.4g", ndigits+1),
01135             i, fMCMCMaximumPoints[probmaxindex * fMCMCNParameters + i]));
01136 
01137    // set flags
01138    fMCMCFlagPreRun = false;
01139    fMCMCFlagRun = true;
01140 
01141    return 1;
01142 }
01143 
01144 // --------------------------------------------------------
01145 void BCEngineMCMC::MCMCResetRunStatistics()
01146 {
01147    for (int j = 0; j < fMCMCNChains; ++j)
01148    {
01149       fMCMCNIterations[j]  = 0;
01150       fMCMCNTrialsTrue[j]  = 0;
01151       fMCMCNTrialsFalse[j] = 0;
01152       fMCMCMean[j]         = 0;
01153       fMCMCVariance[j]     = 0;
01154 
01155       for (int k = 0; k < fMCMCNParameters; ++k)
01156       {
01157          fMCMCNTrialsTrue[j * fMCMCNParameters + k]  = 0;
01158          fMCMCNTrialsFalse[j * fMCMCNParameters + k] = 0;
01159       }
01160    }
01161 
01162    // reset marginalized distributions
01163    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01164       if (fMCMCH1Marginalized[i])
01165          fMCMCH1Marginalized[i] -> Reset();
01166 
01167    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01168       if (fMCMCH2Marginalized[i])
01169          fMCMCH2Marginalized[i] -> Reset();
01170 
01171    fMCMCRValue = 100;
01172 }
01173 
01174 // --------------------------------------------------------
01175 int BCEngineMCMC::MCMCAddParameter(double min, double max)
01176 {
01177    // add the boundaries to the corresponding vectors
01178    fMCMCBoundaryMin.push_back(min);
01179    fMCMCBoundaryMax.push_back(max);
01180 
01181    // set flag for individual parameters
01182    fMCMCFlagsFillHistograms.push_back(true);
01183 
01184    // increase the number of parameters by one
01185    fMCMCNParameters++;
01186 
01187    // return the number of parameters
01188    return fMCMCNParameters;
01189 }
01190 
01191 // --------------------------------------------------------
01192 void BCEngineMCMC::MCMCInitializeMarkovChains()
01193 {
01194    // evaluate function at the starting point
01195    std::vector <double> x0;
01196 
01197    for (int j = 0; j < fMCMCNChains; ++j)
01198    {
01199       x0.clear();
01200       for (int i = 0; i < fMCMCNParameters; ++i)
01201          x0.push_back(fMCMCx[j * fMCMCNParameters + i]);
01202       fMCMCLogProbx[j] = this -> LogEval(x0);
01203    }
01204 
01205    x0.clear();
01206 }
01207 
01208 // --------------------------------------------------------
01209 int BCEngineMCMC::MCMCInitialize()
01210 {
01211    // reset variables
01212    fMCMCNIterations.clear();
01213    fMCMCNTrialsTrue.clear();
01214    fMCMCNTrialsFalse.clear();
01215    fMCMCTrialFunctionScaleFactor.clear();
01216    fMCMCMean.clear();
01217    fMCMCVariance.clear();
01218    fMCMCMeanx.clear();
01219    fMCMCVariancex.clear();
01220    fMCMCx.clear();
01221    fMCMCLogProbx.clear();
01222    fMCMCMaximumPoints.clear();
01223    fMCMCMaximumLogProb.clear();
01224    fMCMCNIterationsConvergenceGlobal = -1;
01225    fMCMCFlagConvergenceGlobal = false;
01226    fMCMCRValueParameters.clear();
01227 
01228    for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
01229       if (fMCMCH1Marginalized[i])
01230          delete fMCMCH1Marginalized[i];
01231 
01232    for (int i = 0; i < int(fMCMCH2Marginalized.size()); ++i)
01233       if (fMCMCH2Marginalized[i])
01234          delete fMCMCH2Marginalized[i];
01235 
01236    // clear plots
01237    fMCMCH1Marginalized.clear();
01238    fMCMCH2Marginalized.clear();
01239 
01240 // free memory for vectors
01241    fMCMCNIterations.assign(fMCMCNChains, 0);
01242    fMCMCMean.assign(fMCMCNChains, 0);
01243    fMCMCVariance.assign(fMCMCNChains, 0);
01244    fMCMCLogProbx.assign(fMCMCNChains, -1.0);
01245    fMCMCMaximumLogProb.assign(fMCMCNChains, -1.0);
01246 
01247    fMCMCNTrialsTrue.assign(fMCMCNChains * fMCMCNParameters, 0);
01248    fMCMCNTrialsFalse.assign(fMCMCNChains * fMCMCNParameters, 0);
01249    fMCMCMaximumPoints.assign(fMCMCNChains * fMCMCNParameters, 0.);
01250    fMCMCMeanx.assign(fMCMCNChains * fMCMCNParameters, 0);
01251    fMCMCVariancex.assign(fMCMCNChains * fMCMCNParameters, 0);
01252 
01253    fMCMCRValueParameters.assign(fMCMCNParameters, 0.);
01254 
01255    if (fMCMCTrialFunctionScaleFactorStart.size() == 0)
01256       fMCMCTrialFunctionScaleFactor.assign(fMCMCNChains * fMCMCNParameters, 1.0);
01257    else
01258       for (int i = 0; i < fMCMCNChains; ++i)
01259          for (int j = 0; j < fMCMCNParameters; ++j)
01260             fMCMCTrialFunctionScaleFactor.push_back(fMCMCTrialFunctionScaleFactorStart.at(j));
01261 
01262    // set initial position
01263    if (fMCMCFlagInitialPosition == 2) // user defined points
01264    {
01265       // define flag
01266       bool flag = true;
01267 
01268       // check the length of the array of initial positions
01269       if (int(fMCMCInitialPosition.size()) != (fMCMCNChains * fMCMCNParameters))
01270       {
01271          BCLog::OutError("BCEngine::MCMCInitialize : Length of vector containing initial positions does not have required length.");
01272          flag = false;
01273       }
01274 
01275       // check the boundaries
01276       if (flag)
01277       {
01278          for (int j = 0; j < fMCMCNChains; ++j)
01279             for (int i = 0; i < fMCMCNParameters; ++i)
01280                if (fMCMCInitialPosition[j * fMCMCNParameters + i] < fMCMCBoundaryMin[i] ||
01281                      fMCMCInitialPosition[j * fMCMCNParameters + i] > fMCMCBoundaryMax[i])
01282                {
01283                   BCLog::OutError("BCEngine::MCMCInitialize : Initial position out of boundaries.");
01284                   flag = false;
01285                }
01286       }
01287 
01288       // check flag
01289       if (!flag)
01290          fMCMCFlagInitialPosition = 1;
01291    }
01292 
01293    if (fMCMCFlagInitialPosition == 0) // center of the interval
01294       for (int j = 0; j < fMCMCNChains; ++j)
01295          for (int i = 0; i < fMCMCNParameters; ++i)
01296             fMCMCx.push_back(fMCMCBoundaryMin[i] + .5 * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01297 
01298    else if (fMCMCFlagInitialPosition == 2) // user defined
01299    {
01300       for (int j = 0; j < fMCMCNChains; ++j)
01301          for (int i = 0; i < fMCMCNParameters; ++i)
01302             fMCMCx.push_back(fMCMCInitialPosition.at(j * fMCMCNParameters + i));
01303    }
01304 
01305    else
01306       for (int j = 0; j < fMCMCNChains; ++j) // random number (default)
01307          for (int i = 0; i < fMCMCNParameters; ++i)
01308             fMCMCx.push_back(fMCMCBoundaryMin[i] + fMCMCRandom->Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01309 
01310    // copy the point of the first chain
01311    fMCMCxLocal.clear();
01312    for (int i = 0; i < fMCMCNParameters; ++i)
01313       fMCMCxLocal.push_back(fMCMCx[i]);
01314 
01315    // define 1-dimensional histograms for marginalization
01316    for(int i = 0; i < fMCMCNParameters; ++i)
01317    {
01318       TH1D * h1 = 0;
01319       if (fMCMCFlagsFillHistograms.at(i))
01320          h1 = new TH1D(TString::Format("h1_%d_parameter_%i", BCLog::GetHIndex() ,i), "",
01321                fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i]);
01322       fMCMCH1Marginalized.push_back(h1);
01323    }
01324 
01325    for(int i = 0; i < fMCMCNParameters; ++i)
01326       for (int k = 0; k < i; ++k)
01327       {
01328          TH2D * h2 = 0;
01329          if (fMCMCFlagsFillHistograms.at(i) && fMCMCFlagsFillHistograms.at(k))
01330             h2 = new TH2D(Form("h2_%d_parameters_%i_vs_%i", BCLog::GetHIndex(), i, k), "",
01331                   fMCMCH1NBins[k], fMCMCBoundaryMin[k], fMCMCBoundaryMax[k],
01332                   fMCMCH1NBins[i], fMCMCBoundaryMin[i], fMCMCBoundaryMax[i] );
01333          fMCMCH2Marginalized.push_back(h2);
01334       }
01335 
01336    fMCMCFlagPreRun = false;
01337    fMCMCFlagRun = false;
01338 
01339    return 1;
01340 }
01341 
01342 // ---------------------------------------------------------
01343 int BCEngineMCMC::SetMarginalized(int index, TH1D * h)
01344 {
01345    if((int)fMCMCH1Marginalized.size()<=index)
01346       return 0;
01347 
01348    if(h==0)
01349       return 0;
01350 
01351    if((int)fMCMCH1Marginalized.size()==index)
01352       fMCMCH1Marginalized.push_back(h);
01353    else
01354       fMCMCH1Marginalized[index]=h;
01355 
01356    return index;
01357 }
01358 
01359 // ---------------------------------------------------------
01360 int BCEngineMCMC::SetMarginalized(int index1, int index2, TH2D * h)
01361 {
01362    int counter = 0;
01363    int index = 0;
01364 
01365    // search for correct combination
01366    for(int i = 0; i < fMCMCNParameters; i++)
01367       for (int j = 0; j < i; ++j)
01368       {
01369          if(j == index1 && i == index2)
01370             index = counter;
01371          counter++;
01372       }
01373 
01374    if((int)fMCMCH2Marginalized.size()<=index)
01375       return 0;
01376 
01377    if(h==0)
01378       return 0;
01379 
01380    if((int)fMCMCH2Marginalized.size()==index)
01381       fMCMCH2Marginalized.push_back(h);
01382    else
01383       fMCMCH2Marginalized[index]=h;
01384 
01385    return index;
01386 }
01387 
01388 // ---------------------------------------------------------
01389 

Generated on Thu Feb 11 11:29:30 2010 for BayesianAnalysisToolkit by  doxygen 1.5.1