00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BAT/BCModel.h"
00011
00012 #include "BAT/BCDataPoint.h"
00013 #include "BAT/BCDataSet.h"
00014 #include "BAT/BCParameter.h"
00015 #include "BAT/BCH1D.h"
00016 #include "BAT/BCH2D.h"
00017 #include "BAT/BCGoFTest.h"
00018 #include "BAT/BCLog.h"
00019 #include "BAT/BCErrorCodes.h"
00020 #include "BAT/BCMath.h"
00021
00022 #include <TCanvas.h>
00023 #include <TPostScript.h>
00024 #include <TDirectory.h>
00025 #include <TFile.h>
00026 #include <TKey.h>
00027 #include <TTree.h>
00028 #include <TMath.h>
00029 #include <TGraph.h>
00030 #include <TH2D.h>
00031
00032 #include <fstream>
00033 #include <iomanip>
00034
00035
00036
00037 BCModel::BCModel(const char * name) : BCIntegrate()
00038 {
00039 fNormalization = -1.0;
00040 fDataSet = 0;
00041 fParameterSet = new BCParameterSet;
00042
00043 fIndex = -1;
00044 fPValue = -1;
00045 fPValueNDoF = -1;
00046 fChi2NDoF = -1;
00047
00048 fName = (char *) name;
00049 flag_ConditionalProbabilityEntry = true;
00050
00051 fDataPointUpperBoundaries = 0;
00052 fDataPointLowerBoundaries = 0;
00053
00054 fErrorBandXY = 0;
00055
00056 fGoFNChains = 5;
00057 fGoFNIterationsMax = 100000;
00058 fGoFNIterationsRun = 2000;
00059 }
00060
00061
00062
00063 BCModel::BCModel() : BCIntegrate()
00064 {
00065 fNormalization = -1.0;
00066 fDataSet = 0;
00067 fParameterSet = new BCParameterSet();
00068
00069 fIndex = -1;
00070 fPValue = -1;
00071 fPValueNDoF = -1;
00072 fChi2NDoF = -1;
00073
00074 fName = "model";
00075 fDataPointUpperBoundaries = 0;
00076 fDataPointLowerBoundaries = 0;
00077
00078 flag_ConditionalProbabilityEntry = true;
00079
00080 fGoFNChains = 5;
00081 fGoFNIterationsMax = 100000;
00082 fGoFNIterationsRun = 2000;
00083 }
00084
00085
00086
00087 BCModel::~BCModel()
00088 {
00089 delete fParameterSet;
00090
00091 if (fDataPointLowerBoundaries)
00092 delete fDataPointLowerBoundaries;
00093
00094 if (fDataPointUpperBoundaries)
00095 delete fDataPointUpperBoundaries;
00096 }
00097
00098
00099
00100 int BCModel::GetNDataPoints()
00101 {
00102 int npoints = 0;
00103 if (fDataSet)
00104 npoints = fDataSet -> GetNDataPoints();
00105 else
00106 {
00107 BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined.");
00108 return ERROR_NOEVENTS;
00109 }
00110
00111 return npoints;
00112 }
00113
00114
00115
00116 BCDataPoint * BCModel::GetDataPoint(unsigned int index)
00117 {
00118 if (fDataSet)
00119 return fDataSet -> GetDataPoint(index);
00120
00121 BCLog::OutWarning("BCModel::GetDataPoint. No data set defined.");
00122 return 0;
00123 }
00124
00125
00126
00127 BCParameter * BCModel::GetParameter(int index)
00128 {
00129 if (!fParameterSet)
00130 return 0;
00131
00132 if (index < 0 || index >= (int)this -> GetNParameters())
00133 {
00134 BCLog::OutWarning(
00135 Form("BCModel::GetParameter. Parameter index %d not within range.", index));
00136 return 0;
00137 }
00138
00139 return fParameterSet -> at(index);
00140 }
00141
00142
00143
00144 BCParameter * BCModel::GetParameter(const char * name)
00145 {
00146 if (!fParameterSet)
00147 return 0;
00148
00149 int index = -1;
00150 for (unsigned int i = 0; i < this->GetNParameters(); i++)
00151 if (name == this -> GetParameter(i) -> GetName())
00152 index = i;
00153
00154 if (index<0)
00155 {
00156 BCLog::OutWarning(Form(
00157 "BCModel::GetParameter : Model %s has no parameter named '%s'",
00158 (this -> GetName()).data(), name
00159 )
00160 );
00161 return 0;
00162 }
00163
00164 return this->GetParameter(index);
00165 }
00166
00167
00168
00169 void BCModel::SetNbins(const char * parname, int nbins)
00170 {
00171 BCParameter * p = this -> GetParameter(parname);
00172 if(!p)
00173 {
00174 BCLog::OutWarning(Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname));
00175 return;
00176 }
00177
00178 this -> BCIntegrate::SetNbins(nbins, p -> GetIndex());
00179 }
00180
00181
00182
00183 std::vector <double> BCModel::GetErrorBand(double level)
00184 {
00185 std::vector <double> errorband;
00186
00187 if (!fErrorBandXY)
00188 return errorband;
00189
00190 int nx = fErrorBandXY -> GetNbinsX();
00191 errorband.assign(nx, 0.0);
00192
00193
00194 for (int ix = 1; ix <= nx; ix++)
00195 {
00196 TH1D * temphist = fErrorBandXY -> ProjectionY("temphist", ix, ix);
00197
00198 int nprobSum = 1;
00199 double q[1];
00200 double probSum[1];
00201 probSum[0] = level;
00202
00203 temphist -> GetQuantiles(nprobSum, q, probSum);
00204
00205 errorband[ix-1] = q[0];
00206 }
00207
00208 return errorband;
00209 }
00210
00211
00212
00213 TGraph * BCModel::GetErrorBandGraph(double level1, double level2)
00214 {
00215 if (!fErrorBandXY)
00216 return 0;
00217
00218
00219 int nx = fErrorBandXY -> GetNbinsX();
00220
00221 TGraph * graph = new TGraph(2 * nx);
00222 graph -> SetFillStyle(1001);
00223 graph -> SetFillColor(kYellow);
00224
00225
00226 std::vector <double> ymin = this -> GetErrorBand(level1);
00227 std::vector <double> ymax = this -> GetErrorBand(level2);
00228
00229 for (int i = 0; i < nx; i++)
00230 {
00231 graph -> SetPoint(i, fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1), ymin.at(i));
00232 graph -> SetPoint(nx + i, fErrorBandXY -> GetXaxis() -> GetBinCenter(nx - i), ymax.at(nx - i - 1));
00233 }
00234
00235 return graph;
00236 }
00237
00238
00239
00240 TH2D * BCModel::GetErrorBandXY_yellow(double level, int nsmooth)
00241 {
00242 if (!fErrorBandXY)
00243 return 0;
00244
00245 int nx = fErrorBandXY -> GetNbinsX();
00246 int ny = fErrorBandXY -> GetNbinsY();
00247
00248
00249 TH2D * hist_tempxy = (TH2D*) fErrorBandXY -> Clone(TString::Format("%s_sub_%f.2",fErrorBandXY->GetName(),level));
00250 hist_tempxy -> Reset();
00251 hist_tempxy -> SetFillColor(kYellow);
00252
00253
00254 for (int ix = 1; ix < nx; ix++)
00255 {
00256 BCH1D * hist_temp = new BCH1D();
00257
00258 TH1D * hproj = fErrorBandXY -> ProjectionY("temphist", ix, ix);
00259 if(nsmooth>0)
00260 hproj->Smooth(nsmooth);
00261
00262 hist_temp -> SetHistogram(hproj);
00263
00264 TH1D * hist_temp_yellow = hist_temp -> GetSmallestIntervalHistogram(level);
00265
00266 for (int iy = 1; iy <= ny; ++iy)
00267 hist_tempxy -> SetBinContent(ix, iy, hist_temp_yellow -> GetBinContent(iy));
00268
00269 delete hist_temp_yellow;
00270 delete hist_temp;
00271 }
00272
00273 return hist_tempxy;
00274 }
00275
00276
00277
00278 TGraph * BCModel::GetFitFunctionGraph(std::vector <double> parameters)
00279 {
00280 if (!fErrorBandXY)
00281 return 0;
00282
00283
00284 int nx = fErrorBandXY -> GetNbinsX();
00285 TGraph * graph = new TGraph(nx);
00286
00287
00288 for (int i = 0; i < nx; i++)
00289 {
00290 double x = fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1);
00291
00292 std::vector <double> xvec;
00293 xvec.push_back(x);
00294 double y = this -> FitFunction(xvec, parameters);
00295 xvec.clear();
00296
00297 graph -> SetPoint(i, x, y);
00298 }
00299
00300 return graph;
00301 }
00302
00303
00304
00305 TGraph * BCModel::GetFitFunctionGraph(std::vector <double> parameters, double xmin, double xmax, int n)
00306 {
00307
00308 TGraph * graph = new TGraph(n+1);
00309
00310 double dx = (xmax-xmin)/(double)n;
00311
00312
00313 for (int i = 0; i <= n; i++)
00314 {
00315 double x = (double)i*dx+xmin;
00316 std::vector <double> xvec;
00317 xvec.push_back(x);
00318 double y = this -> FitFunction(xvec, parameters);
00319
00320 xvec.clear();
00321
00322 graph -> SetPoint(i, x, y);
00323 }
00324
00325 return graph;
00326 }
00327
00328
00329
00330 bool BCModel::GetFlagBoundaries()
00331 {
00332 if (!fDataPointLowerBoundaries)
00333 return false;
00334
00335 if (!fDataPointUpperBoundaries)
00336 return false;
00337
00338 if (fDataPointLowerBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues())
00339 return false;
00340
00341 if (fDataPointUpperBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues())
00342 return false;
00343
00344 return true;
00345 }
00346
00347
00348
00349 void BCModel::SetSingleDataPoint(BCDataPoint * datapoint)
00350 {
00351
00352 BCDataSet * dataset = new BCDataSet();
00353
00354
00355 dataset -> AddDataPoint(datapoint);
00356
00357
00358 this -> SetDataSet(dataset);
00359
00360 }
00361
00362
00363
00364 void BCModel::SetSingleDataPoint(BCDataSet * dataset, unsigned int index)
00365 {
00366 if (index < 0 || index > dataset -> GetNDataPoints())
00367 return;
00368
00369 this -> SetSingleDataPoint(dataset -> GetDataPoint(index));
00370 }
00371
00372
00373
00374 void BCModel::SetDataBoundaries(unsigned int index, double lowerboundary, double upperboundary, bool fixed)
00375 {
00376
00377 if (!fDataSet)
00378 {
00379 BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
00380 return;
00381 }
00382
00383
00384 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
00385 {
00386 BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
00387 return;
00388 }
00389
00390
00391 if (!fDataPointLowerBoundaries)
00392 fDataPointLowerBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues());
00393
00394 if (!fDataPointUpperBoundaries)
00395 fDataPointUpperBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues());
00396
00397 if (fDataFixedValues.size() == 0)
00398 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false);
00399
00400
00401 fDataPointLowerBoundaries -> SetValue(index, lowerboundary);
00402 fDataPointUpperBoundaries -> SetValue(index, upperboundary);
00403 fDataFixedValues[index] = fixed;
00404 }
00405
00406
00407
00408 void BCModel::SetErrorBandContinuous(bool flag)
00409 {
00410 fErrorBandContinuous = flag;
00411
00412 if (flag)
00413 return;
00414
00415
00416 fErrorBandX.clear();
00417
00418
00419 for (unsigned int i = 0; i < fDataSet -> GetNDataPoints(); ++i)
00420 fErrorBandX.push_back(fDataSet -> GetDataPoint(i) -> GetValue(fFitFunctionIndexX));
00421 }
00422
00423
00424
00425 int BCModel::AddParameter(const char * name, double lowerlimit, double upperlimit)
00426 {
00427
00428 BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit);
00429
00430 int flag_ok = this -> AddParameter(parameter);
00431 if (flag_ok)
00432 delete parameter;
00433
00434 return flag_ok;
00435 }
00436
00437
00438
00439 int BCModel::AddParameter(BCParameter * parameter)
00440 {
00441
00442 if (!fParameterSet)
00443 {
00444 BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
00445 return ERROR_PARAMETERSETDOESNOTEXIST;
00446 }
00447
00448
00449 int flag_exists = 0;
00450 for (unsigned int i = 0; i < this -> GetNParameters(); i++)
00451 if (this -> CompareStrings(parameter -> GetName().data(), this -> GetParameter(i) -> GetName().data()) == 0)
00452 flag_exists = -1;
00453
00454 if (flag_exists < 0)
00455 {
00456 BCLog::OutError(
00457 Form("BCModel::AddParameter : Parameter with name %s exists already. ", parameter -> GetName().data()));
00458 return ERROR_PARAMETEREXISTSALREADY;
00459 }
00460
00461
00462 unsigned int index = fParameterSet -> size();
00463 parameter -> SetIndex(index);
00464
00465
00466 fParameterSet -> push_back(parameter);
00467
00468
00469 this -> SetParameters(fParameterSet);
00470
00471 return 0;
00472 }
00473
00474
00475
00476 double BCModel::LogProbabilityNN(std::vector <double> parameters)
00477 {
00478
00479 double logprob = this -> LogLikelihood(parameters);
00480
00481
00482 logprob += this -> LogAPrioriProbability(parameters);
00483
00484 return logprob;
00485 }
00486
00487
00488
00489 double BCModel::LogProbability(std::vector <double> parameters)
00490 {
00491
00492 if (fNormalization<0. || fNormalization==0.)
00493 {
00494 BCLog::Out(BCLog::warning, BCLog::warning,
00495 "BCModel::LogProbability. Normalization not available or zero.");
00496 return 0.;
00497 }
00498
00499 return this -> LogProbabilityNN(parameters) - log(fNormalization);
00500 }
00501
00502
00503
00504 double BCModel::LogLikelihood(std::vector <double> parameters)
00505 {
00506 double logprob = 0.;
00507
00508
00509 for (unsigned int i=0;i<fDataSet -> GetNDataPoints();i++)
00510 {
00511 BCDataPoint * datapoint = this -> GetDataPoint(i);
00512 logprob += this -> LogConditionalProbabilityEntry(datapoint, parameters);
00513 }
00514
00515 return logprob;
00516 }
00517
00518
00519
00520 double BCModel::LogEval(std::vector <double> parameters)
00521 {
00522 return this -> LogProbabilityNN(parameters);
00523 }
00524
00525
00526
00527 double BCModel::EvalSampling(std::vector <double> parameters)
00528 {
00529 return this -> SamplingFunction(parameters);
00530 }
00531
00532
00533
00534 double BCModel::SamplingFunction(std::vector <double> parameters)
00535 {
00536 double probability = 1.0;
00537 for (std::vector<BCParameter*>::const_iterator it = fParameterSet -> begin(); it != fParameterSet -> end(); ++it)
00538 probability *= 1.0 / ((*it) -> GetUpperLimit() - (*it) -> GetLowerLimit());
00539 return probability;
00540 }
00541
00542
00543
00544 double BCModel::Normalize()
00545 {
00546 BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",this->GetName().data()));
00547
00548 unsigned int n = this -> GetNvar();
00549
00550
00551 if (n == 0)
00552 {
00553 this->SetParameters(fParameterSet);
00554 n = this->GetNvar();
00555 }
00556
00557
00558
00559 fNormalization = this -> Integrate();
00560
00561 BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization));
00562
00563 return fNormalization;
00564 }
00565
00566
00567
00568 int BCModel::CheckParameters(std::vector <double> parameters)
00569 {
00570
00571 if (!parameters.size() == fParameterSet -> size())
00572 return ERROR_INVALIDNUMBEROFPARAMETERS;
00573
00574
00575 for (unsigned int i = 0; i < fParameterSet -> size(); i++)
00576 {
00577 BCParameter * modelparameter = fParameterSet -> at(i);
00578
00579 if (modelparameter -> GetLowerLimit() > parameters.at(i) ||
00580 modelparameter -> GetUpperLimit() < parameters.at(i))
00581 {
00582 BCLog::OutError(
00583 Form("BCModel::CheckParameters : Parameter %s not within limits.", fParameterSet -> at(i) -> GetName().data()));
00584 return ERROR_PARAMETERNOTWITHINRANGE;
00585 }
00586 }
00587
00588 return 0;
00589 }
00590
00591
00592
00593 void BCModel::FindMode(std::vector<double> start)
00594 {
00595
00596
00597 BCLog::OutSummary(Form("Model \'%s\': Finding mode", this -> GetName().data()));
00598
00599
00600 this -> SetParameters(fParameterSet);
00601
00602 switch(this -> GetOptimizationMethod())
00603 {
00604 case BCIntegrate::kOptSA:
00605 this -> FindModeSA(start);
00606 return;
00607
00608 case BCIntegrate::kOptMinuit:
00609 {
00610 int printlevel = -1;
00611 if (BCLog::GetLogLevelScreen() <= BCLog::detail)
00612 printlevel = 0;
00613 if (BCLog::GetLogLevelScreen() <= BCLog::debug)
00614 printlevel = 1;
00615 this -> BCIntegrate::FindModeMinuit(start, printlevel);
00616 return;
00617 }
00618
00619 case BCIntegrate::kOptMetropolis:
00620 this -> MarginalizeAll();
00621 return;
00622 }
00623
00624 BCLog::OutError(
00625 Form("BCModel::FindMode : Invalid mode finding method: %d",
00626 this->GetOptimizationMethod()));
00627
00628 return;
00629 }
00630
00631
00632
00633 void BCModel::FindModeMinuit(std::vector<double> start, int printlevel)
00634 {
00635
00636 this -> SetParameters(fParameterSet);
00637
00638 this -> BCIntegrate::FindModeMinuit(start,printlevel);
00639 }
00640
00641
00642
00643 void BCModel::WriteMode(const char * file)
00644 {
00645 ofstream ofi(file);
00646 if(!ofi.is_open())
00647 {
00648 std::cerr<<"Couldn't open file "<<file<<std::endl;
00649 return;
00650 }
00651
00652 int npar = fParameterSet -> size();
00653 for (int i=0; i<npar; i++)
00654 ofi<<fBestFitParameters.at(i)<<std::endl;
00655
00656 ofi<<std::endl;
00657 ofi<<"#######################################################################"<<std::endl;
00658 ofi<<"#"<<std::endl;
00659 ofi<<"# This file was created automatically by BCModel::WriteMode() call."<<std::endl;
00660 ofi<<"# It can be read in by call to BCModel::ReadMode()."<<std::endl;
00661 ofi<<"# Do not modify it unless you know what you're doing."<<std::endl;
00662 ofi<<"#"<<std::endl;
00663 ofi<<"#######################################################################"<<std::endl;
00664 ofi<<"#"<<std::endl;
00665 ofi<<"# Best fit parameters (mode) for model:"<<std::endl;
00666 ofi<<"# \'"<<fName.data()<<"\'"<<std::endl;
00667 ofi<<"#"<<std::endl;
00668 ofi<<"# Number of parameters: "<<npar<<std::endl;
00669 ofi<<"# Parameters ordered as above:"<<std::endl;
00670
00671 for (int i=0; i<npar; i++)
00672 {
00673 ofi<<"# "<<i<<": ";
00674 ofi<<fParameterSet->at(i)->GetName().data()<<" = ";
00675 ofi<<fBestFitParameters.at(i)<<std::endl;
00676 }
00677
00678 ofi<<"#"<<std::endl;
00679 ofi<<"########################################################################"<<std::endl;
00680 }
00681
00682
00683
00684 int BCModel::ReadMode(const char * file)
00685 {
00686 ifstream ifi(file);
00687 if(!ifi.is_open())
00688 {
00689 BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.",file));
00690 return 0;
00691 }
00692
00693 int npar = fParameterSet -> size();
00694 std::vector <double> mode;
00695
00696 int i=0;
00697 while (i<npar && !ifi.eof())
00698 {
00699 double a;
00700 ifi>>a;
00701 mode.push_back(a);
00702 i++;
00703 }
00704
00705 if(i<npar)
00706 {
00707 BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.",file));
00708 BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.",npar,i));
00709 return 0;
00710 }
00711
00712 BCLog::OutSummary(Form("# Read in best fit parameters (mode) for model \'%s\' from file %s:",fName.data(),file));
00713 this->SetMode(mode);
00714 for(int j=0 ; j<npar; j++)
00715 BCLog::OutSummary(Form("# -> Parameter %d : %s = %e", j, fParameterSet->at(j)->GetName().data(), fBestFitParameters[j]));
00716
00717 BCLog::OutWarning("# ! Best fit values obtained before this call will be overwritten !");
00718
00719 return npar;
00720 }
00721
00722
00723
00724 int BCModel::MarginalizeAll()
00725 {
00726 BCLog::OutSummary(Form("Running MCMC for model \'%s\'",this->GetName().data()));
00727
00728
00729 double dx = 0.0;
00730 double dy = 0.0;
00731
00732 if (fFitFunctionIndexX >= 0)
00733 {
00734 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX))
00735 / (double)fErrorBandNbinsX;
00736
00737 dy = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY))
00738 / (double)fErrorBandNbinsY;
00739
00740 fErrorBandXY = new TH2D(
00741 TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "",
00742 fErrorBandNbinsX,
00743 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - .5 * dx,
00744 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + .5 * dx,
00745 fErrorBandNbinsY,
00746 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - .5 * dy,
00747 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + .5 * dy);
00748 fErrorBandXY -> SetStats(kFALSE);
00749
00750 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
00751 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
00752 fErrorBandXY -> SetBinContent(ix, iy, 0.0);
00753 }
00754
00755 this -> MCMCMetropolis();
00756 this -> FindModeMCMC();
00757
00758
00759
00760 return 1;
00761 }
00762
00763
00764
00765
00766 BCH1D * BCModel::GetMarginalized(BCParameter * parameter)
00767 {
00768 if (!parameter)
00769 {
00770 BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist.");
00771 return 0;
00772 }
00773
00774 int index = parameter -> GetIndex();
00775 if (fMCMCFlagsFillHistograms[index] == false)
00776 {
00777 BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' not filled.",parameter->GetName().data()));
00778 return 0;
00779 }
00780
00781
00782 TH1D * hist = this -> MCMCGetH1Marginalized(index);
00783 if(!hist)
00784 return 0;
00785
00786
00787 hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data()));
00788 hist->SetXTitle(parameter->GetName().data());
00789 hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data()));
00790 hist->SetStats(kFALSE);
00791
00792
00793 BCH1D * hprob = new BCH1D();
00794 hprob->SetHistogram(hist);
00795
00796
00797 double bestfit = hprob->GetMode();
00798
00799 if (fBestFitParametersMarginalized.size() == 0)
00800 for (unsigned int i = 0; i < GetNParameters(); i++)
00801 fBestFitParametersMarginalized.push_back(0.);
00802
00803 fBestFitParametersMarginalized[index] = bestfit;
00804
00805 hprob->SetGlobalMode(fBestFitParameters.at(index));
00806
00807 return hprob;
00808 }
00809
00810
00811
00812 int BCModel::ReadMarginalizedFromFile(const char * file)
00813 {
00814 TFile * froot = new TFile(file);
00815 if(!froot->IsOpen())
00816 {
00817 BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.",file));
00818 return 0;
00819 }
00820
00821
00822
00823
00824
00825
00826
00827 this -> MCMCInitialize();
00828
00829 int k=0;
00830 int n=this->GetNParameters();
00831 for(int i=0;i<n;i++)
00832 {
00833 BCParameter * a = this -> GetParameter(i);
00834 TKey * key = froot -> GetKey(TString::Format("hist_%s_%s", this -> GetName().data(), a -> GetName().data()));
00835 if(key)
00836 {
00837 TH1D * h1 = (TH1D*) key -> ReadObjectAny(TH1D::Class());
00838 h1->SetDirectory(0);
00839 if(this->SetMarginalized(i,h1))
00840 k++;
00841 }
00842 else
00843 BCLog::OutWarning(Form(
00844 "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.",
00845 this -> GetName().data(), a -> GetName().data(), file));
00846 }
00847
00848 for(int i=0;i<n-1;i++)
00849 {
00850 for(int j=i+1;j<n;j++)
00851 {
00852 BCParameter * a = this -> GetParameter(i);
00853 BCParameter * b = this -> GetParameter(j);
00854 TKey * key = froot -> GetKey(TString::Format("hist_%s_%s_%s", this -> GetName().data(), a -> GetName().data(), b -> GetName().data()));
00855 if(key)
00856 {
00857 TH2D * h2 = (TH2D*) key -> ReadObjectAny(TH2D::Class());
00858 h2->SetDirectory(0);
00859 if(this->SetMarginalized(i,j,h2))
00860 k++;
00861 }
00862 else
00863 BCLog::OutWarning(Form(
00864 "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.",
00865 this -> GetName().data(), a -> GetName().data(), b -> GetName().data(), file));
00866 }
00867 }
00868
00869 froot->Close();
00870
00871 return k;
00872 }
00873
00874
00875
00876 int BCModel::ReadErrorBandFromFile(const char * file)
00877 {
00878 TFile * froot = new TFile(file);
00879 if(!froot->IsOpen())
00880 {
00881 BCLog::OutWarning(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.",file));
00882 return 0;
00883 }
00884
00885 int r=0;
00886
00887 TH2D * h2 = (TH2D*)froot->Get("errorbandxy");
00888 if(h2)
00889 {
00890 h2->SetDirectory(0);
00891 h2->SetName(TString::Format("errorbandxy_%d",BCLog::GetHIndex()));
00892 this->SetErrorBandHisto(h2);
00893 r=1;
00894 }
00895 else
00896 BCLog::OutWarning(Form(
00897 "BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",
00898 file));
00899
00900 froot->Close();
00901
00902 return r;
00903 }
00904
00905
00906
00907 int BCModel::PrintAllMarginalized1D(const char * filebase)
00908 {
00909 if(fMCMCH1Marginalized.size()==0)
00910 {
00911 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00912 return 0;
00913 }
00914
00915 int n=this->GetNParameters();
00916 for(int i=0;i<n;i++)
00917 {
00918 BCParameter * a = this->GetParameter(i);
00919 if (this -> GetMarginalized(a))
00920 this -> GetMarginalized(a) -> Print(Form("%s_1D_%s.ps",filebase,a->GetName().data()));
00921 }
00922
00923 return n;
00924 }
00925
00926
00927
00928 int BCModel::PrintAllMarginalized2D(const char * filebase)
00929 {
00930 if(fMCMCH2Marginalized.size()==0)
00931 {
00932 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00933 return 0;
00934 }
00935
00936 int k=0;
00937 int n=this->GetNParameters();
00938 for(int i=0;i<n-1;i++)
00939 {
00940 for(int j=i+1;j<n;j++)
00941 {
00942 BCParameter * a = this->GetParameter(i);
00943 BCParameter * b = this->GetParameter(j);
00944
00945 double meana = (a -> GetLowerLimit() + a -> GetUpperLimit()) / 2.0;
00946 double deltaa = (a -> GetUpperLimit() - a -> GetLowerLimit());
00947 if (deltaa <= 1e-7 * meana)
00948 continue;
00949
00950 double meanb = (b -> GetLowerLimit() + b -> GetUpperLimit()) / 2.0;
00951 double deltab = (b -> GetUpperLimit() - b -> GetLowerLimit());
00952 if (deltab <= 1e-7 * meanb)
00953 continue;
00954
00955 if (this -> GetMarginalized(a,b))
00956 this -> GetMarginalized(a,b) -> Print(Form("%s_2D_%s_%s.ps",filebase,a->GetName().data(),b->GetName().data()));
00957 k++;
00958 }
00959 }
00960
00961 return k;
00962 }
00963
00964
00965
00966 int BCModel::PrintAllMarginalized(const char * file, unsigned int hdiv, unsigned int vdiv)
00967 {
00968 if(!fMCMCFlagFillHistograms)
00969 {
00970 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled.");
00971 return 0;
00972 }
00973
00974 int npar = GetNParameters();
00975
00976 if(fMCMCH1Marginalized.size()==0 || (fMCMCH2Marginalized.size()==0 && npar>1))
00977 {
00978 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00979 return 0;
00980 }
00981
00982
00983 if (fMCMCH1Marginalized.size()==1 && fMCMCH2Marginalized.size()==0)
00984 {
00985 BCParameter * a = GetParameter(0);
00986 GetMarginalized(a)->Print(file);
00987 return 1;
00988 }
00989
00990 int c_width=600;
00991 int c_height=600;
00992
00993 int type = 112;
00994
00995 if (hdiv > vdiv)
00996 {
00997 if(hdiv>3)
00998 {
00999 c_width=1000;
01000 c_height=700;
01001 }
01002 else
01003 {
01004 c_width=800;
01005 c_height=600;
01006 }
01007 }
01008 else if(hdiv < vdiv)
01009 {
01010 if(hdiv>3)
01011 {
01012 c_height=1000;
01013 c_width=700;
01014 }
01015 else
01016 {
01017 c_height=800;
01018 c_width=600;
01019 }
01020 type=111;
01021 }
01022
01023
01024 int nplots2d = npar * (npar-1)/2;
01025 int nplots = npar + nplots2d;
01026
01027
01028 BCLog::OutSummary(Form(
01029 "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
01030 npar,nplots2d,nplots,file));
01031 if(nplots>100)
01032 BCLog::OutDetail("This can take a while...");
01033
01034
01035 TCanvas * c = new TCanvas( "c","canvas",c_width,c_height);
01036
01037 TPostScript * ps = new TPostScript(file,type);
01038
01039 if(type==112)
01040 ps->Range(24,16);
01041 else
01042 ps->Range(16,24);
01043
01044
01045 ps->NewPage();
01046 c->cd();
01047 c->Clear();
01048 c->Divide(hdiv,vdiv);
01049
01050 int n=0;
01051 for(int i=0;i<npar;i++)
01052 {
01053
01054 BCParameter * a = GetParameter(i);
01055
01056
01057 if ( !GetMarginalized(a) )
01058 continue;
01059
01060
01061 if ( GetMarginalized(a)->GetHistogram()->Integral() <= 0)
01062 continue;
01063
01064
01065 if(i!=0 && i%(hdiv*vdiv)==0)
01066 {
01067 c->Update();
01068 ps->NewPage();
01069 c->cd();
01070 c->Clear();
01071 c->Divide(hdiv,vdiv);
01072 }
01073
01074
01075 c->cd(i%(hdiv*vdiv)+1);
01076
01077 this -> GetMarginalized(a) -> Draw();
01078 n++;
01079
01080 if(n%100==0)
01081 BCLog::OutDetail(Form(" --> %d plots done",n));
01082 }
01083
01084 if (n > 0)
01085 {
01086 c->Update();
01087
01088
01089 ps->NewPage();
01090 c->cd();
01091 c->Clear();
01092 }
01093
01094 c->Divide(hdiv,vdiv);
01095
01096 int k=0;
01097 for(int i=0;i<npar-1;i++)
01098 {
01099 if(!fMCMCFlagsFillHistograms[i])
01100 continue;
01101 for(int j=i+1;j<npar;j++)
01102 {
01103 if(!fMCMCFlagsFillHistograms[j])
01104 continue;
01105
01106
01107 BCParameter * a = GetParameter(i);
01108 BCParameter * b = GetParameter(j);
01109
01110
01111 if ( !GetMarginalized(a,b) )
01112 continue;
01113
01114
01115 if ( GetMarginalized(a,b)->GetHistogram()->Integral() <= 0 )
01116 continue;
01117
01118
01119 if(k!=0 && k%(hdiv*vdiv)==0)
01120 {
01121 c->Update();
01122 ps->NewPage();
01123 c->cd();
01124 c->Clear();
01125 c->Divide(hdiv,vdiv);
01126 }
01127
01128
01129 c->cd(k%(hdiv*vdiv)+1);
01130
01131 double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
01132 double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
01133 if (deltaa <= 1e-7 * meana)
01134 continue;
01135
01136 double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
01137 double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
01138 if (deltab <= 1e-7 * meanb)
01139 continue;
01140
01141 GetMarginalized(a,b)->Draw(52);
01142 k++;
01143
01144 if((n+k)%100==0)
01145 BCLog::OutDetail(Form(" --> %d plots done",n+k));
01146 }
01147 }
01148
01149 if( (n+k)>100 && (n+k)%100 != 0 )
01150 BCLog::OutDetail(Form(" --> %d plots done",n+k));
01151
01152 c->Update();
01153 ps->Close();
01154
01155 delete c;
01156 delete ps;
01157
01158
01159 return n+k;
01160 }
01161
01162
01163
01164 BCH2D * BCModel::GetMarginalized(BCParameter * par1, BCParameter * par2)
01165 {
01166 int index1 = par1->GetIndex();
01167 int index2 = par2->GetIndex();
01168
01169 if (fMCMCFlagsFillHistograms[index1]==false || fMCMCFlagsFillHistograms[index2]==false )
01170 {
01171 BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.",
01172 par1->GetName().data(),par2->GetName().data()));
01173 return 0;
01174 }
01175
01176 if (index1 == index2)
01177 {
01178 BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
01179 return 0;
01180 }
01181
01182 BCParameter * npar1 = par1;
01183 BCParameter * npar2 = par2;
01184
01185 if (index1>index2)
01186 {
01187 npar1 = par2;
01188 npar2 = par1;
01189
01190 int itmp = index1;
01191 index1 = index2;
01192 index2 = itmp;
01193 }
01194
01195
01196 TH2D * hist = this -> MCMCGetH2Marginalized(index1, index2);
01197
01198 if(hist==0)
01199 return 0;
01200
01201 BCH2D * hprob = new BCH2D();
01202
01203
01204 hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data()));
01205 hist->SetXTitle(Form("%s", npar1->GetName().data()));
01206 hist->SetYTitle(Form("%s", npar2->GetName().data()));
01207 hist->SetStats(kFALSE);
01208
01209 double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
01210 hprob->SetGlobalMode(gmode);
01211
01212
01213 hprob->SetHistogram(hist);
01214
01215 return hprob;
01216 }
01217
01218
01219
01220 double BCModel::GetPvalueFromChi2(std::vector<double> par, int sigma_index)
01221 {
01222 double ll = this -> LogLikelihood(par);
01223 int n = this -> GetNDataPoints();
01224
01225 double sum_sigma=0;
01226 for (int i=0;i<n;i++)
01227 sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index));
01228
01229 double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma);
01230
01231 fPValue = TMath::Prob(chi2,n);
01232
01233 return fPValue;
01234 }
01235
01236
01237
01238 double BCModel::GetPvalueFromChi2NDoF(std::vector<double> par, int sigma_index)
01239 {
01240 double ll = this -> LogLikelihood(par);
01241 int n = this -> GetNDataPoints();
01242 int npar = this -> GetNParameters();
01243
01244 double sum_sigma=0;
01245 for (int i=0;i<n;i++)
01246 sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index));
01247
01248 double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma);
01249
01250 fChi2NDoF = chi2/double(n-npar);
01251 fPValueNDoF = TMath::Prob(chi2,n-npar);
01252
01253 return fPValueNDoF;
01254 }
01255
01256
01257
01258 BCH1D * BCModel::CalculatePValue(std::vector<double> par, bool flag_histogram)
01259 {
01260 BCH1D * hist = 0;
01261
01262
01263 BCLog::OutSummary("Do goodness-of-fit-test");
01264
01265
01266 BCGoFTest * goftest = new BCGoFTest("modeltest");
01267
01268
01269 goftest -> SetTestModel(this);
01270
01271
01272
01273 if (!goftest -> SetTestPoint(par))
01274 return 0;
01275
01276
01277
01278 goftest->MCMCSetFlagFillHistograms(false);
01279
01280
01281 goftest -> MCMCSetNChains(fGoFNChains);
01282 goftest -> MCMCSetNIterationsMax(fGoFNIterationsMax);
01283 goftest -> MCMCSetNIterationsRun(fGoFNIterationsRun);
01284
01285
01286 fPValue = goftest -> GetCalculatedPValue(flag_histogram);
01287
01288
01289 if (flag_histogram)
01290 {
01291 hist = new BCH1D();
01292 hist -> SetHistogram(goftest -> GetHistogramLogProb());
01293 }
01294
01295
01296 delete goftest;
01297
01298
01299 return hist;
01300 }
01301
01302
01303
01304 void BCModel::CorrelateDataPointValues(std::vector<double> &x)
01305 {
01306
01307 }
01308
01309
01310
01311 double BCModel::HessianMatrixElement(BCParameter * par1, BCParameter * par2, std::vector<double> point)
01312 {
01313
01314 if (point.size() != this -> GetNParameters())
01315 {
01316 BCLog::OutWarning("BCModel::HessianMatrixElement. Invalid number of entries in the vector.");
01317 return -1;
01318 }
01319
01320
01321 double nsteps = 1e5;
01322 double dx1 = par1 -> GetRangeWidth() / nsteps;
01323 double dx2 = par2 -> GetRangeWidth() / nsteps;
01324
01325
01326 std::vector<double> xpp = point;
01327 std::vector<double> xpm = point;
01328 std::vector<double> xmp = point;
01329 std::vector<double> xmm = point;
01330
01331 int idx1 = par1 -> GetIndex();
01332 int idx2 = par2 -> GetIndex();
01333
01334 xpp[idx1] += dx1;
01335 xpp[idx2] += dx2;
01336
01337 xpm[idx1] += dx1;
01338 xpm[idx2] -= dx2;
01339
01340 xmp[idx1] -= dx1;
01341 xmp[idx2] += dx2;
01342
01343 xmm[idx1] -= dx1;
01344 xmm[idx2] -= dx2;
01345
01346
01347 double ppp = this -> Likelihood(xpp);
01348 double ppm = this -> Likelihood(xpm);
01349 double pmp = this -> Likelihood(xmp);
01350 double pmm = this -> Likelihood(xmm);
01351
01352
01353 return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
01354 }
01355
01356
01357
01358 void BCModel::FixDataAxis(unsigned int index, bool fixed)
01359 {
01360
01361 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
01362 {
01363 BCLog::OutWarning("BCModel::FixDataAxis. Index out of range.");
01364 return;
01365 }
01366
01367 if (fDataFixedValues.size() == 0)
01368 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false);
01369
01370 fDataFixedValues[index] = fixed;
01371 }
01372
01373
01374
01375 bool BCModel::GetFixedDataAxis(unsigned int index)
01376 {
01377
01378 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
01379 {
01380 BCLog::OutWarning("BCModel::GetFixedDataAxis. Index out of range.");
01381 return false;
01382 }
01383
01384 return fDataFixedValues.at(index);
01385 }
01386
01387
01388
01389 void BCModel::PrintSummary()
01390 {
01391 int nparameters = this -> GetNParameters();
01392
01393
01394 std::cout
01395 << std::endl
01396 << " ---------------------------------" << std::endl
01397 << " Model : " << fName.data() << std::endl
01398 << " ---------------------------------"<< std::endl
01399 << " Index : " << fIndex << std::endl
01400 << " Number of parameters : " << nparameters << std::endl
01401 << std::endl
01402 << " - Parameters : " << std::endl
01403 << std::endl;
01404
01405
01406 for (int i=0; i<nparameters; i++)
01407 fParameterSet -> at(i) -> PrintSummary();
01408
01409
01410 if (this -> GetBestFitParameters().size() > 0)
01411 {
01412 std::cout
01413 << std::endl
01414 << " - Best fit parameters :" << std::endl
01415 << std::endl;
01416
01417 for (int i=0; i<nparameters; i++)
01418 {
01419 std::cout
01420 << " " << fParameterSet -> at(i) -> GetName().data()
01421 << " = " << this -> GetBestFitParameter(i)
01422 << " (overall)" << std::endl;
01423 if ((int)fBestFitParametersMarginalized.size() == nparameters)
01424 std::cout
01425 << " " << fParameterSet -> at(i) -> GetName().data()
01426 << " = " << this -> GetBestFitParameterMarginalized(i)
01427 << " (marginalized)" << std::endl;
01428 }
01429 }
01430
01431 std::cout << std::endl;
01432
01433
01434 if (fPValue >= 0)
01435 {
01436 double likelihood = this -> Likelihood(this -> GetBestFitParameters());
01437 std::cout
01438 << " - Model testing:" << std::endl
01439 << std::endl
01440 << " p(data|lambda*) = " << likelihood << std::endl
01441 << " p-value = " << fPValue << std::endl
01442 << std::endl;
01443 }
01444
01445
01446 if (fNormalization > 0)
01447 std::cout << " Normalization : " << fNormalization << std::endl;
01448 }
01449
01450
01451
01452 void BCModel::PrintResults(const char * file)
01453 {
01454
01455
01456
01457 ofstream ofi(file);
01458
01459
01460 if(!ofi.is_open())
01461 {
01462 std::cerr << "Couldn't open file " << file <<std::endl;
01463 return;
01464 }
01465
01466
01467 int npar = MCMCGetNParameters();
01468 int nchains = MCMCGetNChains();
01469
01470
01471 bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0;
01472
01473 ofi
01474 << std::endl
01475 << " -----------------------------------------------------" << std::endl
01476 << " Summary" << std::endl
01477 << " -----------------------------------------------------" << std::endl
01478 << std::endl;
01479
01480 ofi
01481 << " Model summary" << std::endl
01482 << " =============" << std::endl
01483 << " Model: " << fName.data() << std::endl
01484 << " Number of parameters: " << npar << std::endl
01485 << " List of Parameters and ranges:" << std::endl;
01486 for (int i = 0; i < npar; ++i)
01487 ofi
01488 << " (" << i << ") Parameter \""
01489 << fParameterSet -> at(i) -> GetName().data() << "\""
01490 << ": " << fParameterSet -> at(i) -> GetLowerLimit()
01491 << " - "
01492 << fParameterSet -> at(i) -> GetUpperLimit() << std::endl;
01493 ofi << std::endl;
01494
01495
01496 if (!flag_conv && fMCMCFlagRun)
01497 ofi
01498 << " WARNING: the Markov Chain did not converge!" << std::endl
01499 << " Be cautious using the following results!" << std::endl
01500 << std::endl;
01501
01502
01503 if ( fMCMCFlagRun && fMCMCFlagFillHistograms)
01504 {
01505 ofi
01506 << " Results of the marginalization" << std::endl
01507 << " ==============================" << std::endl
01508 << " List of parameters and properties of the marginalized" << std::endl
01509 << " distributions:" << std::endl;
01510 for (int i = 0; i < npar; ++i)
01511 {
01512 if (!fMCMCFlagsFillHistograms.at(i))
01513 continue;
01514
01515 BCH1D * bch1d = GetMarginalized(fParameterSet -> at(i));
01516
01517 ofi
01518 << " (" << i << ") Parameter \""
01519 << fParameterSet->at(i)->GetName().data() << "\"" << std::endl
01520 << " Mean +- RMS: "
01521 << std::setprecision(4) << bch1d->GetMean()
01522 << " +- "
01523 << std::setprecision(4) << bch1d->GetRMS() << std::endl
01524 << " Median +- sigma: "
01525 << std::setprecision(4) << bch1d->GetMedian()
01526 << " + " << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian()
01527 << " - " << std::setprecision(4) << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl
01528 << " (Marginalized) mode: " << bch1d->GetMode() << std::endl
01529 << " Smallest interval(s) containing 68% and local modes:" << std::endl;
01530
01531 std::vector <double> v;
01532 v = bch1d->GetSmallestIntervals(0.68);
01533 int ninter = int(v.size());
01534
01535 for (int j = 0; j < ninter; j+=5)
01536 ofi
01537 << " " << v.at(j) << " - " << v.at(j+1)
01538 << " (local mode at " << v.at(j+3)
01539 << " with rel. height " << v.at(j+2)
01540 << "; rel. area " << v.at(j+4) << ")"
01541 << std::endl;
01542
01543 ofi
01544 << " 5% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.05) << std::endl
01545 << " 10% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.10) << std::endl
01546 << " 16% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.16) << std::endl
01547 << " 84% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.85) << std::endl
01548 << " 90% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.90) << std::endl
01549 << " 95% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.95) << std::endl
01550 << std::endl;
01551 }
01552 ofi << std::endl;
01553 }
01554
01555 ofi
01556 << " Results of the optimization" << std::endl
01557 << " ===========================" << std::endl
01558 << " Optimization algorithm used: ";
01559 switch(GetOptimizationMethodMode())
01560 {
01561 case BCIntegrate::kOptSA:
01562 ofi << " Simulated Annealing" << std::endl;
01563 break;
01564 case BCIntegrate::kOptMinuit:
01565 ofi << " Minuit" << std::endl;
01566 break;
01567 case BCIntegrate::kOptMetropolis:
01568 ofi << " MCMC" << std::endl;
01569 break;
01570 }
01571
01572 if (int(fBestFitParameters.size()) > 0)
01573 {
01574 ofi << " List of parameters and global mode:" << std::endl;
01575 for (int i = 0; i < npar; ++i)
01576 ofi
01577 << " (" << i << ") Parameter \""
01578 << fParameterSet->at(i)->GetName().data() << "\": "
01579 << fBestFitParameters.at(i) << " +- " << fBestFitParameterErrors.at(i) << std::endl;
01580 ofi << std::endl;
01581 }
01582 else
01583 {
01584 ofi << " No best fit information available." << std::endl;
01585 ofi << std::endl;
01586 }
01587
01588 if (fPValue >= 0.)
01589 {
01590 ofi
01591 << " Results of the model test" << std::endl
01592 << " =========================" << std::endl
01593 << " p-value at global mode: " << fPValue << std::endl
01594 << std::endl;
01595 }
01596
01597 if ( fMCMCFlagRun )
01598 {
01599 ofi
01600 << " Status of the MCMC" << std::endl
01601 << " ==================" << std::endl
01602 << " Convergence reached: " << (flag_conv ? "yes" : "no") << std::endl;
01603
01604 if (flag_conv)
01605 ofi << " Number of iterations until convergence: " << MCMCGetNIterationsConvergenceGlobal() << std::endl;
01606 else
01607 ofi
01608 << " WARNING: the Markov Chain did not converge! Be\n"
01609 << " cautious using the following results!" << std::endl
01610 << std::endl;
01611 ofi
01612 << " Number of chains: " << MCMCGetNChains() << std::endl
01613 << " Number of iterations per chain: " << MCMCGetNIterationsRun() << std::endl
01614 << " Average efficiencies:" << std::endl;
01615
01616 std::vector <double> efficiencies;
01617 efficiencies.assign(npar, 0.);
01618
01619 for (int ipar = 0; ipar < npar; ++ipar)
01620 for (int ichain = 0; ichain < nchains; ++ichain)
01621 {
01622 int index = ichain * npar + ipar;
01623 efficiencies[ipar] += double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index) + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.;
01624 }
01625
01626 for (int ipar = 0; ipar < npar; ++ipar)
01627 ofi
01628 << " (" << ipar << ") Parameter \""
01629 << fParameterSet->at(ipar)->GetName().data() << "\": "
01630 << efficiencies.at(ipar) << "%" << std::endl;
01631 ofi << std::endl;
01632 }
01633
01634 ofi
01635 << " -----------------------------------------------------" << std::endl
01636 << std::endl
01637 << " Notes" << std::endl
01638 << " =====" << std::endl
01639 << " (i) Median +- sigma denotes the median, m, of the" << std::endl
01640 << " marginalized distribution and the intervals from" << std::endl
01641 << " m to the 16% and 84% quantiles." << std::endl
01642 << " -----------------------------------------------------" << std::endl;
01643
01644
01645
01646 }
01647
01648
01649
01650 void BCModel::PrintShortFitSummary(int chi2flag)
01651 {
01652 BCLog::OutSummary("---------------------------------------------------");
01653 BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data()));
01654 BCLog::OutSummary(Form(" Number of parameters: Npar = %i", GetNParameters()));
01655 BCLog::OutSummary(Form(" Number of data points: Ndata = %i", GetNDataPoints()));
01656 BCLog::OutSummary(" Number of degrees of freedom:");
01657 BCLog::OutSummary(Form(" NDoF = Ndata - Npar = %i", GetNDataPoints()-GetNParameters()));
01658
01659 BCLog::OutSummary(" Best fit parameters (global):");
01660 for (unsigned int i = 0; i < GetNParameters(); ++i)
01661 BCLog::OutSummary(Form(" %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i)));
01662
01663 BCLog::OutSummary(" Goodness-of-fit test:");
01664 BCLog::OutSummary(Form(" p-value = %.3g", GetPValue()));
01665 if(chi2flag)
01666 {
01667 BCLog::OutSummary(Form(" p-value corrected for NDoF = %.3g", GetPValueNDoF()));
01668 BCLog::OutSummary(Form(" chi2 / NDoF = %.3g", GetChi2NDoF()));
01669 }
01670 BCLog::OutSummary("---------------------------------------------------");
01671 }
01672
01673
01674
01675 void BCModel::PrintHessianMatrix(std::vector<double> parameters)
01676 {
01677
01678 if (parameters.size() != GetNParameters())
01679 {
01680 BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
01681 return;
01682 }
01683
01684
01685 std::cout
01686 << std::endl
01687 << " Hessian matrix elements: " << std::endl
01688 << " Point: ";
01689
01690 for (int i = 0; i < int(parameters.size()); i++)
01691 std::cout << parameters.at(i) << " ";
01692 std::cout << std::endl;
01693
01694
01695 for (unsigned int i = 0; i < GetNParameters(); i++)
01696 for (unsigned int j = 0; j < i; j++)
01697 {
01698
01699 double hessianmatrixelement = HessianMatrixElement(fParameterSet->at(i),
01700 fParameterSet->at(j), parameters);
01701
01702
01703 std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl;
01704 }
01705 }
01706
01707
01708
01709 BCDataPoint * BCModel::VectorToDataPoint(std::vector<double> data)
01710 {
01711 int sizeofvector = int(data.size());
01712 BCDataPoint * datapoint = new BCDataPoint(sizeofvector);
01713 datapoint -> SetValues(data);
01714 return datapoint;
01715 }
01716
01717
01718
01719 int BCModel::CompareStrings(const char * string1, const char * string2)
01720 {
01721 int flag_same = 0;
01722
01723 if (strlen(string1) != strlen(string2))
01724 return -1;
01725
01726 for (int i = 0; i < int(strlen(string1)); i++)
01727 if (string1[i] != string2[i])
01728 flag_same = -1;
01729
01730 return flag_same;
01731 }
01732
01733
01734