00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BCModel.h"
00011
00012 #include "BCDataPoint.h"
00013 #include "BCDataSet.h"
00014 #include "BCParameter.h"
00015 #include "BCH1D.h"
00016 #include "BCH2D.h"
00017 #include "BCGoFTest.h"
00018 #include "BCLog.h"
00019 #include "BCErrorCodes.h"
00020 #include "BCMath.h"
00021 #include "BCModelOutput.h"
00022
00023 #include <TROOT.h>
00024 #include <TNamed.h>
00025 #include <TCanvas.h>
00026 #include <TPostScript.h>
00027 #include <TDirectory.h>
00028 #include <TFile.h>
00029 #include <TKey.h>
00030 #include <TTree.h>
00031 #include <TMath.h>
00032 #include <TGraph.h>
00033 #include <TH1.h>
00034 #include <TH2D.h>
00035 #include <TH1I.h>
00036 #include <TRandom3.h>
00037 #include <TF1.h>
00038 #include <Math/QuantFuncMathCore.h>
00039
00040 #include <algorithm>
00041 #include <set>
00042
00043 #include <fstream>
00044 #include <iomanip>
00045
00046
00047 BCModel::BCModel(const char * name)
00048 : BCIntegrate()
00049 {
00050 fNormalization = -1.;
00051 fDataSet = 0;
00052 fParameterSet = new BCParameterSet;
00053
00054 fIndex = -1;
00055 fPValue = -1;
00056 fPValueNDoF = -1;
00057 fChi2NDoF = -1;
00058
00059 fName = (char *) name;
00060 flag_ConditionalProbabilityEntry = true;
00061
00062 fDataPointUpperBoundaries = 0;
00063 fDataPointLowerBoundaries = 0;
00064
00065 fErrorBandXY = 0;
00066
00067 fGoFNChains = 5;
00068 fGoFNIterationsMax = 100000;
00069 fGoFNIterationsRun = 2000;
00070
00071 flag_discrete = false;
00072
00073 fPriorConstantAll = false;
00074 fPriorConstantValue = 0;
00075 }
00076
00077
00078 BCModel::BCModel()
00079 : BCIntegrate()
00080 {
00081 fNormalization = -1.;
00082 fDataSet = 0;
00083 fParameterSet = new BCParameterSet();
00084
00085 fIndex = -1;
00086 fPValue = -1;
00087 fPValueNDoF = -1;
00088 fChi2NDoF = -1;
00089
00090 fName = "model";
00091 fDataPointUpperBoundaries = 0;
00092 fDataPointLowerBoundaries = 0;
00093
00094 flag_ConditionalProbabilityEntry = true;
00095
00096 fGoFNChains = 5;
00097 fGoFNIterationsMax = 100000;
00098 fGoFNIterationsRun = 2000;
00099
00100 flag_discrete = false;
00101
00102 fPriorConstantAll = false;
00103 fPriorConstantValue = 0;
00104 }
00105
00106
00107 BCModel::BCModel(const BCModel & bcmodel) : BCIntegrate(bcmodel)
00108 {
00109 fIndex = bcmodel.fIndex;
00110 fName = bcmodel.fName;
00111 fModelAPriori = bcmodel.fModelAPriori;
00112 fModelAPosteriori = bcmodel.fModelAPosteriori;
00113 for (int i = 0; i < int(bcmodel.fParameterSet->size()); ++i) {
00114 if (bcmodel.fParameterSet->at(i)) {
00115 fParameterSet->push_back(new BCParameter(*(bcmodel.fParameterSet->at(i))));
00116
00117 }
00118 else
00119 fParameterSet->push_back(0);
00120 }
00121 if (fDataSet)
00122 fDataSet = bcmodel.fDataSet;
00123 else
00124 fDataSet = 0;
00125 if (bcmodel.fNDataPointsMinimum)
00126 fNDataPointsMinimum = bcmodel.fNDataPointsMinimum;
00127 else
00128 fNDataPointsMinimum = 0;
00129 if (bcmodel.fNDataPointsMaximum)
00130 fNDataPointsMaximum = bcmodel.fNDataPointsMaximum;
00131 else
00132 fNDataPointsMaximum = 0;
00133 flag_ConditionalProbabilityEntry = bcmodel.flag_ConditionalProbabilityEntry;
00134 fPValue = bcmodel.fPValue;
00135 fChi2NDoF = bcmodel.fChi2NDoF;
00136 fPValueNDoF = bcmodel.fPValueNDoF;
00137 flag_discrete = bcmodel.flag_discrete;
00138 fGoFNIterationsMax = bcmodel.fGoFNIterationsMax;
00139 fGoFNIterationsRun = bcmodel.fGoFNIterationsRun;
00140 fGoFNChains = bcmodel.fGoFNChains;
00141 for (int i = 0; i < int(bcmodel.fPriorContainer.size()); ++i) {
00142 if (bcmodel.fPriorContainer.at(i))
00143 fPriorContainer.push_back(new TNamed(*bcmodel.fPriorContainer.at(i)));
00144 else
00145 fPriorContainer.push_back(0);
00146 }
00147 fPriorConstantAll = bcmodel.fPriorConstantAll;
00148 fPriorConstantValue = bcmodel.fPriorConstantValue;
00149 fPriorContainerConstant = bcmodel.fPriorContainerConstant;
00150 fPriorContainerInterpolate = bcmodel.fPriorContainerInterpolate;
00151 fNormalization = bcmodel.fNormalization;
00152 }
00153
00154
00155 BCModel::~BCModel()
00156 {
00157 for (unsigned int i = 0; i < GetNParameters(); ++i)
00158 delete fPriorContainer[i];
00159 fPriorContainer.clear();
00160
00161 delete fParameterSet;
00162
00163 delete fDataPointLowerBoundaries;
00164 delete fDataPointUpperBoundaries;
00165 }
00166
00167
00168 BCModel & BCModel::operator = (const BCModel & bcmodel)
00169 {
00170 BCIntegrate::operator=(bcmodel);
00171 fIndex = bcmodel.fIndex;
00172 fName = bcmodel.fName;
00173 fModelAPriori = bcmodel.fModelAPriori;
00174 fModelAPosteriori = bcmodel.fModelAPosteriori;
00175 for (int i = 0; i < int(bcmodel.fParameterSet->size()); ++i) {
00176 if (bcmodel.fParameterSet->at(i)) {
00177 fParameterSet->push_back(new BCParameter(*(bcmodel.fParameterSet->at(i))));
00178
00179 }
00180 else
00181 fParameterSet->push_back(0);
00182 }
00183 if (fDataSet)
00184 fDataSet = bcmodel.fDataSet;
00185 else
00186 fDataSet = 0;
00187 if (bcmodel.fNDataPointsMinimum)
00188 fNDataPointsMinimum = bcmodel.fNDataPointsMinimum;
00189 else
00190 fNDataPointsMinimum = 0;
00191 if (bcmodel.fNDataPointsMaximum)
00192 fNDataPointsMaximum = bcmodel.fNDataPointsMaximum;
00193 else
00194 fNDataPointsMaximum = 0;
00195 flag_ConditionalProbabilityEntry = bcmodel.flag_ConditionalProbabilityEntry;
00196 fPValue = bcmodel.fPValue;
00197 fChi2NDoF = bcmodel.fChi2NDoF;
00198 fPValueNDoF = bcmodel.fPValueNDoF;
00199 flag_discrete = bcmodel.flag_discrete;
00200 fGoFNIterationsMax = bcmodel.fGoFNIterationsMax;
00201 fGoFNIterationsRun = bcmodel.fGoFNIterationsRun;
00202 fGoFNChains = bcmodel.fGoFNChains;
00203 for (int i = 0; i < int(bcmodel.fPriorContainer.size()); ++i) {
00204 if (bcmodel.fPriorContainer.at(i))
00205 fPriorContainer.push_back(new TNamed(*bcmodel.fPriorContainer.at(i)));
00206 else
00207 fPriorContainer.push_back(0);
00208 }
00209 fPriorConstantAll = bcmodel.fPriorConstantAll;
00210 fPriorConstantValue = bcmodel.fPriorConstantValue;
00211 fPriorContainerConstant = bcmodel.fPriorContainerConstant;
00212 fPriorContainerInterpolate = bcmodel.fPriorContainerInterpolate;
00213 fNormalization = bcmodel.fNormalization;
00214
00215
00216 return *this;
00217 }
00218
00219
00220 int BCModel::GetNDataPoints()
00221 {
00222 int npoints = 0;
00223 if (fDataSet)
00224 npoints = fDataSet->GetNDataPoints();
00225 else {
00226 BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined.");
00227 return ERROR_NOEVENTS;
00228 }
00229
00230 return npoints;
00231 }
00232
00233
00234 BCDataPoint * BCModel::GetDataPoint(unsigned int index)
00235 {
00236 if (fDataSet)
00237 return fDataSet->GetDataPoint(index);
00238
00239 BCLog::OutWarning("BCModel::GetDataPoint : No data set defined.");
00240 return 0;
00241 }
00242
00243
00244 BCParameter * BCModel::GetParameter(int index)
00245 {
00246 if (!fParameterSet)
00247 return 0;
00248
00249 if (index < 0 || index >= (int)GetNParameters()) {
00250 BCLog::OutWarning(
00251 Form("BCModel::GetParameter : Parameter index %d not within range.",index));
00252 return 0;
00253 }
00254
00255 return fParameterSet->at(index);
00256 }
00257
00258
00259 BCParameter * BCModel::GetParameter(const char * name)
00260 {
00261 if (!fParameterSet)
00262 return 0;
00263
00264 int index = -1;
00265 for (unsigned int i = 0; i < GetNParameters(); i++)
00266 if (name == GetParameter(i)->GetName())
00267 index = i;
00268
00269 if (index < 0) {
00270 BCLog::OutWarning(
00271 Form("BCModel::GetParameter : Model %s has no parameter named '%s'",
00272 GetName().data(), name));
00273 return 0;
00274 }
00275
00276 return GetParameter(index);
00277 }
00278
00279
00280 double BCModel::GetBestFitParameter(unsigned int index)
00281 {
00282 if(index >= GetNParameters()) {
00283 BCLog::OutError("BCModel::GetBestFitParameter : Parameter index out of range, returning -1e+111.");
00284 return -1e+111;
00285 }
00286
00287 if(fBestFitParameters.size()==0) {
00288 BCLog::OutError("BCModel::GetBestFitParameter : Mode finding not yet run, returning center of the range.");
00289 return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
00290 }
00291
00292 return fBestFitParameters[index];
00293 }
00294
00295
00296 double BCModel::GetBestFitParameterError(unsigned int index)
00297 {
00298 if(index >= GetNParameters()) {
00299 BCLog::OutError("BCModel::GetBestFitParameterError : Parameter index out of range, returning -1.");
00300 return -1;
00301 }
00302
00303 if(fBestFitParameterErrors.size()==0) {
00304 BCLog::OutError("BCModel::GetBestFitParameterError : Mode finding not yet run, returning -1.");
00305 return -1.;
00306 }
00307
00308 if(fBestFitParameterErrors[index]<0.)
00309 BCLog::OutWarning("BCModel::GetBestFitParameterError : Parameter error not available, returning -1.");
00310
00311 return fBestFitParameterErrors[index];
00312 }
00313
00314
00315 double BCModel::GetBestFitParameterMarginalized(unsigned int index)
00316 {
00317 if(index >= GetNParameters()) {
00318 BCLog::OutError("BCModel::GetBestFitParameterMarginalized : Parameter index out of range, returning -1e+111.");
00319 return -1e+111;
00320 }
00321
00322 if(fBestFitParametersMarginalized.size()==0) {
00323 BCLog::OutError("BCModel::GetBestFitParameterMarginalized : MCMC not yet run, returning center of the range.");
00324 return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
00325 }
00326
00327 return fBestFitParametersMarginalized[index];
00328 }
00329
00330
00331 void BCModel::SetNbins(const char * parname, int nbins)
00332 {
00333 BCParameter * p = GetParameter(parname);
00334 if (!p) {
00335 BCLog::OutWarning(
00336 Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname));
00337 return;
00338 }
00339
00340 BCIntegrate::SetNbins(nbins, p->GetIndex());
00341 }
00342
00343
00344 std::vector<double> BCModel::GetErrorBand(double level)
00345 {
00346 std::vector<double> errorband;
00347
00348 if (!fErrorBandXY)
00349 return errorband;
00350
00351 int nx = fErrorBandXY->GetNbinsX();
00352 errorband.assign(nx, 0.);
00353
00354
00355 for (int ix = 1; ix <= nx; ix++) {
00356 TH1D * temphist = fErrorBandXY->ProjectionY("temphist", ix, ix);
00357
00358 int nprobSum = 1;
00359 double q[1];
00360 double probSum[1];
00361 probSum[0] = level;
00362
00363 temphist->GetQuantiles(nprobSum, q, probSum);
00364
00365 errorband[ix - 1] = q[0];
00366 }
00367
00368 return errorband;
00369 }
00370
00371
00372 TGraph * BCModel::GetErrorBandGraph(double level1, double level2)
00373 {
00374 if (!fErrorBandXY)
00375 return 0;
00376
00377
00378 int nx = fErrorBandXY->GetNbinsX();
00379
00380 TGraph * graph = new TGraph(2 * nx);
00381 graph->SetFillStyle(1001);
00382 graph->SetFillColor(kYellow);
00383
00384
00385 std::vector<double> ymin = GetErrorBand(level1);
00386 std::vector<double> ymax = GetErrorBand(level2);
00387
00388 for (int i = 0; i < nx; i++) {
00389 graph->SetPoint(i, fErrorBandXY->GetXaxis()->GetBinCenter(i + 1), ymin[i]);
00390 graph->SetPoint(nx + i, fErrorBandXY->GetXaxis()->GetBinCenter(nx - i), ymax[nx - i - 1]);
00391 }
00392
00393 return graph;
00394 }
00395
00396
00397 TH2D * BCModel::GetErrorBandXY_yellow(double level, int nsmooth)
00398 {
00399 if (!fErrorBandXY)
00400 return 0;
00401
00402 int nx = fErrorBandXY->GetNbinsX();
00403 int ny = fErrorBandXY->GetNbinsY();
00404
00405
00406 TH2D * hist_tempxy = (TH2D*) fErrorBandXY->Clone(
00407 TString::Format("%s_sub_%f.2", fErrorBandXY->GetName(), level));
00408 hist_tempxy->Reset();
00409 hist_tempxy->SetFillColor(kYellow);
00410
00411
00412 for (int ix = 1; ix < nx; ix++) {
00413 BCH1D * hist_temp = new BCH1D();
00414
00415 TH1D * hproj = fErrorBandXY->ProjectionY("temphist", ix, ix);
00416 if (nsmooth > 0)
00417 hproj->Smooth(nsmooth);
00418
00419 hist_temp->SetHistogram(hproj);
00420
00421 TH1D * hist_temp_yellow = hist_temp->GetSmallestIntervalHistogram(level);
00422
00423 for (int iy = 1; iy <= ny; ++iy)
00424 hist_tempxy->SetBinContent(ix, iy, hist_temp_yellow->GetBinContent(iy));
00425
00426 delete hist_temp_yellow;
00427 delete hist_temp;
00428 }
00429
00430 return hist_tempxy;
00431 }
00432
00433
00434 TGraph * BCModel::GetFitFunctionGraph(std::vector<double> parameters)
00435 {
00436 if (!fErrorBandXY)
00437 return 0;
00438
00439
00440 int nx = fErrorBandXY->GetNbinsX();
00441 TGraph * graph = new TGraph(nx);
00442
00443
00444 for (int i = 0; i < nx; i++) {
00445 double x = fErrorBandXY->GetXaxis()->GetBinCenter(i + 1);
00446
00447 std::vector<double> xvec;
00448 xvec.push_back(x);
00449 double y = FitFunction(xvec, parameters);
00450 xvec.clear();
00451
00452 graph->SetPoint(i, x, y);
00453 }
00454
00455 return graph;
00456 }
00457
00458
00459 TGraph * BCModel::GetFitFunctionGraph(std::vector<double> parameters, double xmin, double xmax, int n)
00460 {
00461
00462 TGraph * graph = new TGraph(n + 1);
00463
00464 double dx = (xmax - xmin) / (double) n;
00465
00466
00467 for (int i = 0; i <= n; i++) {
00468 double x = (double) i * dx + xmin;
00469 std::vector<double> xvec;
00470 xvec.push_back(x);
00471 double y = FitFunction(xvec, parameters);
00472
00473 xvec.clear();
00474
00475 graph->SetPoint(i, x, y);
00476 }
00477
00478 return graph;
00479 }
00480
00481
00482 bool BCModel::GetFlagBoundaries()
00483 {
00484 if (!fDataPointLowerBoundaries)
00485 return false;
00486
00487 if (!fDataPointUpperBoundaries)
00488 return false;
00489
00490 if (fDataPointLowerBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
00491 return false;
00492
00493 if (fDataPointUpperBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
00494 return false;
00495
00496 return true;
00497 }
00498
00499
00500 void BCModel::SetSingleDataPoint(BCDataPoint * datapoint)
00501 {
00502
00503 BCDataSet * dataset = new BCDataSet();
00504
00505
00506 dataset->AddDataPoint(datapoint);
00507
00508
00509 SetDataSet(dataset);
00510 }
00511
00512
00513 void BCModel::SetSingleDataPoint(BCDataSet * dataset, unsigned int index)
00514 {
00515 if (index < 0 || index > dataset->GetNDataPoints())
00516 return;
00517
00518 SetSingleDataPoint(dataset->GetDataPoint(index));
00519 }
00520
00521
00522 void BCModel::SetDataBoundaries(unsigned int index, double lowerboundary, double upperboundary, bool fixed)
00523 {
00524
00525 if (!fDataSet) {
00526 BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
00527 return;
00528 }
00529
00530
00531 if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
00532 BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
00533 return;
00534 }
00535
00536
00537 if (!fDataPointLowerBoundaries)
00538 fDataPointLowerBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());
00539
00540 if (!fDataPointUpperBoundaries)
00541 fDataPointUpperBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());
00542
00543 if (fDataFixedValues.size() == 0)
00544 fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(), false);
00545
00546
00547 fDataPointLowerBoundaries->SetValue(index, lowerboundary);
00548 fDataPointUpperBoundaries->SetValue(index, upperboundary);
00549 fDataFixedValues[index] = fixed;
00550 }
00551
00552
00553 void BCModel::SetErrorBandContinuous(bool flag)
00554 {
00555 fErrorBandContinuous = flag;
00556
00557 if (flag)
00558 return;
00559
00560
00561 fErrorBandX.clear();
00562
00563
00564 for (unsigned int i = 0; i < fDataSet->GetNDataPoints(); ++i)
00565 fErrorBandX.push_back(fDataSet->GetDataPoint(i)->GetValue(fFitFunctionIndexX));
00566 }
00567
00568
00569 int BCModel::AddParameter(const char * name, double lowerlimit, double upperlimit)
00570 {
00571
00572 BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit);
00573
00574 int flag_ok = AddParameter(parameter);
00575 if (flag_ok)
00576 delete parameter;
00577
00578 return flag_ok;
00579 }
00580
00581
00582 int BCModel::AddParameter(BCParameter * parameter)
00583 {
00584
00585 if (!fParameterSet) {
00586 BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
00587 return ERROR_PARAMETERSETDOESNOTEXIST;
00588 }
00589
00590
00591 int flag_exists = 0;
00592 for (unsigned int i = 0; i < GetNParameters(); i++)
00593 if (CompareStrings(parameter->GetName().data(), GetParameter(i)->GetName().data()) == 0)
00594 flag_exists = -1;
00595
00596 if (flag_exists < 0) {
00597 BCLog::OutError(Form(
00598 "BCModel::AddParameter : Parameter with name %s exists already. ",
00599 parameter->GetName().data()));
00600 return ERROR_PARAMETEREXISTSALREADY;
00601 }
00602
00603
00604 unsigned int index = fParameterSet->size();
00605 parameter->SetIndex(index);
00606
00607
00608 fParameterSet->push_back(parameter);
00609
00610
00611 fPriorContainer.push_back(0);
00612
00613
00614 fPriorContainerInterpolate.push_back(false);
00615
00616
00617 fPriorContainerConstant.push_back(false);
00618
00619
00620 SetParameters(fParameterSet);
00621
00622
00623 ResetResults();
00624
00625 return 0;
00626 }
00627
00628
00629 double BCModel::LogProbabilityNN(std::vector<double> parameters)
00630 {
00631
00632 double logprob = LogLikelihood(parameters);
00633
00634
00635 if(fPriorConstantAll)
00636 logprob += fPriorConstantValue;
00637 else
00638 logprob += LogAPrioriProbability(parameters);
00639
00640 return logprob;
00641 }
00642
00643
00644 double BCModel::LogProbability(std::vector<double> parameters)
00645 {
00646
00647 if (fNormalization < 0. || fNormalization == 0.) {
00648 BCLog::OutError("BCModel::LogProbability. Normalization not available or zero.");
00649 return 0.;
00650 }
00651
00652 return LogProbabilityNN(parameters) - log(fNormalization);
00653 }
00654
00655
00656 double BCModel::LogLikelihood(std::vector<double> parameters)
00657 {
00658 double logprob = 0.;
00659
00660
00661 for (unsigned int i = 0; i < fDataSet->GetNDataPoints(); i++) {
00662 BCDataPoint * datapoint = GetDataPoint(i);
00663 logprob += LogConditionalProbabilityEntry(datapoint, parameters);
00664 }
00665
00666 return logprob;
00667 }
00668
00669
00670 double BCModel::LogAPrioriProbability(std::vector<double> parameters)
00671 {
00672 double logprob = 0;
00673
00674
00675 int npar = GetNParameters();
00676
00677
00678 for (int i = 0; i < npar; ++i) {
00679 if (fPriorContainer[i]) {
00680
00681 TF1 * f = dynamic_cast<TF1*>(fPriorContainer[i]);
00682 TH1 * h = dynamic_cast<TH1*>(fPriorContainer[i]);
00683
00684 if (f)
00685 logprob += log(f->Eval(parameters[i]));
00686 else if (h) {
00687 if(fPriorContainerInterpolate[i])
00688 logprob += log(h->Interpolate(parameters[i]));
00689 else
00690 logprob += log(h->GetBinContent(h->FindBin(parameters[i])));
00691 }
00692 else
00693 BCLog::OutError(Form(
00694 "BCModel::LogAPrioriProbability : Prior for parameter %s "
00695 "is defined but not recodnized.",
00696 GetParameter(i)->GetName().c_str()));
00697 }
00698
00699 else if (!fPriorContainerConstant[i]) {
00700 BCLog::OutWarning(Form(
00701 "BCModel::LogAPrioriProbability: Prior for parameter %s "
00702 "is undefined. Using constant prior to proceed.",
00703 GetParameter(i)->GetName().c_str()));
00704 logprob -= log(GetParameter(i)->GetRangeWidth());
00705 }
00706 }
00707
00708
00709 logprob += fPriorConstantValue;
00710
00711 return logprob;
00712 }
00713
00714
00715 double BCModel::LogEval(std::vector<double> parameters)
00716 {
00717 return LogProbabilityNN(parameters);
00718 }
00719
00720
00721 double BCModel::EvalSampling(std::vector<double> parameters)
00722 {
00723 return SamplingFunction(parameters);
00724 }
00725
00726
00727 double BCModel::SamplingFunction(std::vector<double> parameters)
00728 {
00729 double probability = 1.;
00730 for (std::vector<BCParameter *>::const_iterator it = fParameterSet->begin(); it != fParameterSet->end(); ++it)
00731 probability *= 1. / ((*it)->GetUpperLimit() - (*it)->GetLowerLimit());
00732 return probability;
00733 }
00734
00735
00736 double BCModel::Normalize()
00737 {
00738 if(fParameterSet->size()<1) {
00739 BCLog::OutError(Form("Normalize : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
00740 return -1.;
00741 }
00742
00743 BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",GetName().data()));
00744
00745 unsigned int n = GetNvar();
00746
00747
00748 if (n == 0) {
00749 SetParameters(fParameterSet);
00750 n = GetNvar();
00751 }
00752
00753
00754
00755 fNormalization = Integrate();
00756
00757 BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization));
00758
00759 return fNormalization;
00760 }
00761
00762
00763 int BCModel::CheckParameters(std::vector<double> parameters)
00764 {
00765
00766 if (!parameters.size() == fParameterSet->size())
00767 return ERROR_INVALIDNUMBEROFPARAMETERS;
00768
00769
00770 for (unsigned int i = 0; i < fParameterSet->size(); i++) {
00771 BCParameter * modelparameter = fParameterSet->at(i);
00772
00773 if (modelparameter->GetLowerLimit() > parameters.at(i)
00774 || modelparameter->GetUpperLimit() < parameters.at(i)) {
00775 BCLog::OutError(Form(
00776 "BCModel::CheckParameters : Parameter %s not within limits.",
00777 fParameterSet-> at(i)-> GetName().data()));
00778 return ERROR_PARAMETERNOTWITHINRANGE;
00779 }
00780 }
00781
00782 return 0;
00783 }
00784
00785
00786 void BCModel::FindMode(std::vector<double> start)
00787 {
00788 if(fParameterSet->size()<1) {
00789 BCLog::OutError(Form("FindMode : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
00790 return;
00791 }
00792
00793
00794
00795 BCLog::OutSummary(Form("Model \'%s\': Finding mode using %s",GetName().data(), DumpOptimizationMethod().c_str()));
00796
00797
00798 SetParameters(fParameterSet);
00799
00800 switch (GetOptimizationMethod()) {
00801 case BCIntegrate::kOptSA:
00802 FindModeSA(start);
00803 return;
00804
00805 case BCIntegrate::kOptMinuit: {
00806 int printlevel = -1;
00807 if (BCLog::GetLogLevelScreen() <= BCLog::detail)
00808 printlevel = 0;
00809 if (BCLog::GetLogLevelScreen() <= BCLog::debug)
00810 printlevel = 1;
00811 BCIntegrate::FindModeMinuit(start, printlevel);
00812 return;
00813 }
00814
00815 case BCIntegrate::kOptMetropolis:
00816 MarginalizeAll();
00817 return;
00818
00819 default:
00820 BCLog::OutError(Form("BCModel::FindMode : Invalid mode finding method: %d",GetOptimizationMethod()));
00821 }
00822
00823
00824 return;
00825 }
00826
00827
00828 void BCModel::FindModeMinuit(std::vector<double> start, int printlevel)
00829 {
00830 if(fParameterSet->size()<1) {
00831 BCLog::OutError(Form("FindModeMinuit : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
00832 return;
00833 }
00834
00835
00836 SetParameters(fParameterSet);
00837 BCIntegrate::FindModeMinuit(start, printlevel);
00838 }
00839
00840
00841 void BCModel::WriteMode(const char * file)
00842 {
00843 ofstream ofi(file);
00844 if (!ofi.is_open()) {
00845 std::cerr << "Couldn't open file " << file << std::endl;
00846 return;
00847 }
00848
00849 int npar = fParameterSet->size();
00850 for (int i = 0; i < npar; i++)
00851 ofi << fBestFitParameters.at(i) << std::endl;
00852
00853 ofi << std::endl;
00854 ofi << "#######################################################################"
00855 << std::endl;
00856 ofi << "#" << std::endl;
00857 ofi << "# This file was created automatically by BCModel::WriteMode() call."
00858 << std::endl;
00859 ofi << "# It can be read in by call to BCModel::ReadMode()." << std::endl;
00860 ofi << "# Do not modify it unless you know what you're doing." << std::endl;
00861 ofi << "#" << std::endl;
00862 ofi << "#######################################################################"
00863 << std::endl;
00864 ofi << "#" << std::endl;
00865 ofi << "# Best fit parameters (mode) for model:" << std::endl;
00866 ofi << "# \'" << fName.data() << "\'" << std::endl;
00867 ofi << "#" << std::endl;
00868 ofi << "# Number of parameters: " << npar << std::endl;
00869 ofi << "# Parameters ordered as above:" << std::endl;
00870
00871 for (int i = 0; i < npar; i++) {
00872 ofi << "# " << i << ": ";
00873 ofi << fParameterSet->at(i)->GetName().data() << " = ";
00874 ofi << fBestFitParameters.at(i) << std::endl;
00875 }
00876
00877 ofi << "#" << std::endl;
00878 ofi << "########################################################################"
00879 << std::endl;
00880 }
00881
00882
00883 int BCModel::ReadMode(const char * file)
00884 {
00885 ifstream ifi(file);
00886 if (!ifi.is_open()) {
00887 BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.", file));
00888 return 0;
00889 }
00890
00891 int npar = fParameterSet->size();
00892 std::vector<double> mode;
00893
00894 int i = 0;
00895 while (i < npar && !ifi.eof()) {
00896 double a;
00897 ifi >> a;
00898 mode.push_back(a);
00899 i++;
00900 }
00901
00902 if (i < npar) {
00903 BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.", file));
00904 BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.", npar, i));
00905 return 0;
00906 }
00907
00908 BCLog::OutSummary(
00909 Form("# Read in best fit parameters (mode) for model \'%s\' from file %s:",
00910 fName.data(), file));
00911 SetMode(mode);
00912 for (int j = 0; j < npar; j++)
00913 BCLog::OutSummary(Form("# ->Parameter %d : %s = %e", j,
00914 fParameterSet->at(j)->GetName().data(), fBestFitParameters[j]));
00915
00916 BCLog::OutWarning("# ! Best fit values obtained before this call will be overwritten !");
00917
00918 return npar;
00919 }
00920
00921
00922 int BCModel::MarginalizeAll()
00923 {
00924 if(fParameterSet->size()<1) {
00925 BCLog::OutError(Form("MarginalizeAll : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
00926 return 0;
00927 }
00928
00929 BCLog::OutSummary(Form("Running MCMC for model \'%s\'",GetName().data()));
00930
00931
00932 double dx = 0.;
00933 double dy = 0.;
00934
00935 if (fFitFunctionIndexX >= 0) {
00936 dx = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX)
00937 - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX))
00938 / (double) fErrorBandNbinsX;
00939
00940 dy = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY)
00941 - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY))
00942 / (double) fErrorBandNbinsY;
00943
00944 fErrorBandXY
00945 = new TH2D(TString::Format("errorbandxy_%d", BCLog::GetHIndex()), "",
00946 fErrorBandNbinsX,
00947 fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX) - .5 * dx,
00948 fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX) + .5 * dx,
00949 fErrorBandNbinsY,
00950 fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY) - .5 * dy,
00951 fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY) + .5 * dy);
00952 fErrorBandXY->SetStats(kFALSE);
00953
00954 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
00955 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
00956 fErrorBandXY->SetBinContent(ix, iy, 0.);
00957 }
00958
00959
00960 MCMCMetropolis();
00961
00962
00963 StoreMode();
00964
00965 return 1;
00966 }
00967
00968
00969 BCH1D * BCModel::GetMarginalized(BCParameter * parameter)
00970 {
00971 if (!parameter) {
00972
00973
00974 return 0;
00975 }
00976
00977 int index = parameter->GetIndex();
00978 if (fMCMCFlagsFillHistograms[index] == false) {
00979
00980
00981
00982
00983 return 0;
00984 }
00985
00986
00987 TH1D * hist = MCMCGetH1Marginalized(index);
00988 if (!hist)
00989 return 0;
00990
00991
00992 hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data()));
00993 hist->SetXTitle(parameter->GetName().data());
00994 hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data()));
00995 hist->SetStats(kFALSE);
00996
00997
00998 BCH1D * hprob = new BCH1D();
00999 hprob->SetHistogram(hist);
01000
01001
01002 double bestfit = hprob->GetMode();
01003
01004 if (fBestFitParametersMarginalized.empty())
01005 fBestFitParametersMarginalized.assign(GetNParameters(), 0.0);
01006
01007 fBestFitParametersMarginalized[index] = bestfit;
01008
01009 hprob->SetGlobalMode(fBestFitParameters.at(index));
01010
01011 return hprob;
01012 }
01013
01014
01015 int BCModel::ReadMarginalizedFromFile(const char * file)
01016 {
01017 TFile * froot = new TFile(file);
01018 if (!froot->IsOpen()) {
01019 BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.", file));
01020 return 0;
01021 }
01022
01023
01024
01025
01026
01027
01028
01029 MCMCInitialize();
01030
01031 int k = 0;
01032 int n = GetNParameters();
01033 for (int i = 0; i < n; i++) {
01034 BCParameter * a = GetParameter(i);
01035 TKey * key = froot->GetKey(TString::Format("hist_%s_%s",GetName().data(), a->GetName().data()));
01036 if (key) {
01037 TH1D * h1 = (TH1D*) key->ReadObjectAny(TH1D::Class());
01038 h1->SetDirectory(0);
01039 if (SetMarginalized(i, h1))
01040 k++;
01041 }
01042 else
01043 BCLog::OutWarning(
01044 Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.",
01045 GetName().data(), a->GetName().data(), file));
01046 }
01047
01048 for (int i = 0; i < n - 1; i++) {
01049 for (int j = i + 1; j < n; j++) {
01050 BCParameter * a = GetParameter(i);
01051 BCParameter * b = GetParameter(j);
01052 TKey * key = froot->GetKey(
01053 TString::Format("hist_%s_%s_%s",GetName().data(), a->GetName().data(),b->GetName().data()));
01054 if (key) {
01055 TH2D * h2 = (TH2D*) key->ReadObjectAny(TH2D::Class());
01056 h2->SetDirectory(0);
01057 if (SetMarginalized(i, j, h2))
01058 k++;
01059 }
01060 else
01061 BCLog::OutWarning(
01062 Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.",
01063 GetName().data(), a->GetName().data(), b->GetName().data(), file));
01064 }
01065 }
01066
01067 froot->Close();
01068
01069 return k;
01070 }
01071
01072
01073 int BCModel::ReadErrorBandFromFile(const char * file)
01074 {
01075 TFile * froot = new TFile(file);
01076 if (!froot->IsOpen()) {
01077 BCLog::OutError(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.", file));
01078 return 0;
01079 }
01080
01081 int r = 0;
01082
01083 TH2D * h2 = (TH2D*) froot->Get("errorbandxy");
01084 if (h2) {
01085 h2->SetDirectory(0);
01086 h2->SetName(TString::Format("errorbandxy_%d", BCLog::GetHIndex()));
01087 SetErrorBandHisto(h2);
01088 r = 1;
01089 }
01090 else
01091 BCLog::OutWarning(
01092 Form("BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",file));
01093
01094 froot->Close();
01095
01096 return r;
01097 }
01098
01099
01100 int BCModel::PrintAllMarginalized1D(const char * filebase)
01101 {
01102 if (fMCMCH1Marginalized.size() == 0) {
01103 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
01104 return 0;
01105 }
01106
01107 int n = GetNParameters();
01108 for (int i = 0; i < n; i++) {
01109 BCParameter * a = GetParameter(i);
01110 if (GetMarginalized(a))
01111 GetMarginalized(a)->Print(Form("%s_1D_%s.ps", filebase, a->GetName().data()));
01112 }
01113
01114 return n;
01115 }
01116
01117
01118 int BCModel::PrintAllMarginalized2D(const char * filebase)
01119 {
01120 if (fMCMCH2Marginalized.size() == 0) {
01121 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
01122 return 0;
01123 }
01124
01125 int k = 0;
01126 int n = GetNParameters();
01127 for (int i = 0; i < n - 1; i++) {
01128 for (int j = i + 1; j < n; j++) {
01129 BCParameter * a = GetParameter(i);
01130 BCParameter * b = GetParameter(j);
01131
01132 double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
01133 double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
01134 if (deltaa <= 1e-7 * meana)
01135 continue;
01136
01137 double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
01138 double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
01139 if (deltab <= 1e-7 * meanb)
01140 continue;
01141
01142 if (GetMarginalized(a, b))
01143 GetMarginalized(a, b)->Print(
01144 Form("%s_2D_%s_%s.ps",filebase, a->GetName().data(), b->GetName().data()));
01145 k++;
01146 }
01147 }
01148
01149 return k;
01150 }
01151
01152
01153 int BCModel::PrintAllMarginalized(const char * file, unsigned int hdiv, unsigned int vdiv)
01154 {
01155 if (!fMCMCFlagFillHistograms) {
01156 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled.");
01157 return 0;
01158 }
01159
01160 int npar = GetNParameters();
01161
01162 if (fMCMCH1Marginalized.size() == 0 || (fMCMCH2Marginalized.size() == 0
01163 && npar > 1)) {
01164 BCLog::OutError(
01165 "BCModel::PrintAllMarginalized : Marginalized distributions not available.");
01166 return 0;
01167 }
01168
01169
01170 if (fMCMCH1Marginalized.size() == 1 && fMCMCH2Marginalized.size() == 0) {
01171 BCParameter * a = GetParameter(0);
01172 if (GetMarginalized(a))
01173 GetMarginalized(a)->Print(file);
01174 return 1;
01175 }
01176
01177 int c_width = 600;
01178 int c_height = 600;
01179
01180 int type = 112;
01181
01182 if (hdiv > vdiv) {
01183 if (hdiv > 3) {
01184 c_width = 1000;
01185 c_height = 700;
01186 }
01187 else {
01188 c_width = 800;
01189 c_height = 600;
01190 }
01191 }
01192 else if (hdiv < vdiv) {
01193 if (hdiv > 3) {
01194 c_height = 1000;
01195 c_width = 700;
01196 }
01197 else {
01198 c_height = 800;
01199 c_width = 600;
01200 }
01201 type = 111;
01202 }
01203
01204
01205 int nplots2d = npar * (npar - 1) / 2;
01206 int nplots = npar + nplots2d;
01207
01208
01209 BCLog::OutSummary(
01210 Form(
01211 "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
01212 npar, nplots2d, nplots, file));
01213 if (nplots > 100)
01214 BCLog::OutDetail("This can take a while...");
01215
01216
01217 TCanvas * c = new TCanvas("c", "canvas", c_width, c_height);
01218
01219 TPostScript * ps = new TPostScript(file, type);
01220
01221 if (type == 112)
01222 ps->Range(24, 16);
01223 else
01224 ps->Range(16, 24);
01225
01226
01227 ps->NewPage();
01228 c->cd();
01229 c->Clear();
01230 c->Divide(hdiv, vdiv);
01231
01232 int n = 0;
01233 for (int i = 0; i < npar; i++) {
01234
01235 BCParameter * a = GetParameter(i);
01236
01237
01238 if (!GetMarginalized(a))
01239 continue;
01240
01241
01242 if (GetMarginalized(a)->GetHistogram()->Integral() <= 0)
01243 continue;
01244
01245
01246 if (i != 0 && i % (hdiv * vdiv) == 0) {
01247 c->Update();
01248 ps->NewPage();
01249 c->cd();
01250 c->Clear();
01251 c->Divide(hdiv, vdiv);
01252 }
01253
01254
01255 c->cd(i % (hdiv * vdiv) + 1);
01256
01257 GetMarginalized(a)->Draw();
01258 n++;
01259
01260 if (n % 100 == 0)
01261 BCLog::OutDetail(Form(" --> %d plots done", n));
01262 }
01263
01264 if (n > 0) {
01265 c->Update();
01266
01267
01268 ps->NewPage();
01269 c->cd();
01270 c->Clear();
01271 }
01272
01273 c->Divide(hdiv, vdiv);
01274
01275 int k = 0;
01276 for (int i = 0; i < npar - 1; i++) {
01277 if (!fMCMCFlagsFillHistograms[i])
01278 continue;
01279 for (int j = i + 1; j < npar; j++) {
01280 if (!fMCMCFlagsFillHistograms[j])
01281 continue;
01282
01283
01284 BCParameter * a = GetParameter(i);
01285 BCParameter * b = GetParameter(j);
01286
01287
01288 if (!GetMarginalized(a, b))
01289 continue;
01290
01291
01292 if (GetMarginalized(a, b)->GetHistogram()->Integral() <= 0)
01293 continue;
01294
01295
01296 if (k != 0 && k % (hdiv * vdiv) == 0) {
01297 c->Update();
01298 ps->NewPage();
01299 c->cd();
01300 c->Clear();
01301 c->Divide(hdiv, vdiv);
01302 }
01303
01304
01305 c->cd(k % (hdiv * vdiv) + 1);
01306
01307 double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
01308 double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
01309 if (deltaa <= 1e-7 * meana)
01310 continue;
01311
01312 double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
01313 double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
01314 if (deltab <= 1e-7 * meanb)
01315 continue;
01316
01317 GetMarginalized(a, b)->Draw(52);
01318 k++;
01319
01320 if ((n + k) % 100 == 0)
01321 BCLog::OutDetail(Form(" --> %d plots done", n + k));
01322 }
01323 }
01324
01325 if ((n + k) > 100 && (n + k) % 100 != 0)
01326 BCLog::OutDetail(Form(" --> %d plots done", n + k));
01327
01328 c->Update();
01329 ps->Close();
01330
01331 delete c;
01332 delete ps;
01333
01334
01335 return n + k;
01336 }
01337
01338
01339 BCH2D * BCModel::GetMarginalized(BCParameter * par1, BCParameter * par2)
01340 {
01341 int index1 = par1->GetIndex();
01342 int index2 = par2->GetIndex();
01343
01344 if (fMCMCFlagsFillHistograms[index1] == false || fMCMCFlagsFillHistograms[index2] == false) {
01345
01346
01347
01348
01349 return 0;
01350 }
01351
01352 if (index1 == index2) {
01353
01354
01355 return 0;
01356 }
01357
01358 BCParameter * npar1 = par1;
01359 BCParameter * npar2 = par2;
01360
01361 if (index1 > index2) {
01362 npar1 = par2;
01363 npar2 = par1;
01364
01365 int itmp = index1;
01366 index1 = index2;
01367 index2 = itmp;
01368 }
01369
01370
01371 TH2D * hist = MCMCGetH2Marginalized(index1, index2);
01372
01373 if (hist == 0)
01374 return 0;
01375
01376 BCH2D * hprob = new BCH2D();
01377
01378
01379 hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data()));
01380 hist->SetXTitle(Form("%s", npar1->GetName().data()));
01381 hist->SetYTitle(Form("%s", npar2->GetName().data()));
01382 hist->SetStats(kFALSE);
01383
01384 double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
01385 hprob->SetGlobalMode(gmode);
01386
01387
01388 hprob->SetHistogram(hist);
01389
01390 return hprob;
01391 }
01392
01393
01394 double BCModel::GetPvalueFromChi2(std::vector<double> par, int sigma_index)
01395 {
01396 double ll = LogLikelihood(par);
01397 int n = GetNDataPoints();
01398
01399 double sum_sigma = 0;
01400 for (int i = 0; i < n; i++)
01401 sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));
01402
01403 double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);
01404
01405 fPValue = TMath::Prob(chi2, n);
01406
01407 return fPValue;
01408 }
01409
01410
01411 std::vector<double> BCModel::GetChi2Runs(int dataIndex, int sigmaIndex)
01412 {
01413 std::vector<double> x;
01414 return x;
01415 }
01416
01417
01418 double BCModel::GetPvalueFromChi2Johnson(std::vector<double> par)
01419 {
01420 double chi2 = GetChi2Johnson(par);
01421
01422 fPValue = TMath::Prob(chi2, NumberBins() - 1);
01423 return fPValue;
01424 }
01425
01426
01427 double BCModel::GetChi2Johnson(std::vector<double> par, int nBins)
01428 {
01429 typedef unsigned int uint;
01430
01431
01432 int n = GetNDataPoints();
01433
01434 if (nBins < 0)
01435 nBins = NumberBins();
01436
01437
01438 std::vector<double> a;
01439 for (int i = 0; i <= nBins; i++)
01440 a.push_back(i / double(nBins));
01441
01442
01443 TH1I * hist = new TH1I("h1", "h1 title", nBins, 0., 1.);
01444
01445
01446
01447 if (flag_discrete) {
01448
01449 for (int j = 0; j < n; j++) {
01450
01451 double CDFval = CDF(par, j, false);
01452
01453 double CDFlower = CDF(par, j, true);
01454
01455
01456 int qMax = 1;
01457 for (int i = 0; i < nBins; i++) {
01458 if (CDFval > a[i])
01459 qMax = i + 1;
01460 else
01461 break;
01462 }
01463 int qMin = 1;
01464 for (int i = 0; i < nBins; i++) {
01465 if (CDFlower > a[i])
01466 qMin = i + 1;
01467 else
01468 break;
01469 }
01470
01471
01472 if (qMin == qMax) {
01473
01474 if (CDFval < 1)
01475 hist->Fill(CDFval);
01476 else
01477 hist->AddBinContent(qMax);
01478
01479 continue;
01480 }
01481
01482
01483
01484
01485
01486
01487
01488 std::vector<double> prob;
01489
01490 double norm = 1 / double(CDFval - CDFlower);
01491
01492 for (int i = 0; i < (qMax - qMin + 1); i++) {
01493 if (i == 0) {
01494 prob.push_back(norm * (a[qMin] - CDFlower));
01495 continue;
01496 }
01497 if (i == (qMax - qMin)) {
01498 prob.push_back(norm * (CDFval - a[qMax - 1]));
01499 continue;
01500 }
01501
01502 prob.push_back(norm * (a[i] - a[i - 1]));
01503 }
01504
01505 double U = fRandom->Rndm();
01506
01507 for (uint i = 1; i < prob.size(); i++)
01508 prob[i] += prob[i - 1];
01509
01510 for (uint i = 0; i < prob.size(); i++) {
01511
01512 if (U <= prob[i]) {
01513 hist->Fill((a[qMin + i - 1] + a[qMin + i]) / 2.);
01514 break;
01515 }
01516 }
01517 }
01518 }
01519 else {
01520 for (int j = 0; j < n; j++)
01521 hist->Fill(CDF(par, j, false));
01522 }
01523
01524
01525 double chi2 = 0.;
01526 double mk, pk;
01527 double N = double(n);
01528 for (int i = 1; i <= nBins; i++) {
01529 mk = hist->GetBinContent(i);
01530 pk = a[i] - a[i - 1];
01531 chi2 += (mk - N * pk) * (mk - N * pk) / (N * pk);
01532 }
01533
01534 delete hist;
01535
01536 return chi2;
01537 }
01538
01539
01540 double BCModel::GetAvalueFromChi2Johnson(TTree * tree, TH1D * distribution)
01541 {
01542
01543 int nPar = (int)GetNParameters();
01544 std::vector<double> param(nPar);
01545
01546
01547 int nParBranches = -1;
01548 tree->SetBranchAddress("fNParameters", &nParBranches);
01549
01550
01551 tree->GetEntry(0);
01552 if (nParBranches != nPar) {
01553 BCLog::OutError(Form(
01554 "Cannot compute A value: number of parameters in tree (%d)"
01555 "doesn't match number of parameters in model (%d)",
01556 nParBranches, nPar));
01557 return -1.;
01558 }
01559
01560
01561 char * branchName = new char[10 + nPar];
01562
01563
01564 for (int i = 0; i < (int) nPar; i++) {
01565 +sprintf(branchName, "fParameter%d", i);
01566 tree->SetBranchAddress(branchName, ¶m[i]);
01567 }
01568
01569
01570 long nEntries = tree->GetEntries();
01571
01572
01573 std::vector<double> randoms(nEntries);
01574 int K = NumberBins();
01575 BCMath::RandomChi2(randoms, K - 1);
01576
01577
01578 int nBigger = 0;
01579 for (int i = 0; i < nEntries; i++) {
01580 tree->GetEntry(i);
01581 double chi2 = GetChi2Johnson(param);
01582 if (distribution != 0)
01583 distribution->Fill(chi2);
01584
01585 if (chi2 > randoms[i])
01586 nBigger++;
01587 }
01588
01589 return nBigger / double(nEntries);
01590 }
01591
01592
01593 double BCModel::GetPvalueFromChi2NDoF(std::vector<double> par, int sigma_index)
01594 {
01595 double ll = LogLikelihood(par);
01596 int n = GetNDataPoints();
01597 int npar = GetNParameters();
01598
01599 double sum_sigma = 0;
01600 for (int i = 0; i < n; i++)
01601 sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));
01602
01603 double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);
01604
01605 fChi2NDoF = chi2 / double(n - npar);
01606 fPValueNDoF = TMath::Prob(chi2, n - npar);
01607
01608 return fPValueNDoF;
01609 }
01610
01611
01612 double BCModel::GetPvalueFromKolmogorov(const std::vector<double>& par,int index)
01613 {
01614 if (flag_discrete) {
01615 BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
01616 "test defined only for continuous distributions."));
01617 return -1.;
01618 }
01619
01620
01621 std::vector<double> yData = fDataSet->GetDataComponents(index);
01622 TH1D * ECDF = BCMath::ECDF(yData);
01623
01624 int N = GetNDataPoints();
01625
01626
01627 std::set<double> uniqueObservations;
01628 for (int i = 0; i < N; i++)
01629 uniqueObservations.insert(CDF(par, i, false));
01630
01631 int nUnique = uniqueObservations.size();
01632 if (nUnique != ECDF->GetNbinsX() + 1) {
01633 BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
01634 "Number of unique data points doesn't match (%d vs %d)", nUnique,
01635 ECDF->GetNbinsX() + 1));
01636 return -1.;
01637 }
01638
01639
01640 double distMax = 0.;
01641
01642
01643 double dist = 0.;
01644
01645 std::set<double>::const_iterator iter = uniqueObservations.begin();
01646 for (int iBin = 0; iBin < nUnique; ++iBin) {
01647
01648 dist = TMath::Abs(*iter - ECDF->GetBinContent(iBin + 1));
01649
01650 distMax = TMath::Max(dist, distMax);
01651
01652
01653
01654
01655
01656
01657 ++iter;
01658 }
01659
01660
01661
01662 double z = distMax * sqrt(N);
01663
01664 fPValue = TMath::KolmogorovProb(z);
01665
01666
01667
01668
01669
01670 delete ECDF;
01671
01672 return fPValue;
01673 }
01674
01675
01676 BCH1D * BCModel::CalculatePValue(std::vector<double> par, bool flag_histogram)
01677 {
01678 BCH1D * hist = 0;
01679
01680
01681 BCLog::OutSummary("Do goodness-of-fit-test");
01682
01683
01684 BCGoFTest * goftest = new BCGoFTest("modeltest");
01685
01686
01687 goftest->SetTestModel(this);
01688
01689
01690
01691 if (!goftest->SetTestPoint(par))
01692 return 0;
01693
01694
01695
01696 goftest->MCMCSetFlagFillHistograms(false);
01697
01698
01699 goftest->MCMCSetNChains(fGoFNChains);
01700 goftest->MCMCSetNIterationsMax(fGoFNIterationsMax);
01701 goftest->MCMCSetNIterationsRun(fGoFNIterationsRun);
01702
01703
01704 fPValue = goftest->GetCalculatedPValue(flag_histogram);
01705
01706
01707 if (flag_histogram) {
01708 hist = new BCH1D();
01709 hist->SetHistogram(goftest->GetHistogramLogProb());
01710 }
01711
01712
01713 delete goftest;
01714
01715
01716 return hist;
01717 }
01718
01719
01720 void BCModel::CorrelateDataPointValues(std::vector<double> &x)
01721 {
01722
01723 }
01724
01725
01726 double BCModel::HessianMatrixElement(BCParameter * par1, BCParameter * par2, std::vector<double> point)
01727 {
01728
01729 if (point.size() != GetNParameters()) {
01730 BCLog::OutError("BCModel::HessianMatrixElement : Invalid number of entries in the vector.");
01731 return -1;
01732 }
01733
01734
01735 double nsteps = 1e5;
01736 double dx1 = par1->GetRangeWidth() / nsteps;
01737 double dx2 = par2->GetRangeWidth() / nsteps;
01738
01739
01740 std::vector<double> xpp = point;
01741 std::vector<double> xpm = point;
01742 std::vector<double> xmp = point;
01743 std::vector<double> xmm = point;
01744
01745 int idx1 = par1->GetIndex();
01746 int idx2 = par2->GetIndex();
01747
01748 xpp[idx1] += dx1;
01749 xpp[idx2] += dx2;
01750
01751 xpm[idx1] += dx1;
01752 xpm[idx2] -= dx2;
01753
01754 xmp[idx1] -= dx1;
01755 xmp[idx2] += dx2;
01756
01757 xmm[idx1] -= dx1;
01758 xmm[idx2] -= dx2;
01759
01760
01761 double ppp = Likelihood(xpp);
01762 double ppm = Likelihood(xpm);
01763 double pmp = Likelihood(xmp);
01764 double pmm = Likelihood(xmm);
01765
01766
01767 return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
01768 }
01769
01770
01771 void BCModel::FixDataAxis(unsigned int index, bool fixed)
01772 {
01773
01774 if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
01775 BCLog::OutError("BCModel::FixDataAxis : Index out of range.");
01776 return;
01777 }
01778
01779 if (fDataFixedValues.size() == 0)
01780 fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(),
01781 false);
01782
01783 fDataFixedValues[index] = fixed;
01784 }
01785
01786
01787 bool BCModel::GetFixedDataAxis(unsigned int index)
01788 {
01789
01790 if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
01791 BCLog::OutError("BCModel::GetFixedDataAxis : Index out of range.");
01792 return false;
01793 }
01794
01795 return fDataFixedValues[index];
01796 }
01797
01798
01799 int BCModel::SetPrior(int index, TF1 * f)
01800 {
01801
01802 if (index < 0 && index >= int(GetNParameters())) {
01803 BCLog::OutError("BCModel::SetPrior : Index out of range.");
01804 return 0;
01805 }
01806
01807 if (fPriorContainer[index])
01808 delete fPriorContainer[index];
01809
01810
01811 fPriorContainer[index] = new TF1(*f);
01812
01813 fPriorContainerConstant[index] = false;
01814
01815 RecalculatePriorConstant();
01816
01817
01818 ResetResults();
01819
01820
01821 return 1;
01822 }
01823
01824
01825 int BCModel::SetPrior(const char * name, TF1 * f)
01826 {
01827
01828 int index = -1;
01829 for (unsigned int i = 0; i < GetNParameters(); i++)
01830 if (name == GetParameter(i)->GetName())
01831 index = i;
01832
01833
01834 return SetPrior(index, f);
01835 }
01836
01837
01838 int BCModel::SetPriorDelta(int index, double value)
01839 {
01840
01841 SetParameterRange(index, value, value);
01842
01843
01844 return SetPriorConstant(index);
01845 }
01846
01847
01848 int BCModel::SetPriorDelta(const char* name, double value)
01849 {
01850
01851 int index = -1;
01852 for (unsigned int i = 0; i < GetNParameters(); i++)
01853 if (name == GetParameter(i)->GetName())
01854 index = i;
01855
01856
01857 return SetPriorDelta(index, value);
01858 }
01859
01860
01861 int BCModel::SetPriorGauss(int index, double mean, double sigma)
01862 {
01863
01864 if (index < 0 || index >= int(GetNParameters())) {
01865 BCLog::OutError("BCModel::SetPriorGauss : Index out of range.");
01866 return 0;
01867 }
01868
01869
01870 TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
01871 "1./sqrt(2.*TMath::Pi())/[1] * exp(- (x-[0])*(x-[0])/2./[1]/[1])",
01872 GetParameter(index)->GetLowerLimit(),
01873 GetParameter(index)->GetUpperLimit());
01874 f->SetParameter(0, mean);
01875 f->SetParameter(1, sigma);
01876
01877
01878 return SetPrior(index, f);
01879 }
01880
01881
01882 int BCModel::SetPriorGauss(const char* name, double mean, double sigma)
01883 {
01884
01885 int index = -1;
01886 for (unsigned int i = 0; i < GetNParameters(); i++)
01887 if (name == GetParameter(i)->GetName())
01888 index = i;
01889
01890
01891 return SetPriorGauss(index, mean, sigma);
01892 }
01893
01894
01895 int BCModel::SetPriorGauss(int index, double mean, double sigmadown, double sigmaup)
01896 {
01897
01898 if (index < 0 || index >= int(GetNParameters())) {
01899 BCLog::OutError("BCModel::SetPriorGauss : Index out of range.");
01900 return 0;
01901 }
01902
01903
01904 TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
01905 BCMath::SplitGaussian,
01906 GetParameter(index)->GetLowerLimit(),
01907 GetParameter(index)->GetUpperLimit(),
01908 3);
01909 f->SetParameter(0, mean);
01910 f->SetParameter(1, sigmadown);
01911 f->SetParameter(2, sigmaup);
01912
01913
01914 return SetPrior(index, f);
01915
01916 return 0;
01917 }
01918
01919
01920 int BCModel::SetPriorGauss(const char * name, double mean, double sigmadown, double sigmaup)
01921 {
01922
01923 int index = -1;
01924 for (unsigned int i = 0; i < GetNParameters(); i++)
01925 if (name == GetParameter(i)->GetName())
01926 index = i;
01927
01928
01929 return SetPriorGauss(index, mean, sigmadown, sigmaup);
01930 }
01931
01932
01933 int BCModel::SetPrior(int index, TH1 * h, bool interpolate)
01934 {
01935
01936 if (index < 0 && index >= int(GetNParameters())) {
01937 BCLog::OutError("BCModel::SetPrior : Index out of range.");
01938 return 0;
01939 }
01940
01941
01942 if(h) {
01943
01944
01945 if (h->GetDimension() != 1) {
01946 BCLog::OutError(Form("BCModel::SetPrior : Histogram given for parameter %d is not 1D.",index));
01947 return 0;
01948 }
01949
01950
01951 h->Scale(1./h->Integral("width"));
01952
01953 if(fPriorContainer[index])
01954 delete fPriorContainer[index];
01955
01956
01957 fPriorContainer[index] = new TH1(*h);
01958
01959 if (interpolate)
01960 fPriorContainerInterpolate[index] = true;
01961
01962 fPriorContainerConstant[index] = false;
01963 }
01964
01965 RecalculatePriorConstant();
01966
01967
01968 ResetResults();
01969
01970
01971 return 1;
01972 }
01973
01974
01975 int BCModel::SetPrior(const char * name, TH1 * h, bool interpolate)
01976 {
01977
01978 int index = -1;
01979 for (unsigned int i = 0; i < GetNParameters(); i++)
01980 if (name == GetParameter(i)->GetName())
01981 index = i;
01982
01983
01984 return SetPrior(index, h, interpolate);
01985 }
01986
01987
01988 int BCModel::SetPriorConstant(int index)
01989 {
01990
01991 if (index < 0 && index >= int(GetNParameters())) {
01992 BCLog::OutError("BCModel::SetPriorConstant : Index out of range.");
01993 return 0;
01994 }
01995
01996 if(fPriorContainer[index]) {
01997 delete fPriorContainer[index];
01998 fPriorContainer[index] = 0;
01999 }
02000
02001
02002 fPriorContainerConstant[index] = true;
02003
02004 RecalculatePriorConstant();
02005
02006
02007 ResetResults();
02008
02009
02010 return 1;
02011 }
02012
02013
02014 int BCModel::SetPriorConstant(const char* name)
02015 {
02016
02017 int index = -1;
02018 for (unsigned int i = 0; i < GetNParameters(); i++) {
02019 if (name == GetParameter(i)->GetName()) {
02020 index = i;
02021 break;
02022 }
02023 }
02024
02025 if (index == -1) {
02026 BCLog::OutError(Form(
02027 "BCModel::SetPriorConstant : parameter '%s' doesn't exist.", name));
02028 return 0;
02029 }
02030
02031 return SetPriorConstant(index);
02032 }
02033
02034
02035 int BCModel::SetPriorConstantAll()
02036 {
02037
02038 int nPar = GetNParameters();
02039
02040 if (nPar == 0)
02041 BCLog::OutWarning("BCModel::SetPriorConstantAll : No parameters defined.");
02042
02043
02044 for (int i = 0; i < nPar; ++i) {
02045 if (fPriorContainer[i]) {
02046 delete fPriorContainer[i];
02047 fPriorContainer[i]=0;
02048 }
02049 fPriorContainerConstant[i] = true;
02050 }
02051
02052 RecalculatePriorConstant();
02053
02054
02055 ResetResults();
02056
02057
02058 return 1;
02059 }
02060
02061
02062 void BCModel::RecalculatePriorConstant()
02063 {
02064 fPriorConstantValue = 0.;
02065
02066
02067 int npar = GetNParameters();
02068
02069 int nconstant = 0;
02070
02071 for (int i=0; i<npar; ++i)
02072 if (fPriorContainerConstant[i]) {
02073 fPriorConstantValue -= log(GetParameter(i)->GetRangeWidth());
02074 ++nconstant;
02075 }
02076
02077 if (nconstant == npar)
02078 fPriorConstantAll = true;
02079 else
02080 fPriorConstantAll = false;
02081 }
02082
02083
02084 int BCModel::SetParameterRange(int index, double parmin, double parmax)
02085 {
02086
02087 if (index < 0 || index >= int(GetNParameters())) {
02088 BCLog::OutError("BCModel::SetParameterRange : Index out of range.");
02089 return 0;
02090 }
02091
02092
02093 GetParameter(index)->SetLowerLimit(parmin);
02094 GetParameter(index)->SetUpperLimit(parmax);
02095 fMCMCBoundaryMin[index] = parmin;
02096 fMCMCBoundaryMax[index] = parmax;
02097
02098
02099 ResetResults();
02100
02101
02102 return 1;
02103 }
02104
02105
02106 int BCModel::ResetResults()
02107 {
02108 BCIntegrate::IntegrateResetResults();
02109
02110 BCEngineMCMC::MCMCResetResults();
02111
02112
02113 return 1;
02114 }
02115
02116
02117 void BCModel::PrintSummary()
02118 {
02119 int nparameters = GetNParameters();
02120
02121
02122 std::cout << std::endl
02123 << " ---------------------------------" << std::endl
02124 << " Model : " << fName.data() << std::endl
02125 << " ---------------------------------" << std::endl
02126 << " Index : " << fIndex << std::endl
02127 << " Number of parameters : " << nparameters << std::endl << std::endl
02128 << " - Parameters : " << std::endl << std::endl;
02129
02130
02131 for (int i = 0; i < nparameters; i++)
02132 fParameterSet->at(i)->PrintSummary();
02133
02134
02135 if (GetBestFitParameters().size() > 0) {
02136 std::cout << std::endl << " - Best fit parameters :" << std::endl << std::endl;
02137
02138 for (int i = 0; i < nparameters; i++) {
02139 std::cout << " " << fParameterSet->at(i)->GetName().data()
02140 << " = " << GetBestFitParameter(i) << " (overall)" << std::endl;
02141 if ((int) fBestFitParametersMarginalized.size() == nparameters)
02142 std::cout << " "
02143 << fParameterSet->at(i)->GetName().data() << " = "
02144 << GetBestFitParameterMarginalized(i)
02145 << " (marginalized)" << std::endl;
02146 }
02147 }
02148
02149 std::cout << std::endl;
02150
02151
02152 if (fPValue >= 0) {
02153 double likelihood = Likelihood(GetBestFitParameters());
02154 std::cout << " - Model testing:" << std::endl << std::endl
02155 << " p(data|lambda*) = " << likelihood << std::endl
02156 << " p-value = " << fPValue << std::endl << std::endl;
02157 }
02158
02159
02160 if (fNormalization > 0)
02161 std::cout << " Normalization : " << fNormalization << std::endl;
02162 }
02163
02164
02165 void BCModel::PrintResults(const char * file)
02166 {
02167
02168
02169
02170 ofstream ofi(file);
02171
02172
02173 if (!ofi.is_open()) {
02174 std::cerr << "Couldn't open file " << file << std::endl;
02175 return;
02176 }
02177
02178
02179 int npar = MCMCGetNParameters();
02180 int nchains = MCMCGetNChains();
02181
02182
02183 bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0;
02184
02185 ofi << std::endl
02186 << " -----------------------------------------------------" << std::endl
02187 << " Summary" << std::endl
02188 << " -----------------------------------------------------" << std::endl
02189 << std::endl;
02190
02191 ofi << " Model summary" << std::endl << " =============" << std::endl
02192 << " Model: " << fName.data() << std::endl
02193 << " Number of parameters: " << npar << std::endl
02194 << " List of Parameters and ranges:" << std::endl;
02195 for (int i = 0; i < npar; ++i)
02196 ofi << " (" << i << ") Parameter \""
02197 << fParameterSet->at(i)->GetName().data() << "\"" << ": "
02198 << "(" << fParameterSet->at(i)->GetLowerLimit() << ", "
02199 << fParameterSet->at(i)->GetUpperLimit() << ")" << std::endl;
02200 ofi << std::endl;
02201
02202 ofi << " Results of the optimization" << std::endl
02203 << " ===========================" << std::endl
02204 << " Optimization algorithm used: "
02205 << DumpUsedOptimizationMethod()<< std::endl;
02206
02207 if (int(fBestFitParameters.size()) > 0) {
02208 ofi << " List of parameters and global mode:" << std::endl;
02209 for (int i = 0; i < npar; ++i) {
02210 ofi << " (" << i << ") Parameter \""
02211 << fParameterSet->at(i)->GetName().data() << "\": "
02212 << fBestFitParameters[i];
02213 if (int(fBestFitParameterErrors.size()) == npar)
02214 if(fBestFitParameterErrors[i]>=0.)
02215 ofi << " +- " << fBestFitParameterErrors[i];
02216 ofi << std::endl;
02217 }
02218 ofi << std::endl;
02219 }
02220 else {
02221 ofi << " No best fit information available." << std::endl;
02222 ofi << std::endl;
02223 }
02224
02225 if (fPValue >= 0.) {
02226 ofi << " Results of the model test" << std::endl
02227 << " =========================" << std::endl
02228 << " p-value at global mode: " << fPValue << std::endl << std::endl;
02229 }
02230
02231 if (fNormalization >= 0.) {
02232 ofi << " Results of the normalization" << std::endl
02233 << " ============================" << std::endl
02234 << " Integration method used:"
02235 << DumpIntegrationMethod() << std::endl;
02236 ofi << " Normalization factor: " << fNormalization << std::endl << std::endl;
02237 }
02238
02239
02240 if (!flag_conv && fMCMCFlagRun)
02241 ofi << " WARNING: the Markov Chain did not converge!" << std::endl
02242 << " Be cautious using the following results!" << std::endl
02243 << std::endl;
02244
02245
02246 if (fMCMCFlagRun && fMCMCFlagFillHistograms) {
02247 ofi << " Results of the marginalization" << std::endl
02248 << " ==============================" << std::endl
02249 << " List of parameters and properties of the marginalized"
02250 << std::endl << " distributions:" << std::endl;
02251 for (int i = 0; i < npar; ++i) {
02252 if (!fMCMCFlagsFillHistograms[i])
02253 continue;
02254
02255 BCH1D * bch1d = GetMarginalized(fParameterSet->at(i));
02256
02257 ofi << " (" << i << ") Parameter \""
02258 << fParameterSet->at(i)->GetName().data() << "\":" << std::endl
02259
02260 << " Mean +- sqrt(V): " << std::setprecision(4)
02261 << bch1d->GetMean() << " +- " << std::setprecision(4)
02262 << bch1d->GetRMS() << std::endl
02263
02264 << " Median +- central 68% interval: "
02265 << std::setprecision(4) << bch1d->GetMedian() << " + "
02266 << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian()
02267 << " - " << std::setprecision(4)
02268 << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl
02269
02270 << " (Marginalized) mode: " << bch1d->GetMode() << std::endl;
02271
02272 ofi << " 5% quantile: " << std::setprecision(4)
02273 << bch1d->GetQuantile(0.05) << std::endl
02274 << " 10% quantile: " << std::setprecision(4)
02275 << bch1d->GetQuantile(0.10) << std::endl
02276 << " 16% quantile: " << std::setprecision(4)
02277 << bch1d->GetQuantile(0.16) << std::endl
02278 << " 84% quantile: " << std::setprecision(4)
02279 << bch1d->GetQuantile(0.85) << std::endl
02280 << " 90% quantile: " << std::setprecision(4)
02281 << bch1d->GetQuantile(0.90) << std::endl
02282 << " 95% quantile: " << std::setprecision(4)
02283 << bch1d->GetQuantile(0.95) << std::endl;
02284
02285 std::vector<double> v;
02286 v = bch1d->GetSmallestIntervals(0.68);
02287 int ninter = int(v.size());
02288
02289 ofi << " Smallest interval(s) containing 68% and local modes:"
02290 << std::endl;
02291 for (int j = 0; j < ninter; j += 5)
02292 ofi << " (" << v[j] << ", " << v[j + 1]
02293 << ") (local mode at " << v[j + 3] << " with rel. height "
02294 << v[j + 2] << "; rel. area " << v[j + 4] << ")"
02295 << std::endl;
02296 ofi << std::endl;
02297 }
02298 }
02299 if (fMCMCFlagRun) {
02300 ofi << " Status of the MCMC" << std::endl << " ==================" << std::endl
02301 << " Convergence reached: " << (flag_conv ? "yes" : "no")
02302 << std::endl;
02303
02304 if (flag_conv)
02305 ofi << " Number of iterations until convergence: "
02306 << MCMCGetNIterationsConvergenceGlobal() << std::endl;
02307 else
02308 ofi << " WARNING: the Markov Chain did not converge! Be\n"
02309 << " cautious using the following results!" << std::endl
02310 << std::endl;
02311 ofi << " Number of chains: " << MCMCGetNChains() << std::endl
02312 << " Number of iterations per chain: " << MCMCGetNIterationsRun() << std::endl
02313 << " Average efficiencies:" << std::endl;
02314
02315 std::vector<double> efficiencies;
02316 efficiencies.assign(npar, 0.);
02317
02318 for (int ipar = 0; ipar < npar; ++ipar)
02319 for (int ichain = 0; ichain < nchains; ++ichain) {
02320 int index = ichain * npar + ipar;
02321 efficiencies[ipar] +=
02322 double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index)
02323 + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.;
02324 }
02325
02326 for (int ipar = 0; ipar < npar; ++ipar)
02327 ofi << " (" << ipar << ") Parameter \""
02328 << fParameterSet->at(ipar)->GetName().data() << "\": "
02329 << efficiencies.at(ipar) << "%" << std::endl;
02330 ofi << std::endl;
02331 }
02332
02333 ofi << " -----------------------------------------------------" << std::endl
02334 << " Notation:" << std::endl
02335 << " Mean : mean value of the marg. pdf" << std::endl
02336 << " Median : maximum of the marg. pdf" << std::endl
02337 << " Marg. mode : most probable value of the marg. pdf" << std::endl
02338 << " V : Variance of the marg. pdf" << std::endl
02339 << " Quantiles : most commonly used quantiles" <<std::endl
02340 << " -----------------------------------------------------" << std::endl
02341 << std::endl;
02342
02343
02344
02345 }
02346
02347
02348 void BCModel::PrintShortFitSummary(int chi2flag)
02349 {
02350 BCLog::OutSummary("---------------------------------------------------");
02351 BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data()));
02352 BCLog::OutSummary(Form(" Number of parameters: Npar = %i", GetNParameters()));
02353 BCLog::OutSummary(Form(" Number of data points: Ndata = %i", GetNDataPoints()));
02354 BCLog::OutSummary(" Number of degrees of freedom:");
02355 BCLog::OutSummary(Form(" NDoF = Ndata - Npar = %i", GetNDataPoints() - GetNParameters()));
02356
02357 BCLog::OutSummary(" Best fit parameters (global):");
02358 for (unsigned int i = 0; i < GetNParameters(); ++i)
02359 BCLog::OutSummary(Form(" %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i)));
02360
02361 BCLog::OutSummary(" Goodness-of-fit test:");
02362 BCLog::OutSummary(Form(" p-value = %.3g", GetPValue()));
02363 if (chi2flag) {
02364 BCLog::OutSummary(Form(" p-value corrected for NDoF = %.3g", GetPValueNDoF()));
02365 BCLog::OutSummary(Form(" chi2 / NDoF = %.3g", GetChi2NDoF()));
02366 }
02367 BCLog::OutSummary("---------------------------------------------------");
02368 }
02369
02370
02371 void BCModel::PrintHessianMatrix(std::vector<double> parameters)
02372 {
02373
02374 if (parameters.size() != GetNParameters()) {
02375 BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
02376 return;
02377 }
02378
02379
02380 std::cout << std::endl << " Hessian matrix elements: " << std::endl << " Point: ";
02381
02382 for (int i = 0; i < int(parameters.size()); i++)
02383 std::cout << parameters.at(i) << " ";
02384 std::cout << std::endl;
02385
02386
02387 for (unsigned int i = 0; i < GetNParameters(); i++)
02388 for (unsigned int j = 0; j < i; j++) {
02389
02390 double hessianmatrixelement = HessianMatrixElement(
02391 fParameterSet->at(i), fParameterSet->at(j), parameters);
02392
02393
02394 std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl;
02395 }
02396 }
02397
02398
02399 BCDataPoint * BCModel::VectorToDataPoint(std::vector<double> data)
02400 {
02401 int sizeofvector = int(data.size());
02402 BCDataPoint * datapoint = new BCDataPoint(sizeofvector);
02403 datapoint->SetValues(data);
02404 return datapoint;
02405 }
02406
02407
02408 int BCModel::CompareStrings(const char * string1, const char * string2)
02409 {
02410 int flag_same = 0;
02411
02412 if (strlen(string1) != strlen(string2))
02413 return -1;
02414
02415 for (int i = 0; i < int(strlen(string1)); i++)
02416 if (string1[i] != string2[i])
02417 flag_same = -1;
02418
02419 return flag_same;
02420 }
02421
02422 void BCModel::StoreMode()
02423 {
02424
02425 double probmax = fMCMCprobMax.at(0);
02426 int probmaxindex = 0;
02427
02428
02429 for (int i = 1; i < fMCMCNChains; ++i) {
02430 if (fMCMCprobMax.at(i) > probmax) {
02431 probmax = fMCMCprobMax.at(i);
02432 probmaxindex = i;
02433 }
02434 }
02435
02436
02437 if (fBestFitParameters.empty() || probmax > LogEval(fBestFitParameters)) {
02438 fBestFitParameters.assign(fMCMCNParameters, 0.0);
02439 for (int i = 0; i < fMCMCNParameters; ++i)
02440 fBestFitParameters[i] = fMCMCxMax[probmaxindex * fMCMCNParameters + i];
02441
02442 SetOptimizationMethodMode(BCIntegrate::kOptMetropolis);
02443 }
02444 }
02445
02446
02447