00001
00002
00003
00004
00005
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
00028 this -> MCMCSetValuesDefault();
00029
00030
00031 fMCMCRandom = new TRandom3(0);
00032 }
00033
00034
00035 BCEngineMCMC::BCEngineMCMC(int n)
00036 {
00037
00038 fMCMCNChains = n;
00039
00040
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
00101 if (fMCMCRandom)
00102 delete fMCMCRandom;
00103
00104
00105 for (int i = 0; i < int(fMCMCH1Marginalized.size()); ++i)
00106 if (fMCMCH1Marginalized[i])
00107 delete fMCMCH1Marginalized[i];
00108 fMCMCH1Marginalized.clear();
00109
00110
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
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
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
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
00174 std::vector <double> x;
00175
00176
00177 if (i < 0 || i >= fMCMCNChains)
00178 return x;
00179
00180
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
00193 this -> MCMCInitialize();
00194 }
00195
00196
00197 void BCEngineMCMC::MCMCSetInitialPositions(std::vector<double> x0s)
00198 {
00199
00200 fMCMCInitialPosition.clear();
00201
00202
00203 int n = int(x0s.size());
00204
00205 for (int i = 0; i < n; ++i)
00206 fMCMCInitialPosition.push_back(x0s.at(i));
00207
00208
00209 this -> MCMCSetFlagInitialPosition(2);
00210 }
00211
00212
00213 void BCEngineMCMC::MCMCSetInitialPositions(std::vector< std::vector<double> > x0s)
00214 {
00215
00216 std::vector <double> y0s;
00217
00218
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
00239 if (index < 0 || index > fMCMCNParameters)
00240 {
00241 BCLog::OutWarning("BCEngineMCMC :MCMCSetFlagFillHistograms. Index out of range.");
00242 return;
00243 }
00244
00245
00246 fMCMCFlagsFillHistograms[index] = flag;
00247 }
00248
00249
00250 void BCEngineMCMC::MCMCSetMarkovChainTrees(std::vector <TTree *> trees)
00251 {
00252
00253 fMCMCTrees.clear();
00254
00255
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
00268
00269
00270
00271
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
00280
00281
00282
00283
00284
00285 x[iparameter] = fMCMCRandom -> BreitWigner(0.0, fMCMCTrialFunctionScaleFactor[ichain * fMCMCNParameters + iparameter]);
00286 }
00287
00288
00289 std::vector <double> BCEngineMCMC::MCMCGetTrialFunctionScaleFactor(int ichain)
00290 {
00291
00292 std::vector <double> x;
00293
00294
00295 if (ichain < 0 || ichain >= fMCMCNChains)
00296 return x;
00297
00298
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
00309 if (ichain < 0 || ichain >= fMCMCNChains)
00310 return 0;
00311
00312
00313 if (ipar < 0 || ipar >= fMCMCNParameters)
00314 return 0;
00315
00316
00317 return fMCMCTrialFunctionScaleFactor.at(ichain * fMCMCNChains + ipar);
00318 }
00319
00320
00321 std::vector <double> BCEngineMCMC::MCMCGetx(int ichain)
00322 {
00323
00324 std::vector <double> x;
00325
00326
00327 if (ichain < 0 || ichain >= fMCMCNChains)
00328 return x;
00329
00330
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
00341 if (ichain < 0 || ichain >= fMCMCNChains)
00342 return 0;
00343
00344
00345 if (ipar < 0 || ipar >= fMCMCNParameters)
00346 return 0;
00347
00348
00349 return fMCMCx.at(ichain * fMCMCNParameters + ipar);
00350 }
00351
00352
00353 double BCEngineMCMC::MCMCGetLogProbx(int ichain)
00354 {
00355
00356 if (ichain < 0 || ichain >= fMCMCNChains)
00357 return -1;
00358
00359
00360 return fMCMCLogProbx.at(ichain);
00361 }
00362
00363
00364 bool BCEngineMCMC::MCMCGetProposalPointMetropolis(int chain, std::vector <double> &x)
00365 {
00366
00367 this -> MCMCTrialFunction(chain, x);
00368
00369
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
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
00385
00386 this -> MCMCTrialFunctionSingle(chain, parameter, x);
00387
00388
00389 for (int i = 0; i < fMCMCNParameters; ++i)
00390 if (i != parameter)
00391 x[i] = fMCMCx[chain * fMCMCNParameters + i];
00392
00393
00394 x[parameter] = fMCMCx[chain * fMCMCNParameters + parameter] + x[parameter] * (fMCMCBoundaryMax.at(parameter) - fMCMCBoundaryMin.at(parameter));
00395
00396
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
00408 int index = chain * fMCMCNParameters;
00409
00410
00411 fMCMCNIterations[chain]++;
00412
00413
00414 if (!this -> MCMCGetProposalPointMetropolis(chain, parameter, fMCMCxLocal))
00415 {
00416
00417 fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00418
00419 return false;
00420 }
00421
00422
00423 double p0 = fMCMCLogProbx[chain];
00424 double p1 = this -> LogEval(fMCMCxLocal);
00425
00426
00427 bool accept = false;
00428
00429
00430 if (p1 >= p0)
00431 accept = true;
00432
00433 else
00434 {
00435 double r = log(fMCMCRandom -> Rndm());
00436
00437 if(r < p1 - p0)
00438 accept = true;
00439 }
00440
00441
00442 if(accept)
00443 {
00444
00445 fMCMCNTrialsTrue[chain * fMCMCNParameters + parameter]++;
00446
00447
00448 for(int i = 0; i < fMCMCNParameters; ++i)
00449 {
00450
00451 fMCMCx[index + i] = fMCMCxLocal[i];
00452
00453
00454 fMCMCLogProbx[chain] = p1;
00455 }
00456 }
00457 else
00458 {
00459
00460 fMCMCNTrialsFalse[chain * fMCMCNParameters + parameter]++;
00461 }
00462
00463 return accept;
00464 }
00465
00466
00467 bool BCEngineMCMC::MCMCGetNewPointMetropolis(int chain)
00468 {
00469
00470 int index = chain * fMCMCNParameters;
00471
00472
00473 fMCMCNIterations[chain]++;
00474
00475
00476 if (!this -> MCMCGetProposalPointMetropolis(chain, fMCMCxLocal))
00477 {
00478
00479 for (int i = 0; i < fMCMCNParameters; ++i)
00480 fMCMCNTrialsFalse[chain * fMCMCNParameters + i]++;
00481
00482 return false;
00483 }
00484
00485
00486 double p0 = fMCMCLogProbx[chain];
00487 double p1 = this -> LogEval(fMCMCxLocal);
00488
00489
00490 bool accept = false;
00491
00492
00493 if (p1 >= p0)
00494 accept = true;
00495
00496
00497 else
00498 {
00499 double r = log(fMCMCRandom -> Rndm());
00500
00501 if(r < p1 - p0)
00502 accept = true;
00503 }
00504
00505
00506 if(accept)
00507 {
00508
00509 for (int i = 0; i < fMCMCNParameters; ++i)
00510 fMCMCNTrialsTrue[chain * fMCMCNParameters + i]++;
00511
00512
00513 for(int i = 0; i < fMCMCNParameters; ++i)
00514 {
00515
00516 fMCMCx[index + i] = fMCMCxLocal[i];
00517
00518
00519 fMCMCLogProbx[chain] = p1;
00520 }
00521 }
00522 else
00523 {
00524
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
00536 for (int i = 0; i < fMCMCNChains; ++i)
00537 {
00538
00539 if (fMCMCLogProbx[i] > fMCMCMaximumLogProb[i] || fMCMCNIterations[i] == 1)
00540 {
00541
00542 fMCMCMaximumLogProb[i] = fMCMCLogProbx[i];
00543
00544
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
00555 if (!fMCMCFlagFillHistograms)
00556 return;
00557
00558
00559 for (int i = 0; i < fMCMCNChains; ++i)
00560 {
00561
00562 for (int j = 0; j < fMCMCNParameters; ++j)
00563 if (fMCMCFlagsFillHistograms.at(j))
00564 fMCMCH1Marginalized[j] -> Fill(fMCMCx[i * fMCMCNParameters + j]);
00565
00566
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
00583 int npoints = fMCMCNTrialsTrue[0] + fMCMCNTrialsFalse[0];
00584
00585 if (fMCMCNChains > 1 && npoints > 1)
00586 {
00587
00588 bool flag_convergence = true;
00589
00590
00591 for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00592 {
00593 double sum = 0;
00594 double sum2 = 0;
00595 double sumv = 0;
00596
00597
00598 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00599 {
00600
00601 int index = ichains * fMCMCNParameters + iparameters;
00602
00603
00604 fMCMCMeanx[index] += (fMCMCx[index] - fMCMCMeanx[index]) / double(npoints);
00605
00606
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
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
00628 if (! ((fMCMCRValueParameters[iparameters]-1.0) < fMCMCRValueParametersCriterion))
00629 flag_convergence = false;
00630 }
00631
00632
00633 double sum = 0;
00634 double sum2 = 0;
00635 double sumv = 0;
00636
00637
00638 for (int i = 0; i < fMCMCNChains; ++i)
00639 {
00640
00641 fMCMCMean[i] += (fMCMCLogProbx[i] - fMCMCMean[i]) / double(npoints);
00642
00643
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
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
00666 if (! ((fMCMCRValue-1.0) < fMCMCRValueCriterion))
00667 flag_convergence = false;
00668
00669
00670 if (fMCMCNIterationsConvergenceGlobal == -1 && flag_convergence == true)
00671 fMCMCNIterationsConvergenceGlobal = fMCMCNIterations[0] / fMCMCNParameters;
00672 }
00673 }
00674
00675
00676 void BCEngineMCMC::MCMCUpdateStatisticsWriteChains()
00677 {
00678
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
00687
00688 return 0.0;
00689 }
00690
00691
00692 int BCEngineMCMC::MCMCMetropolisPreRun()
00693 {
00694
00695 BCLog::OutSummary("Pre-run Metropolis MCMC...");
00696
00697
00698 this -> MCMCInitialize();
00699 this -> MCMCInitializeMarkovChains();
00700
00701
00702 int ndigits = (int)log10(fMCMCNParameters) +1;
00703 if(ndigits<4)
00704 ndigits=4;
00705
00706
00707 if (fMCMCNIterationsBurnIn > 0)
00708 BCLog::OutDetail(Form(" --> Start burn-in run with %i iterations", fMCMCNIterationsBurnIn));
00709
00710
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
00719 else
00720 {
00721 for (int i = 0; i < fMCMCNIterationsBurnIn; ++i)
00722
00723 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00724
00725 this -> MCMCGetNewPointMetropolis(ichains);
00726 }
00727
00728
00729 this -> MCMCResetRunStatistics();
00730 fMCMCNIterationsConvergenceGlobal = -1;
00731
00732
00733 BCLog::OutSummary(Form(" --> Perform MCMC pre-run with %i chains, each with maximum %i iterations", fMCMCNChains, fMCMCNIterationsMax));
00734
00735
00736
00737 bool tempflag_writetofile = fMCMCFlagWriteChainToFile;
00738 fMCMCFlagWriteChainToFile = false;
00739
00740
00741 int counter = 0;
00742 int counterupdate = 0;
00743 bool convergence = false;
00744 bool flagefficiency = false;
00745
00746
00747 std::vector <double> efficiency;
00748 efficiency.assign(fMCMCNParameters * fMCMCNChains, 0.0);
00749
00750
00751 int updateLimit = ( fMCMCNIterationsUpdateMax<fMCMCNIterationsUpdate*(fMCMCNParameters+1) && fMCMCNIterationsUpdateMax>0 ) ?
00752 fMCMCNIterationsUpdateMax : fMCMCNIterationsUpdate*(fMCMCNParameters+1);
00753
00754
00755
00756
00757
00758
00759 while (counter < fMCMCNIterationsPreRunMin || (counter < fMCMCNIterationsMax && !(convergence && flagefficiency)))
00760 {
00761
00762 convergence = false;
00763
00764
00765 fMCMCNIterationsConvergenceGlobal = -1;
00766
00767
00768 if (fMCMCFlagOrderParameters)
00769 {
00770
00771 for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
00772 {
00773
00774 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00775 this -> MCMCGetNewPointMetropolis(ichains, iparameters);
00776
00777
00778 this -> MCMCUpdateStatisticsCheckMaximum();
00779 }
00780 }
00781
00782
00783 else
00784 {
00785
00786 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00787
00788 this -> MCMCGetNewPointMetropolis(ichains);
00789
00790
00791 this -> MCMCUpdateStatisticsCheckMaximum();
00792 }
00793
00794
00795 if ( counter > 0 && counter % fMCMCNIterationsUpdate == 0 )
00796 BCLog::OutDetail(Form(" --> Iteration %i", fMCMCNIterations[0]/fMCMCNParameters - 1));
00797
00798
00799
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
00808
00809
00810
00811
00812 fMCMCNIterationsConvergenceGlobal = -1;
00813
00814 this -> MCMCUpdateStatisticsTestConvergenceAllChains();
00815
00816
00817 if (fMCMCNIterationsConvergenceGlobal > 0)
00818 convergence = true;
00819
00820
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
00853
00854
00855
00856 flagefficiency = true;
00857
00858 bool flagprintefficiency = true;
00859
00860
00861 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
00862 {
00863
00864 for (int iparameter = 0; iparameter < fMCMCNParameters; ++iparameter)
00865 {
00866
00867 int ipc = ichains * fMCMCNParameters + iparameter;
00868
00869
00870 efficiency[ipc] = double(fMCMCNTrialsTrue[ipc]) / double(fMCMCNTrialsTrue[ipc] + fMCMCNTrialsFalse[ipc]);
00871
00872
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
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
00908 if ((efficiency[ipc] < fMCMCEfficiencyMin && fMCMCTrialFunctionScaleFactor[ipc] > .01)
00909 || (efficiency[ipc] > fMCMCEfficiencyMax && fMCMCTrialFunctionScaleFactor[ipc] < 1.))
00910 flagefficiency = false;
00911
00912
00913 counterupdate = 0;
00914 fMCMCNTrialsTrue[ipc] = 0;
00915 fMCMCNTrialsFalse[ipc] = 0;
00916 }
00917
00918
00919 fMCMCMean[ichains] = 0;
00920 fMCMCVariance[ichains] = 0;
00921 }
00922
00923
00924 if (flagefficiency)
00925 BCLog::OutDetail(Form(" * Efficiency status: Efficiencies within pre-defined ranges."));
00926
00927 }
00928
00929
00930 counter++;
00931 counterupdate++;
00932 }
00933
00934
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
00944
00945
00946
00947 if (fMCMCNIterationsConvergenceGlobal > 0)
00948 fMCMCFlagConvergenceGlobal = true;
00949 else
00950 fMCMCFlagConvergenceGlobal = false;
00951
00952
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
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
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
01006 fMCMCFlagWriteChainToFile = tempflag_writetofile;
01007
01008
01009 fMCMCFlagPreRun = true;
01010
01011
01012 return 1;
01013 }
01014
01015
01016 int BCEngineMCMC::MCMCMetropolis()
01017 {
01018
01019 if (!fMCMCFlagPreRun)
01020 this -> MCMCMetropolisPreRun();
01021
01022
01023 BCLog::OutSummary( "Run Metropolis MCMC...");
01024
01025
01026 this -> MCMCResetRunStatistics();
01027
01028
01029 BCLog::OutSummary(Form(" --> Perform MCMC run with %i chains, each with %i iterations.", fMCMCNChains, fMCMCNIterationsRun));
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
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
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
01058 if (fMCMCFlagOrderParameters)
01059 {
01060
01061 for (int iparameters = 0; iparameters < fMCMCNParameters; ++iparameters)
01062 {
01063
01064 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01065 this -> MCMCGetNewPointMetropolis(ichains, iparameters);
01066
01067
01068 this -> MCMCUpdateStatisticsCheckMaximum();
01069
01070
01071 if ( (fMCMCNParameters * iiterations + iparameters) % (fMCMCNLag * fMCMCNParameters) == 0)
01072 {
01073
01074 this -> MCMCUpdateStatisticsFillHistograms();
01075
01076
01077 if (fMCMCFlagWriteChainToFile)
01078 this -> MCMCUpdateStatisticsWriteChains();
01079
01080
01081 this -> MCMCIterationInterface();
01082 }
01083
01084 }
01085 }
01086
01087 else
01088 {
01089
01090 for (int ichains = 0; ichains < fMCMCNChains; ++ichains)
01091
01092 this -> MCMCGetNewPointMetropolis(ichains);
01093
01094
01095 this -> MCMCUpdateStatisticsCheckMaximum();
01096
01097
01098 if (iiterations % fMCMCNLag == 0)
01099 {
01100
01101 this -> MCMCUpdateStatisticsFillHistograms();
01102
01103
01104 if (fMCMCFlagWriteChainToFile)
01105 this -> MCMCUpdateStatisticsWriteChains();
01106
01107
01108 this -> MCMCIterationInterface();
01109 }
01110 }
01111
01112 }
01113
01114
01115 BCLog::OutSummary(Form(" --> Markov chains ran for %i iterations.", fMCMCNIterationsRun));
01116
01117
01118
01119
01120 double probmax = fMCMCMaximumLogProb.at(0);
01121 int probmaxindex = 0;
01122
01123
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
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
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
01178 fMCMCBoundaryMin.push_back(min);
01179 fMCMCBoundaryMax.push_back(max);
01180
01181
01182 fMCMCFlagsFillHistograms.push_back(true);
01183
01184
01185 fMCMCNParameters++;
01186
01187
01188 return fMCMCNParameters;
01189 }
01190
01191
01192 void BCEngineMCMC::MCMCInitializeMarkovChains()
01193 {
01194
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
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
01237 fMCMCH1Marginalized.clear();
01238 fMCMCH2Marginalized.clear();
01239
01240
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
01263 if (fMCMCFlagInitialPosition == 2)
01264 {
01265
01266 bool flag = true;
01267
01268
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
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
01289 if (!flag)
01290 fMCMCFlagInitialPosition = 1;
01291 }
01292
01293 if (fMCMCFlagInitialPosition == 0)
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)
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)
01307 for (int i = 0; i < fMCMCNParameters; ++i)
01308 fMCMCx.push_back(fMCMCBoundaryMin[i] + fMCMCRandom->Rndm() * (fMCMCBoundaryMax[i] - fMCMCBoundaryMin[i]));
01309
01310
01311 fMCMCxLocal.clear();
01312 for (int i = 0; i < fMCMCNParameters; ++i)
01313 fMCMCxLocal.push_back(fMCMCx[i]);
01314
01315
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
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