• Main Page
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

BCModel.cxx

Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 2008-2011, Daniel Kollar and Kevin Kroeninger.
00003  * All rights reserved.
00004  *
00005  * For the licensing terms see doc/COPYING.
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    // return this
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    // loop over x and y bins
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    // define new graph
00378    int nx = fErrorBandXY->GetNbinsX();
00379 
00380    TGraph * graph = new TGraph(2 * nx);
00381    graph->SetFillStyle(1001);
00382    graph->SetFillColor(kYellow);
00383 
00384    // get error bands
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    // copy existing histogram
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    // loop over x bins
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    // define new graph
00440    int nx = fErrorBandXY->GetNbinsX();
00441    TGraph * graph = new TGraph(nx);
00442 
00443    // loop over x values
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    // define new graph
00462    TGraph * graph = new TGraph(n + 1);
00463 
00464    double dx = (xmax - xmin) / (double) n;
00465 
00466    // loop over x values
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    // create new data set consisting of a single data point
00503    BCDataSet * dataset = new BCDataSet();
00504 
00505    // add the data point
00506    dataset->AddDataPoint(datapoint);
00507 
00508    // set this new data set
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    // check if data set exists
00525    if (!fDataSet) {
00526       BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
00527       return;
00528    }
00529 
00530    // check if index is within range
00531    if (index < 0 || index > fDataSet->GetDataPoint(0)->GetNValues()) {
00532       BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
00533       return;
00534    }
00535 
00536    // check if boundary data points exist
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    // set boundaries
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    // clear x-values
00561    fErrorBandX.clear();
00562 
00563    // copy data x-values
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    // create new parameter
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    // check if parameter set exists
00585    if (!fParameterSet) {
00586       BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
00587       return ERROR_PARAMETERSETDOESNOTEXIST;
00588    }
00589 
00590    // check if parameter with same name exists
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    // define index of new parameter
00604    unsigned int index = fParameterSet->size();
00605    parameter->SetIndex(index);
00606 
00607    // add parameter to parameter container
00608    fParameterSet->push_back(parameter);
00609 
00610    // add empty object to prior container
00611    fPriorContainer.push_back(0); 
00612 
00613    // don't interpolate the prior histogram by default
00614    fPriorContainerInterpolate.push_back(false);
00615 
00616    // prior assumed to be non-constant in general case
00617    fPriorContainerConstant.push_back(false);
00618 
00619    // add parameters to integration methods
00620    SetParameters(fParameterSet);
00621 
00622    // reset results
00623    ResetResults();
00624 
00625    return 0;
00626 }
00627 
00628 // ---------------------------------------------------------
00629 double BCModel::LogProbabilityNN(std::vector<double> parameters)
00630 {
00631    // add log of conditional probability
00632    double logprob = LogLikelihood(parameters);
00633 
00634    // add log of prior probability
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    // check if normalized
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    // add log of conditional probabilities event-by-event
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    // get number of parameters
00675    int npar = GetNParameters();
00676 
00677    // loop over all 1-d priors
00678    for (int i = 0; i < npar; ++i) {
00679       if (fPriorContainer[i]) {
00680          // check what type of object is stored
00681          TF1 * f = dynamic_cast<TF1*>(fPriorContainer[i]);
00682          TH1 * h = dynamic_cast<TH1*>(fPriorContainer[i]);
00683 
00684          if (f) // TF1
00685             logprob += log(f->Eval(parameters[i]));
00686          else if (h) { // TH1
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())); // this should never happen
00697       }
00698       // use constant only if user has defined it
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    // add the contribution from constant priors in one step
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    // initialize BCIntegrate if not done already
00748    if (n == 0) {
00749       SetParameters(fParameterSet);
00750       n = GetNvar();
00751    }
00752 
00753    // integrate and get best fit parameters
00754    // maybe we have to remove the mode finding from here in the future
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    // check if vectors are of equal size
00766    if (!parameters.size() == fParameterSet->size())
00767       return ERROR_INVALIDNUMBEROFPARAMETERS;
00768 
00769    // check if parameters are within limits
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    // this implementation is CLEARLY not good we have to work on this.
00794 
00795    BCLog::OutSummary(Form("Model \'%s\': Finding mode using %s",GetName().data(), DumpOptimizationMethod().c_str()));
00796 
00797    // synchronize parameters in BCIntegrate
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    // synchronize parameters in BCIntegrate
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    // prepare function fitting
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    // run the Markov chains
00960    MCMCMetropolis();
00961 
00962    // store the mode found by the chains
00963    StoreMode();
00964 
00965    return 1;
00966 }
00967 
00968 // ---------------------------------------------------------
00969 BCH1D * BCModel::GetMarginalized(BCParameter * parameter)
00970 {
00971    if (!parameter) {
00972 // don't print any error message, should be done upstream
00973 //      BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist.");
00974       return 0;
00975    }
00976 
00977    int index = parameter->GetIndex();
00978    if (fMCMCFlagsFillHistograms[index] == false) {
00979 // don't print any error message, should be done upstream
00980 //      BCLog::OutError(Form(
00981 //            "BCModel::GetMarginalized : Distribuion for '%s' not filled.",
00982 //            parameter->GetName().data()));
00983       return 0;
00984    }
00985 
00986    // get histogram
00987    TH1D * hist = MCMCGetH1Marginalized(index);
00988    if (!hist)
00989       return 0;
00990 
00991    // set axis labels
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    // set histogram
00998    BCH1D * hprob = new BCH1D();
00999    hprob->SetHistogram(hist);
01000 
01001    // set best fit parameter
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    // We reset the MCMCEngine here for the moment.
01024    // In the future maybe we only want to do this if the engine
01025    // wans't initialized at all or when there were some changes
01026    // in the model.
01027    // But maybe we want reset everything since we're overwriting
01028    // the marginalized distributions anyway.
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    // if there's only one parameter, we just want to call Print()
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; // default canvas width
01178    int c_height = 600; // default canvas height
01179 
01180    int type = 112; // landscape
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    // calculate number of plots
01205    int nplots2d = npar * (npar - 1) / 2;
01206    int nplots = npar + nplots2d;
01207 
01208    // give out warning if too many plots
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    // setup the canvas and postscript file
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    // draw all 1D distributions
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       // get corresponding parameter
01235       BCParameter * a = GetParameter(i);
01236 
01237       // check if histogram exists
01238       if (!GetMarginalized(a))
01239          continue;
01240 
01241       // check if histogram is filled
01242       if (GetMarginalized(a)->GetHistogram()->Integral() <= 0)
01243          continue;
01244 
01245       // if current page is full, switch to new page
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       // go to next pad
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       // draw all the 2D distributions
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          // get corresponding parameters
01284          BCParameter * a = GetParameter(i);
01285          BCParameter * b = GetParameter(j);
01286 
01287          // check if histogram exists
01288          if (!GetMarginalized(a, b))
01289             continue;
01290 
01291          // check if histogram is filled
01292          if (GetMarginalized(a, b)->GetHistogram()->Integral() <= 0)
01293             continue;
01294 
01295          // if current page is full, switch to new page
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          // go to next pad
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    // return total number of drawn histograms
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 // don't print any error message, should be done upstream
01346 //      BCLog::OutError(
01347 //            Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.",
01348 //                  par1->GetName().data(), par2->GetName().data()));
01349       return 0;
01350    }
01351 
01352    if (index1 == index2) {
01353 // don't print any error message, should be done upstream
01354 //      BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
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    // get histogram
01371    TH2D * hist = MCMCGetH2Marginalized(index1, index2);
01372 
01373    if (hist == 0)
01374       return 0;
01375 
01376    BCH2D * hprob = new BCH2D();
01377 
01378    // set axis labels
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    // set histogram
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    // look up corresponding p value
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    // number of observations
01432    int n = GetNDataPoints();
01433 
01434    if (nBins < 0)
01435       nBins = NumberBins();
01436 
01437    // fixed width quantiles, including final point!
01438    std::vector<double> a;
01439    for (int i = 0; i <= nBins; i++)
01440       a.push_back(i / double(nBins));
01441 
01442    // determine the bin counts and fill the histogram with data using the CDF
01443    TH1I * hist = new TH1I("h1", "h1 title", nBins, 0., 1.);
01444 
01445    // discrete case requires randomization to allocate counts of bins that cover more
01446    // than one quantile
01447    if (flag_discrete) {
01448       // loop over observations, each may have different likelihood and CDF
01449       for (int j = 0; j < n; j++) {
01450          // actual value
01451             double CDFval = CDF(par, j, false);
01452          // for the bin just before
01453          double CDFlower = CDF(par, j, true);
01454 
01455          // what quantiles q are covered, count from q_1 to q_{nBins}
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          // simplest case: observation bin entirely contained in one quantile
01472          if (qMin == qMax) {
01473             // watch out for overflow because CDF exactly 1
01474             if (CDFval < 1)
01475                hist->Fill(CDFval);
01476             else
01477                hist->AddBinContent(qMax);
01478 
01479             continue; // this observation finished
01480          }
01481 
01482          // if more than quantile is covered need more work:
01483          // determine probabilities of this observation to go for each quantile covered
01484          // as follows: If each quantile has size 0.25 and the CDF(integral of likelihood)
01485          // for current observation gives gives 0.27, but for observation-1 we would have
01486          // 0.20, then 5/7 of the 7% go for first quantile and 2/7 for the second.
01487          // This extend to bins covering more than two quantiles
01488          std::vector<double> prob;
01489          // normalization
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             // default case
01502             prob.push_back(norm * (a[i] - a[i - 1]));
01503          }
01504          // have distribution, use inverse-transform method
01505          double U = fRandom->Rndm();
01506          // build up the integral (CDF)
01507          for (uint i = 1; i < prob.size(); i++)
01508             prob[i] += prob[i - 1];
01509          // and search with linear comput. complexity
01510          for (uint i = 0; i < prob.size(); i++) {
01511             // we finally allocate the count, as center of quantile
01512             if (U <= prob[i]) {
01513                hist->Fill((a[qMin + i - 1] + a[qMin + i]) / 2.);
01514                break;
01515             }
01516          }
01517       }
01518    }
01519    else { //continuous case is simple
01520       for (int j = 0; j < n; j++)
01521             hist->Fill(CDF(par, j, false));
01522    }
01523 
01524    // calculate chi^2
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    // model parameters
01543    int nPar = (int)GetNParameters();
01544    std::vector<double> param(nPar);
01545 
01546    //parameters saved in branches should be the same
01547    int nParBranches = -1;
01548    tree->SetBranchAddress("fNParameters", &nParBranches);
01549 
01550    //assume all events have same number of parameters, so check only first
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    // buffer to construct correct branchnames for parameters, e.g. "fParameter3"
01561    char * branchName = new char[10 + nPar];
01562    // set up variables filled for each sample of parameters
01563    // assume same order as in model
01564    for (int i = 0; i < (int) nPar; i++) {
01565       +sprintf(branchName, "fParameter%d", i);
01566       tree->SetBranchAddress(branchName, &param[i]);
01567    }
01568 
01569    // get the p value from Johson's definition for each param from posterior
01570    long nEntries = tree->GetEntries();
01571 
01572    // RN ~ chi2 with K-1 DoF needed for comparison
01573    std::vector<double> randoms(nEntries);
01574    int K = NumberBins();
01575    BCMath::RandomChi2(randoms, K - 1);
01576 
01577    // number of Johnson chi2 values bigger than reference
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       // compare to set of chi2 variables
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    //calculate the ECDF from the 1D data
01621    std::vector<double> yData = fDataSet->GetDataComponents(index);
01622    TH1D * ECDF = BCMath::ECDF(yData);
01623 
01624    int N = GetNDataPoints();
01625 
01626    // calculated expected CDF for unique points only
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    // find maximum distance
01640    double distMax = 0.;
01641 
01642    // current distance
01643    double dist = 0.;
01644 
01645    std::set<double>::const_iterator iter = uniqueObservations.begin();
01646    for (int iBin = 0; iBin < nUnique; ++iBin) {
01647       // distance between current points
01648       dist = TMath::Abs(*iter - ECDF->GetBinContent(iBin + 1));
01649       // update maximum if necessary
01650       distMax = TMath::Max(dist, distMax);
01651 
01652 //      BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
01653 //         "expected vs empirical (%f vs %f)", *iter, ECDF->GetBinContent(iBin
01654 //            + 1)));
01655 
01656       // advance to next entry in the set
01657       ++iter;
01658    }
01659 
01660    // correct for total #points, not unique #points.
01661    // would need sqrt(n1*n2/(n1+n2)) if we had two experimental datasets
01662    double z = distMax * sqrt(N);
01663 
01664    fPValue = TMath::KolmogorovProb(z);
01665 
01666 //   BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
01667 //      "max distance vs corrected (%f vs %f)", distMax, z));
01668 
01669    // clean up
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    // print log
01681    BCLog::OutSummary("Do goodness-of-fit-test");
01682 
01683    // create model test
01684    BCGoFTest * goftest = new BCGoFTest("modeltest");
01685 
01686    // set this model as the model to be tested
01687    goftest->SetTestModel(this);
01688 
01689    // set the point in parameter space which is tested an initialize
01690    // the model testing
01691    if (!goftest->SetTestPoint(par))
01692       return 0;
01693 
01694    // disable the creation of histograms to save _a lot_ of memory
01695    // (histograms are not needed for p-value calculation)
01696    goftest->MCMCSetFlagFillHistograms(false);
01697 
01698    // set parameters of the MCMC for the GoFTest
01699    goftest->MCMCSetNChains(fGoFNChains);
01700    goftest->MCMCSetNIterationsMax(fGoFNIterationsMax);
01701    goftest->MCMCSetNIterationsRun(fGoFNIterationsRun);
01702 
01703    // get p-value
01704    fPValue = goftest->GetCalculatedPValue(flag_histogram);
01705 
01706    // get histogram
01707    if (flag_histogram) {
01708       hist = new BCH1D();
01709       hist->SetHistogram(goftest->GetHistogramLogProb());
01710    }
01711 
01712    // delete model test
01713    delete goftest;
01714 
01715    // return histogram
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    // check number of entries in vector
01729    if (point.size() != GetNParameters()) {
01730       BCLog::OutError("BCModel::HessianMatrixElement : Invalid number of entries in the vector.");
01731       return -1;
01732    }
01733 
01734    // define steps
01735    double nsteps = 1e5;
01736    double dx1 = par1->GetRangeWidth() / nsteps;
01737    double dx2 = par2->GetRangeWidth() / nsteps;
01738 
01739    // define points at which to evaluate
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    // calculate probability at these points
01761    double ppp = Likelihood(xpp);
01762    double ppm = Likelihood(xpm);
01763    double pmp = Likelihood(xmp);
01764    double pmm = Likelihood(xmm);
01765 
01766    // return derivative
01767    return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
01768 }
01769 
01770 // ---------------------------------------------------------
01771 void BCModel::FixDataAxis(unsigned int index, bool fixed)
01772 {
01773    // check if index is within range
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    // check if index is within range
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    // check index range
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    // copy function
01811    fPriorContainer[index] = new TF1(*f);
01812 
01813    fPriorContainerConstant[index] = false;
01814 
01815    RecalculatePriorConstant();
01816 
01817    // reset all results
01818    ResetResults();
01819 
01820    // no error 
01821    return 1; 
01822 }
01823 
01824 // ---------------------------------------------------------
01825 int BCModel::SetPrior(const char * name, TF1 * f)
01826 {
01827    // find index
01828    int index = -1;
01829    for (unsigned int i = 0; i < GetNParameters(); i++)
01830       if (name == GetParameter(i)->GetName())
01831         index = i;
01832 
01833    // set prior
01834    return SetPrior(index, f);
01835 }
01836 
01837 // ---------------------------------------------------------
01838 int BCModel::SetPriorDelta(int index, double value)
01839 {
01840    // set range to value
01841    SetParameterRange(index, value, value);
01842 
01843    // set prior
01844    return SetPriorConstant(index);
01845 }
01846 
01847 // ---------------------------------------------------------
01848 int BCModel::SetPriorDelta(const char* name, double value)
01849 {
01850   // find index
01851    int index = -1;
01852    for (unsigned int i = 0; i < GetNParameters(); i++)
01853       if (name == GetParameter(i)->GetName())
01854         index = i;
01855 
01856    // set prior
01857     return SetPriorDelta(index, value);
01858 }
01859 
01860 // ---------------------------------------------------------
01861 int BCModel::SetPriorGauss(int index, double mean, double sigma)
01862 {
01863    // check index range
01864    if (index < 0 || index >= int(GetNParameters())) {
01865       BCLog::OutError("BCModel::SetPriorGauss : Index out of range."); 
01866       return 0;
01867    }
01868 
01869    // create new function
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    // set prior
01878    return SetPrior(index, f); 
01879 }
01880 
01881 // ---------------------------------------------------------
01882 int BCModel::SetPriorGauss(const char* name, double mean, double sigma)
01883 {
01884    // find index
01885    int index = -1;
01886    for (unsigned int i = 0; i < GetNParameters(); i++)
01887       if (name == GetParameter(i)->GetName())
01888          index = i;
01889 
01890    // set prior
01891    return SetPriorGauss(index, mean, sigma); 
01892 }
01893 
01894 // ---------------------------------------------------------
01895 int BCModel::SetPriorGauss(int index, double mean, double sigmadown, double sigmaup)
01896 {
01897    // check index range
01898    if (index < 0 || index >= int(GetNParameters())) {
01899       BCLog::OutError("BCModel::SetPriorGauss : Index out of range."); 
01900       return 0;
01901    }
01902 
01903    // create new function
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    // set prior
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    // find index
01923    int index = -1;
01924    for (unsigned int i = 0; i < GetNParameters(); i++)
01925       if (name == GetParameter(i)->GetName())
01926         index = i;
01927 
01928    // set prior
01929    return SetPriorGauss(index, mean, sigmadown, sigmaup); 
01930 }
01931 
01932 // ---------------------------------------------------------
01933 int BCModel::SetPrior(int index, TH1 * h, bool interpolate)
01934 {
01935    // check index range
01936    if (index < 0 && index >= int(GetNParameters())) {
01937       BCLog::OutError("BCModel::SetPrior : Index out of range."); 
01938       return 0;
01939    }
01940 
01941    // if the histogram exists
01942    if(h) {
01943 
01944       // check if histogram is 1d
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       // normalize the histogram
01951       h->Scale(1./h->Integral("width"));
01952 
01953       if(fPriorContainer[index])
01954          delete fPriorContainer[index];
01955 
01956       // set function
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    // reset all results
01968    ResetResults();
01969 
01970    // no error 
01971    return 1; 
01972 }
01973 
01974 // ---------------------------------------------------------
01975 int BCModel::SetPrior(const char * name, TH1 * h, bool interpolate)
01976 {
01977    // find index
01978    int index = -1;
01979    for (unsigned int i = 0; i < GetNParameters(); i++)
01980       if (name == GetParameter(i)->GetName())
01981         index = i;
01982 
01983    // set prior
01984    return SetPrior(index, h, interpolate);
01985 }
01986 
01987 // ---------------------------------------------------------
01988 int BCModel::SetPriorConstant(int index)
01989 {
01990    // check index range
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    // set prior to a constant
02002    fPriorContainerConstant[index] = true;
02003 
02004    RecalculatePriorConstant();
02005 
02006    // reset all results
02007    ResetResults();
02008 
02009    // no error
02010    return 1;
02011 }
02012 
02013 // ---------------------------------------------------------
02014 int BCModel::SetPriorConstant(const char* name)
02015 {
02016    // find index
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    // get number of parameters
02038    int nPar = GetNParameters();
02039 
02040    if (nPar == 0)
02041       BCLog::OutWarning("BCModel::SetPriorConstantAll : No parameters defined.");
02042 
02043    // loop over all 1-d priors
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    // reset all results
02055    ResetResults();
02056 
02057    // no error
02058    return 1;
02059 }
02060 
02061 // ---------------------------------------------------------
02062 void BCModel::RecalculatePriorConstant()
02063 {
02064    fPriorConstantValue = 0.;
02065 
02066    // get number of parameters
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    // check index
02087    if (index < 0 || index >= int(GetNParameters())) {
02088       BCLog::OutError("BCModel::SetParameterRange : Index out of range.");
02089       return 0; 
02090    }
02091 
02092    // set parameter ranges in BAT 
02093    GetParameter(index)->SetLowerLimit(parmin); 
02094    GetParameter(index)->SetUpperLimit(parmax); 
02095    fMCMCBoundaryMin[index] = parmin; 
02096    fMCMCBoundaryMax[index] = parmax; 
02097 
02098    // reset results
02099    ResetResults();
02100 
02101    // no error 
02102    return 1; 
02103 }
02104 
02105 // ---------------------------------------------------------
02106 int BCModel::ResetResults()
02107 {
02108    BCIntegrate::IntegrateResetResults();
02109   
02110    BCEngineMCMC::MCMCResetResults();
02111 
02112    // no error
02113    return 1;
02114 }
02115 
02116 // ---------------------------------------------------------
02117 void BCModel::PrintSummary()
02118 {
02119    int nparameters = GetNParameters();
02120 
02121    // model summary
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    // parameter summary
02131    for (int i = 0; i < nparameters; i++)
02132       fParameterSet->at(i)->PrintSummary();
02133 
02134    // best fit parameters
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    // model testing
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    // normalization
02160    if (fNormalization > 0)
02161       std::cout << "     Normalization        : " << fNormalization << std::endl;
02162 }
02163 
02164 // ---------------------------------------------------------
02165 void BCModel::PrintResults(const char * file)
02166 {
02167    // print summary of Markov Chain Monte Carlo
02168 
02169    // open file
02170    ofstream ofi(file);
02171 
02172    // check if file is open
02173    if (!ofi.is_open()) {
02174       std::cerr << "Couldn't open file " << file << std::endl;
02175       return;
02176    }
02177 
02178    // number of parameters and chains
02179    int npar = MCMCGetNParameters();
02180    int nchains = MCMCGetNChains();
02181 
02182    // check convergence
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    // give warning if MCMC did not converge
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    // print results of marginalization (if MCMC was run)
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    // close file
02344    //   ofi.close;
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    // check number of entries in vector
02374    if (parameters.size() != GetNParameters()) {
02375       BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
02376       return;
02377    }
02378 
02379    // print to screen
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    // loop over all parameter pairs
02387    for (unsigned int i = 0; i < GetNParameters(); i++)
02388       for (unsigned int j = 0; j < i; j++) {
02389          // calculate Hessian matrix element
02390          double hessianmatrixelement = HessianMatrixElement(
02391                fParameterSet->at(i), fParameterSet->at(j), parameters);
02392 
02393          // print to screen
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    //start with max. of posterior at first chain
02425    double probmax = fMCMCprobMax.at(0);
02426    int probmaxindex = 0;
02427 
02428    // loop over all chains and find the maximum point
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    // save if improved the log posterior
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 

Generated by  doxygen 1.7.1