Protected Attributes | Private Member Functions | Private Attributes

BCModel Class Reference

The base class for all user-defined models. More...

#include <BCModel.h>

Inherits BCIntegrate.

Inherited by BCEfficiencyFitter, BCGoFTest, BCGraphFitter, BCHistogramFitter, BCMTF, BCRooInterface, BCSummaryPriorModel, and BCTemplateFitter.

Collaboration diagram for BCModel:
Collaboration graph
[legend]

List of all members.

Public Member Functions

Constructors and destructors

 BCModel ()
 BCModel (const char *name)
 BCModel (const BCModel &bcmodel)
virtual ~BCModel ()
Assignment operators

BCModel & operator= (const BCModel &bcmodel)
Member functions (get)

std::string GetName ()
int GetIndex ()
double GetModelAPrioriProbability ()
double GetModelAPosterioriProbability ()
double GetNormalization ()
BCDataSet * GetDataSet ()
BCDataPoint * GetDataPointLowerBoundaries ()
BCDataPoint * GetDataPointUpperBoundaries ()
double GetDataPointLowerBoundary (unsigned int index)
double GetDataPointUpperBoundary (unsigned int index)
bool GetFlagBoundaries ()
int GetNDataPoints ()
BCDataPoint * GetDataPoint (unsigned int index)
unsigned int GetNDataPointsMinimum ()
unsigned int GetNDataPointsMaximum ()
unsigned int GetNParameters ()
BCParameter * GetParameter (int index)
BCParameter * GetParameter (const char *name)
BCParameterSet * GetParameterSet ()
double GetBestFitParameter (unsigned int index)
double GetBestFitParameterError (unsigned int index)
std::vector< double > GetBestFitParameters ()
std::vector< double > GetBestFitParameterErrors ()
double GetBestFitParameterMarginalized (unsigned int index)
std::vector< double > GetBestFitParametersMarginalized ()
TH2D * GetErrorBandXY ()
TH2D * GetErrorBandXY_yellow (double level=.68, int nsmooth=0)
std::vector< double > GetErrorBand (double level)
TGraph * GetErrorBandGraph (double level1, double level2)
TGraph * GetFitFunctionGraph (const std::vector< double > &parameters)
TGraph * GetFitFunctionGraph ()
TGraph * GetFitFunctionGraph (const std::vector< double > &parameters, double xmin, double xmax, int n=1000)
bool GetFixedDataAxis (unsigned int index)
Member functions (set)

void SetName (const char *name)
void SetIndex (int index)
void SetParameterSet (BCParameterSet *parset)
int SetParameterRange (int index, double parmin, double parmax)
void SetModelAPrioriProbability (double probability)
void SetModelAPosterioriProbability (double probability)
void SetNormalization (double norm)
void SetDataSet (BCDataSet *dataset)
void SetSingleDataPoint (BCDataPoint *datapoint)
void SetSingleDataPoint (BCDataSet *dataset, unsigned int index)
void SetNDataPointsMinimum (unsigned int minimum)
void SetNDataPointsMaximum (unsigned int maximum)
void SetDataBoundaries (unsigned int index, double lowerboundary, double upperboundary, bool fixed=false)
void SetErrorBandContinuous (bool flag)
void SetNbins (const char *parname, int nbins)
int SetPrior (int index, TF1 *f)
int SetPrior (const char *name, TF1 *f)
int SetPriorDelta (int index, double value)
int SetPriorDelta (const char *name, double value)
int SetPriorGauss (int index, double mean, double sigma)
int SetPriorGauss (const char *name, double mean, double sigma)
int SetPriorGauss (int index, double mean, double sigmadown, double sigmaup)
int SetPriorGauss (const char *name, double mean, double sigmadown, double sigmaup)
int SetPrior (int index, TH1 *h, bool flag=false)
int SetPrior (const char *name, TH1 *h, bool flag=false)
int SetPriorConstant (int index)
int SetPriorConstant (const char *name)
int SetPriorConstantAll ()
Member functions (miscellaneous methods)

int AddParameter (const char *name, double lowerlimit, double upperlimit)
int AddParameter (BCParameter *parameter)
double APrioriProbability (const std::vector< double > &parameters)
virtual double LogAPrioriProbability (const std::vector< double > &parameters)
double Likelihood (const std::vector< double > &params)
virtual double LogLikelihood (const std::vector< double > &params)=0
double ProbabilityNN (const std::vector< double > &params)
double LogProbabilityNN (const std::vector< double > &parameter)
double Probability (const std::vector< double > &parameter)
double LogProbability (const std::vector< double > &parameter)
virtual double SamplingFunction (const std::vector< double > &parameters)
double Eval (const std::vector< double > &parameters)
double LogEval (const std::vector< double > &parameters)
double EvalSampling (const std::vector< double > &parameters)
double Normalize ()
int CheckParameters (const std::vector< double > &parameters)
void FindMode (std::vector< double > start=std::vector< double >(0))
void FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1)
void WriteMode (const char *file)
int ReadMode (const char *file)
int ReadMarginalizedFromFile (const char *file)
int ReadErrorBandFromFile (const char *file)
int MarginalizeAll ()
BCH1D * GetMarginalized (BCParameter *parameter)
BCH1D * GetMarginalized (const char *name)
BCH2D * GetMarginalized (BCParameter *parameter1, BCParameter *parameter2)
BCH2D * GetMarginalized (const char *name1, const char *name2)
int PrintAllMarginalized1D (const char *filebase)
int PrintAllMarginalized2D (const char *filebase)
int PrintAllMarginalized (const char *file, unsigned int hdiv=1, unsigned int ndiv=1)
virtual void CorrelateDataPointValues (std::vector< double > &x)
double GetPvalueFromChi2 (const std::vector< double > &par, int sigma_index)
double GetPvalueFromChi2Johnson (std::vector< double > par)
double GetPvalueFromKolmogorov (const std::vector< double > &par, int index)
double GetChi2Johnson (std::vector< double > par, const int nBins=-1)
double GetAvalueFromChi2Johnson (TTree *tree, TH1D *distribution=0)
double GetPvalueFromChi2NDoF (std::vector< double > par, int sigma_index)
BCH1D * CalculatePValue (std::vector< double > par, bool flag_histogram=false)
double GetPValue ()
double GetPValueNDoF ()
double GetChi2NDoF ()
std::vector< double > GetChi2Runs (int dataIndex, int sigmaIndex)
void SetGoFNIterationsMax (int n)
void SetGoFNIterationsRun (int n)
void SetGoFNChains (int n)
double HessianMatrixElement (BCParameter *parameter1, BCParameter *parameter2, std::vector< double > point)
void PrintSummary ()
void PrintResults (const char *file)
void PrintShortFitSummary (int chi2flag=0)
void PrintHessianMatrix (std::vector< double > parameters)
void FixDataAxis (unsigned int index, bool fixed)
virtual double CDF (const std::vector< double > &, int, bool)
int ResetResults ()

Protected Attributes

int fIndex
std::string fName
double fModelAPriori
double fModelAPosteriori
BCParameterSet * fParameterSet
BCDataSet * fDataSet
unsigned int fNDataPointsMinimum
unsigned int fNDataPointsMaximum
double fPValue
double fChi2NDoF
double fPValueNDoF
bool flag_discrete
int fGoFNIterationsMax
int fGoFNIterationsRun
int fGoFNChains
std::vector< TNamed * > fPriorContainer
bool fPriorConstantAll
double fPriorConstantValue
std::vector< bool > fPriorContainerConstant
std::vector< bool > fPriorContainerInterpolate

Private Member Functions

int CompareStrings (const char *string1, const char *string2)
int NumberBins ()
void RecalculatePriorConstant ()
void StoreMode ()
BCDataPoint * VectorToDataPoint (const std::vector< double > &data)

Private Attributes

double fNormalization

Detailed Description

The base class for all user-defined models.

Author:
Daniel Kollar
Kevin Kröninger
Version:
1.0
Date:
08.2008 This class represents a model. It contains a container of parameters, their prior distributions and the likelihood for the parameters. The methods which implement the prior and the likelihood have to be overloaded by the user in the user defined model class derived from this class.

Definition at line 52 of file BCModel.h.

Constructor & Destructor Documentation

BCModel::BCModel (  ) 

The default constructor.

Definition at line 77 of file BCModel.cxx.

   : BCIntegrate()
{
   fNormalization = -1.;
   fDataSet = 0;
   fParameterSet = new BCParameterSet();

   fIndex = -1;
   fPValue = -1;
   fPValueNDoF = -1;
   fChi2NDoF = -1;

   fName = "model";
   fDataPointUpperBoundaries = 0;
   fDataPointLowerBoundaries = 0;

   fGoFNChains = 5;
   fGoFNIterationsMax = 100000;
   fGoFNIterationsRun = 2000;

   flag_discrete = false;

   fPriorConstantAll = false;
   fPriorConstantValue = 0;
}

BCModel::BCModel ( const char *  name  ) 

A constructor.

Parameters:
name The name of the model

Definition at line 47 of file BCModel.cxx.

   : BCIntegrate()
{
   fNormalization = -1.;
   fDataSet = 0;
   fParameterSet = new BCParameterSet;

   fIndex = -1;
   fPValue = -1;
   fPValueNDoF = -1;
   fChi2NDoF = -1;

   fName = (char *) name;

   fDataPointUpperBoundaries = 0;
   fDataPointLowerBoundaries = 0;

   fErrorBandXY = 0;

   fGoFNChains = 5;
   fGoFNIterationsMax = 100000;
   fGoFNIterationsRun = 2000;

   flag_discrete = false;

   fPriorConstantAll = false;
   fPriorConstantValue = 0;
}

BCModel::BCModel ( const BCModel &  bcmodel  ) 

The copy constructor.

Definition at line 104 of file BCModel.cxx.

                                        : BCIntegrate(bcmodel)
{
   fIndex                           = bcmodel.fIndex;
   fName                            = bcmodel.fName;
   fModelAPriori                    = bcmodel.fModelAPriori;
   fModelAPosteriori                = bcmodel.fModelAPosteriori;
   for (int i = 0; i < int(bcmodel.fParameterSet->size()); ++i) {
      if (bcmodel.fParameterSet->at(i)) {
         fParameterSet->push_back(new BCParameter(*(bcmodel.fParameterSet->at(i))));

      }
      else
         fParameterSet->push_back(0);
   }
   if (fDataSet)
      fDataSet = bcmodel.fDataSet;
   else
      fDataSet = 0;
   if (bcmodel.fNDataPointsMinimum)
      fNDataPointsMinimum = bcmodel.fNDataPointsMinimum;
   else
      fNDataPointsMinimum = 0;
   if (bcmodel.fNDataPointsMaximum)
      fNDataPointsMaximum = bcmodel.fNDataPointsMaximum;
   else
      fNDataPointsMaximum = 0;
   fPValue                          = bcmodel.fPValue;
   fChi2NDoF                        = bcmodel.fChi2NDoF;
   fPValueNDoF                      = bcmodel.fPValueNDoF;
   flag_discrete                    = bcmodel.flag_discrete;
   fGoFNIterationsMax               = bcmodel.fGoFNIterationsMax;
   fGoFNIterationsRun               = bcmodel.fGoFNIterationsRun;
   fGoFNChains                      = bcmodel.fGoFNChains;
   for (int i = 0; i < int(bcmodel.fPriorContainer.size()); ++i) {
      if (bcmodel.fPriorContainer.at(i))
         fPriorContainer.push_back(new TNamed(*bcmodel.fPriorContainer.at(i)));
      else
         fPriorContainer.push_back(0);
   }
   fPriorConstantAll                = bcmodel.fPriorConstantAll;
   fPriorConstantValue              = bcmodel.fPriorConstantValue;
   fPriorContainerConstant          = bcmodel.fPriorContainerConstant;
   fPriorContainerInterpolate       = bcmodel.fPriorContainerInterpolate;
   fNormalization                   = bcmodel.fNormalization;
}

BCModel::~BCModel (  )  [virtual]

The default destructor.

Definition at line 151 of file BCModel.cxx.

{
   for (unsigned int i = 0; i < GetNParameters(); ++i)
      delete fPriorContainer[i];
   fPriorContainer.clear();

   delete fParameterSet;

   delete fDataPointLowerBoundaries;
   delete fDataPointUpperBoundaries;
}

Member Function Documentation

int BCModel::AddParameter ( const char *  name,
double  lowerlimit,
double  upperlimit 
)

Adds a parameter to the parameter set

Parameters:
name The name of the parameter
lowerlimit The lower limit of the parameter values
upperlimit The upper limit of the parameter values
See also:
AddParameter(BCParameter* parameter);

Definition at line 564 of file BCModel.cxx.

{
   // create new parameter
   BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit);

   int flag_ok = AddParameter(parameter);
   if (flag_ok)
      delete parameter;

   return flag_ok;
}

int BCModel::AddParameter ( BCParameter *  parameter  ) 

Adds a parameter to the model.

Parameters:
parameter A model parameter
See also:
AddParameter(const char * name, double lowerlimit, double upperlimit);

Definition at line 577 of file BCModel.cxx.

{
   // check if parameter set exists
   if (!fParameterSet) {
      BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
      return ERROR_PARAMETERSETDOESNOTEXIST;
   }

   // check if parameter with same name exists
   int flag_exists = 0;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (CompareStrings(parameter->GetName().data(), GetParameter(i)->GetName().data()) == 0)
         flag_exists = -1;

   if (flag_exists < 0) {
      BCLog::OutError(Form(
            "BCModel::AddParameter : Parameter with name %s exists already. ",
            parameter->GetName().data()));
      return ERROR_PARAMETEREXISTSALREADY;
   }

   // define index of new parameter
   unsigned int index = fParameterSet->size();
   parameter->SetIndex(index);

   // add parameter to parameter container
   fParameterSet->push_back(parameter);

   // add empty object to prior container
   fPriorContainer.push_back(0);

   // don't interpolate the prior histogram by default
   fPriorContainerInterpolate.push_back(false);

   // prior assumed to be non-constant in general case
   fPriorContainerConstant.push_back(false);

   // add parameters to integration methods
   SetParameters(fParameterSet);

   // reset results
   ResetResults();

   return 0;
}

double BCModel::APrioriProbability ( const std::vector< double > &  parameters  )  [inline]

Returns the prior probability.

Parameters:
parameters A set of parameter values
Returns:
The prior probability p(parameters)
See also:
GetPrior(std::vector<double> parameters)

Definition at line 474 of file BCModel.h.

         { return exp( this->LogAPrioriProbability(parameters) ); }

BCH1D * BCModel::CalculatePValue ( std::vector< double >  par,
bool  flag_histogram = false 
)

Definition at line 1655 of file BCModel.cxx.

{
   BCH1D * hist = 0;

   // print log
   BCLog::OutSummary("Do goodness-of-fit-test");

   // create model test
   BCGoFTest * goftest = new BCGoFTest("modeltest");

   // set this model as the model to be tested
   goftest->SetTestModel(this);

   // set the point in parameter space which is tested an initialize
   // the model testing
   if (!goftest->SetTestPoint(par))
      return 0;

   // disable the creation of histograms to save _a lot_ of memory
   // (histograms are not needed for p-value calculation)
   goftest->MCMCSetFlagFillHistograms(false);

   // set parameters of the MCMC for the GoFTest
   goftest->MCMCSetNChains(fGoFNChains);
   goftest->MCMCSetNIterationsMax(fGoFNIterationsMax);
   goftest->MCMCSetNIterationsRun(fGoFNIterationsRun);

   // get p-value
   fPValue = goftest->GetCalculatedPValue(flag_histogram);

   // get histogram
   if (flag_histogram) {
      hist = new BCH1D();
      hist->SetHistogram(goftest->GetHistogramLogProb());
   }

   // delete model test
   delete goftest;

   // return histogram
   return hist;
}

virtual double BCModel::CDF ( const std::vector< double > &  ,
int  ,
bool   
) [inline, virtual]

1dim cumulative distribution function of the probability to get the data f(x_i|param) for a single measurement, assumed to be of identical functional form for all measurements

Parameters:
parameters The parameter values at which point to compute the cdf
index The data point index starting at 0,1...N-1
lower only needed for discrete distributions! Return the CDF for the count one less than actually observed, e.g. in Poisson process, if 3 actually observed, then CDF(2) is returned

Reimplemented in BCGraphFitter, and BCHistogramFitter.

Definition at line 761 of file BCModel.h.

      {return 0.0;}

int BCModel::CheckParameters ( const std::vector< double > &  parameters  ) 

Checks if a set of parameters values is within the given range.

Parameters:
parameters A set of parameter values
Returns:
Error code (0: OK, -1 length of parameters not correct, -2 values not within range)

Definition at line 744 of file BCModel.cxx.

{
   // check if vectors are of equal size
   if (!parameters.size() == fParameterSet->size())
      return ERROR_INVALIDNUMBEROFPARAMETERS;

   // check if parameters are within limits
   for (unsigned int i = 0; i < fParameterSet->size(); i++) {
      BCParameter * modelparameter = fParameterSet->at(i);

      if (modelparameter->GetLowerLimit() > parameters.at(i)
          || modelparameter->GetUpperLimit() < parameters.at(i)) {
         BCLog::OutError(Form(
               "BCModel::CheckParameters : Parameter %s not within limits.",
               fParameterSet-> at(i)-> GetName().data()));
         return ERROR_PARAMETERNOTWITHINRANGE;
      }
   }

   return 0;
}

int BCModel::CompareStrings ( const char *  string1,
const char *  string2 
) [private]

Compares to strings

Definition at line 2389 of file BCModel.cxx.

{
   int flag_same = 0;

   if (strlen(string1) != strlen(string2))
      return -1;

   for (int i = 0; i < int(strlen(string1)); i++)
      if (string1[i] != string2[i])
         flag_same = -1;

   return flag_same;
}

void BCModel::CorrelateDataPointValues ( std::vector< double > &  x  )  [virtual]

Constrains a data point

Parameters:
x A vector of double

Definition at line 1699 of file BCModel.cxx.

{
   // ...
}

double BCModel::Eval ( const std::vector< double > &  parameters  )  [inline, virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 535 of file BCModel.h.

         { return exp( this->LogEval(parameters) ); }

double BCModel::EvalSampling ( const std::vector< double > &  parameters  )  [virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 702 of file BCModel.cxx.

{
   return SamplingFunction(parameters);
}

void BCModel::FindMode ( std::vector< double >  start = std::vector<double>(0)  ) 

Do the mode finding using a method set via SetOptimizationMethod. Default is Minuit. The mode can be extracted using the GetBestFitParameters() method.

A starting point for the mode finding can be specified for Minuit. If not specified, Minuit default will be used (center of the parameter space).

If running mode finding after the MCMC it is a good idea to specify the mode obtained from MCMC as a starting point for the Minuit minimization. MCMC will find the location of the global mode and Minuit will converge to the mode precisely. The commands are: 

         model -> MarginalizeAll();
         model -> FindMode( model -> GetBestFitParameters() );

startinf point of Minuit minimization

Definition at line 767 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("FindMode : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return;
   }

   // this implementation is CLEARLY not good we have to work on this.

   BCLog::OutSummary(Form("Model \'%s\': Finding mode using %s",GetName().data(), DumpOptimizationMethod().c_str()));

   // synchronize parameters in BCIntegrate
   SetParameters(fParameterSet);

   switch (GetOptimizationMethod()) {
      case BCIntegrate::kOptSA:
         FindModeSA(start);
         return;

      case BCIntegrate::kOptMinuit: {
         int printlevel = -1;
         if (BCLog::GetLogLevelScreen() <= BCLog::detail)
            printlevel = 0;
         if (BCLog::GetLogLevelScreen() <= BCLog::debug)
            printlevel = 1;
         BCIntegrate::FindModeMinuit(start, printlevel);
         return;
      }

      case BCIntegrate::kOptMetropolis:
          MarginalizeAll();
         return;

      default:
         BCLog::OutError(Form("BCModel::FindMode : Invalid mode finding method: %d",GetOptimizationMethod()));
         break;
   }

   return;
}

void BCModel::FindModeMinuit ( std::vector< double >  start = std::vector<double>(0),
int  printlevel = 1 
) [virtual]

Does the mode finding using Minuit. If starting point is not specified, finding will start from the center of the parameter space.

Parameters:
start point in parameter space from which the mode finding is started.
printlevel The print level.

Reimplemented from BCIntegrate.

Definition at line 809 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("FindModeMinuit : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return;
   }

   // synchronize parameters in BCIntegrate
   SetParameters(fParameterSet);
   BCIntegrate::FindModeMinuit(start, printlevel);
}

void BCModel::FixDataAxis ( unsigned int  index,
bool  fixed 
)

Definition at line 1750 of file BCModel.cxx.

{
   // check if index is within range
   if (index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::FixDataAxis : Index out of range.");
      return;
   }

   if (fDataFixedValues.size() == 0)
      fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(),
            false);

   fDataFixedValues[index] = fixed;
}

double BCModel::GetAvalueFromChi2Johnson ( TTree *  tree,
TH1D *  distribution = 0 
)

Calculate the A-value, a summary statistic. It computes the frequency that a Chi2 value determined from the data by Johnson's binning prescription is larger than a value sampled from the reference chi2 distribution. They out from one chain is used. A=1/2 provides no evidence against the null hypothesis. Compare Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:
par tree contains the samples of posterior of the parameters
par histogram filled by function with distribution of p values

Definition at line 1519 of file BCModel.cxx.

{
   // model parameters
   int nPar = (int)GetNParameters();
   std::vector<double> param(nPar);

   //parameters saved in branches should be the same
   int nParBranches = -1;
   tree->SetBranchAddress("fNParameters", &nParBranches);

   //assume all events have same number of parameters, so check only first
   tree->GetEntry(0);
   if (nParBranches != nPar) {
      BCLog::OutError(Form(
            "Cannot compute A value: number of parameters in tree (%d)"
               "doesn't match  number of parameters in model (%d)",
            nParBranches, nPar));
      return -1.;
   }

   // buffer to construct correct branchnames for parameters, e.g. "fParameter3"
   char * branchName = new char[10 + nPar];
   // set up variables filled for each sample of parameters
   // assume same order as in model
   for (int i = 0; i < (int) nPar; i++) {
      sprintf(branchName, "fParameter%d", i);
      tree->SetBranchAddress(branchName, &param[i]);
   }

   // get the p value from Johson's definition for each param from posterior
   long nEntries = tree->GetEntries();

   // RN ~ chi2 with K-1 DoF needed for comparison
   std::vector<double> randoms(nEntries);
   int K = NumberBins();
   BCMath::RandomChi2(randoms, K - 1);

   // number of Johnson chi2 values bigger than reference
   int nBigger = 0;
   for (int i = 0; i < nEntries; i++) {
      tree->GetEntry(i);
      double chi2 = GetChi2Johnson(param);
      if (distribution != 0)
         distribution->Fill(chi2);
      // compare to set of chi2 variables
      if (chi2 > randoms[i])
         nBigger++;
   }

   return nBigger / double(nEntries);
}

double BCModel::GetBestFitParameter ( unsigned int  index  ) 

Returns the value of a parameter (defined by index) at the global mode of the posterior pdf.

Parameters:
index index of the parameter.
Returns:
best fit value of the parameter or -1e+111 on error or center of the range if mode finding not yer run

Definition at line 275 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameter : Parameter index out of range, returning -1e+111.");
      return -1e+111;
   }

   if(fBestFitParameters.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameter : Mode finding not yet run, returning center of the range.");
      return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
   }

   return fBestFitParameters[index];
}

double BCModel::GetBestFitParameterError ( unsigned int  index  ) 

Returns the error on the value of a parameter (defined by index) at the global mode of the posterior pdf.

Parameters:
index index of the parameter.
Returns:
error on the best fit value of the parameter or -1 if undefined

Definition at line 291 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameterError : Parameter index out of range, returning -1.");
      return -1;
   }

   if(fBestFitParameterErrors.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameterError : Mode finding not yet run, returning -1.");
      return -1.;
   }

   if(fBestFitParameterErrors[index]<0.)
      BCLog::OutWarning("BCModel::GetBestFitParameterError : Parameter error not available, returning -1.");

   return fBestFitParameterErrors[index];
}

std::vector<double> BCModel::GetBestFitParameterErrors (  )  [inline]

Definition at line 206 of file BCModel.h.

         { return fBestFitParameterErrors; }

double BCModel::GetBestFitParameterMarginalized ( unsigned int  index  ) 

Returns the value of a particular parameter (defined by index) at the modes of the marginalized posterior pdfs.

Parameters:
index index of the parameter.
Returns:
best fit parameter or -1e+111 on error or center of the range if marginalization not yet run

Definition at line 310 of file BCModel.cxx.

{
   if(index >= GetNParameters()) {
      BCLog::OutError("BCModel::GetBestFitParameterMarginalized : Parameter index out of range, returning -1e+111.");
      return -1e+111;
   }

   if(fBestFitParametersMarginalized.size()==0) {
      BCLog::OutError("BCModel::GetBestFitParameterMarginalized : MCMC not yet run, returning center of the range.");
      return (GetParameter(index)->GetUpperLimit() + GetParameter(index)->GetLowerLimit() ) / 2.;
   }

   return fBestFitParametersMarginalized[index];
}

std::vector<double> BCModel::GetBestFitParameters (  )  [inline]

Returns the set of values of the parameters at the global mode of the posterior pdf.

Returns:
The best fit parameters

Definition at line 203 of file BCModel.h.

         { return fBestFitParameters; }

std::vector<double> BCModel::GetBestFitParametersMarginalized (  )  [inline]

Returns the set of values of the parameters at the modes of the marginalized posterior pdfs.

Returns:
best fit parameters

Definition at line 220 of file BCModel.h.

         { return fBestFitParametersMarginalized; }

double BCModel::GetChi2Johnson ( std::vector< double >  par,
const int  nBins = -1 
)

Calculate Chi2 (also called R^{B}) for arbitrary problems with binned data using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:
par Parameter set for the calculation of the likelihood
nBins how many bins to use for the data, for negative an adapted rule \ of thumb by Wald(1942) is used, with at least three bins

Definition at line 1406 of file BCModel.cxx.

{
   typedef unsigned int uint;

   // number of observations
   int n = GetNDataPoints();

   if (nBins < 0)
      nBins = NumberBins();

   // fixed width quantiles, including final point!
   std::vector<double> a;
   for (int i = 0; i <= nBins; i++)
      a.push_back(i / double(nBins));

   // determine the bin counts and fill the histogram with data using the CDF
   TH1I * hist = new TH1I("h1", "h1 title", nBins, 0., 1.);

   // discrete case requires randomization to allocate counts of bins that cover more
   // than one quantile
   if (flag_discrete) {
      // loop over observations, each may have different likelihood and CDF
      for (int j = 0; j < n; j++) {
         // actual value
            double CDFval = CDF(par, j, false);
         // for the bin just before
         double CDFlower = CDF(par, j, true);

         // what quantiles q are covered, count from q_1 to q_{nBins}
         int qMax = 1;
         for (int i = 0; i < nBins; i++) {
            if (CDFval > a[i])
               qMax = i + 1;
            else
               break;
         }
         int qMin = 1;
         for (int i = 0; i < nBins; i++) {
            if (CDFlower > a[i])
               qMin = i + 1;
            else
               break;
         }

         // simplest case: observation bin entirely contained in one quantile
         if (qMin == qMax) {
            // watch out for overflow because CDF exactly 1
            if (CDFval < 1)
               hist->Fill(CDFval);
            else
               hist->AddBinContent(qMax);

            continue; // this observation finished
         }

         // if more than quantile is covered need more work:
         // determine probabilities of this observation to go for each quantile covered
         // as follows: If each quantile has size 0.25 and the CDF(integral of likelihood)
         // for current observation gives gives 0.27, but for observation-1 we would have
         // 0.20, then 5/7 of the 7% go for first quantile and 2/7 for the second.
         // This extend to bins covering more than two quantiles
         std::vector<double> prob;
         // normalization
         double norm = 1 / double(CDFval - CDFlower);

         for (int i = 0; i < (qMax - qMin + 1); i++) {
            if (i == 0) {
               prob.push_back(norm * (a[qMin] - CDFlower));
               continue;
            }
            if (i == (qMax - qMin)) {
               prob.push_back(norm * (CDFval - a[qMax - 1]));
               continue;
            }
            // default case
            prob.push_back(norm * (a[i] - a[i - 1]));
         }
         // have distribution, use inverse-transform method
         double U = fRandom->Rndm();
         // build up the integral (CDF)
         for (uint i = 1; i < prob.size(); i++)
            prob[i] += prob[i - 1];
         // and search with linear comput. complexity
         for (uint i = 0; i < prob.size(); i++) {
            // we finally allocate the count, as center of quantile
            if (U <= prob[i]) {
               hist->Fill((a[qMin + i - 1] + a[qMin + i]) / 2.);
               break;
            }
         }
      }
   }
   else { //continuous case is simple
      for (int j = 0; j < n; j++)
            hist->Fill(CDF(par, j, false));
   }

   // calculate chi^2
   double chi2 = 0.;
   double mk, pk;
   double N = double(n);
   for (int i = 1; i <= nBins; i++) {
      mk = hist->GetBinContent(i);
      pk = a[i] - a[i - 1];
      chi2 += (mk - N * pk) * (mk - N * pk) / (N * pk);
   }

   delete hist;

   return chi2;
}

double BCModel::GetChi2NDoF (  )  [inline]

Definition at line 697 of file BCModel.h.

         { return fChi2NDoF; }

std::vector< double > BCModel::GetChi2Runs ( int  dataIndex,
int  sigmaIndex 
)

For a Gaussian problem, calculate the chi2 of the longest run of consecutive values above/below the expected values

Parameters:
dataIndex component of datapoint with the observed value
sigmaIndex component of datapoint with uncertainty

Definition at line 1390 of file BCModel.cxx.

{
   std::vector<double> x;
   return x;
}

BCDataPoint * BCModel::GetDataPoint ( unsigned int  index  ) 
Parameters:
index The index of the data point.
Returns:
The data point in the current data set at index

Definition at line 229 of file BCModel.cxx.

{
   if (fDataSet)
      return fDataSet->GetDataPoint(index);

   BCLog::OutWarning("BCModel::GetDataPoint : No data set defined.");
   return 0;
}

BCDataPoint* BCModel::GetDataPointLowerBoundaries (  )  [inline]
Returns:
The lower boundaries of possible data values.

Definition at line 121 of file BCModel.h.

         { return fDataPointLowerBoundaries; }

double BCModel::GetDataPointLowerBoundary ( unsigned int  index  )  [inline]
Parameters:
index The index of the variable.
Returns:
The lower boundary of possible data values for a particular variable.

Definition at line 132 of file BCModel.h.

         { return fDataPointLowerBoundaries -> GetValue(index); }

BCDataPoint* BCModel::GetDataPointUpperBoundaries (  )  [inline]
Returns:
The upper boundaries of possible data values.

Definition at line 126 of file BCModel.h.

         { return fDataPointUpperBoundaries; }

double BCModel::GetDataPointUpperBoundary ( unsigned int  index  )  [inline]
Parameters:
index The index of the variable.
Returns:
The upper boundary of possible data values for a particular variable.

Definition at line 138 of file BCModel.h.

         { return fDataPointUpperBoundaries -> GetValue(index); }

BCDataSet* BCModel::GetDataSet (  )  [inline]
Returns:
The data set.

Definition at line 116 of file BCModel.h.

         { return fDataSet; }

std::vector< double > BCModel::GetErrorBand ( double  level  ) 

Returns a vector of y-values at a certain probability level.

Parameters:
level The level of probability
Returns:
vector of y-values

Definition at line 339 of file BCModel.cxx.

{
   std::vector<double> errorband;

   if (!fErrorBandXY)
      return errorband;

   int nx = fErrorBandXY->GetNbinsX();
   errorband.assign(nx, 0.);

   // loop over x and y bins
   for (int ix = 1; ix <= nx; ix++) {
      TH1D * temphist = fErrorBandXY->ProjectionY("temphist", ix, ix);

      int nprobSum = 1;
      double q[1];
      double probSum[1];
      probSum[0] = level;

      temphist->GetQuantiles(nprobSum, q, probSum);

      errorband[ix - 1] = q[0];
   }

   return errorband;
}

TGraph * BCModel::GetErrorBandGraph ( double  level1,
double  level2 
)

Definition at line 367 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   // define new graph
   int nx = fErrorBandXY->GetNbinsX();

   TGraph * graph = new TGraph(2 * nx);
   graph->SetFillStyle(1001);
   graph->SetFillColor(kYellow);

   // get error bands
   std::vector<double> ymin = GetErrorBand(level1);
   std::vector<double> ymax = GetErrorBand(level2);

   for (int i = 0; i < nx; i++) {
      graph->SetPoint(i, fErrorBandXY->GetXaxis()->GetBinCenter(i + 1), ymin[i]);
      graph->SetPoint(nx + i, fErrorBandXY->GetXaxis()->GetBinCenter(nx - i), ymax[nx - i - 1]);
   }

   return graph;
}

TH2D* BCModel::GetErrorBandXY (  )  [inline]
Returns:
The 2-d histogram of the error band.

Definition at line 225 of file BCModel.h.

         { return fErrorBandXY; }

TH2D * BCModel::GetErrorBandXY_yellow ( double  level = .68,
int  nsmooth = 0 
)

Definition at line 392 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   int nx = fErrorBandXY->GetNbinsX();
   int ny = fErrorBandXY->GetNbinsY();

   // copy existing histogram
   TH2D * hist_tempxy = (TH2D*) fErrorBandXY->Clone(
         TString::Format("%s_sub_%f.2", fErrorBandXY->GetName(), level));
   hist_tempxy->Reset();
   hist_tempxy->SetFillColor(kYellow);

   // loop over x bins
   for (int ix = 1; ix < nx; ix++) {
      BCH1D * hist_temp = new BCH1D();

      TH1D * hproj = fErrorBandXY->ProjectionY("temphist", ix, ix);
      if (nsmooth > 0)
         hproj->Smooth(nsmooth);

      hist_temp->SetHistogram(hproj);

      TH1D * hist_temp_yellow = hist_temp->GetSmallestIntervalHistogram(level);

      for (int iy = 1; iy <= ny; ++iy)
         hist_tempxy->SetBinContent(ix, iy, hist_temp_yellow->GetBinContent(iy));

      delete hist_temp_yellow;
      delete hist_temp;
   }

   return hist_tempxy;
}

TGraph* BCModel::GetFitFunctionGraph (  )  [inline]

Definition at line 240 of file BCModel.h.

         { return this -> GetFitFunctionGraph(this -> GetBestFitParameters()); }

TGraph * BCModel::GetFitFunctionGraph ( const std::vector< double > &  parameters,
double  xmin,
double  xmax,
int  n = 1000 
)

Definition at line 454 of file BCModel.cxx.

{
   // define new graph
   TGraph * graph = new TGraph(n + 1);

   double dx = (xmax - xmin) / (double) n;

   // loop over x values
   for (int i = 0; i <= n; i++) {
      double x = (double) i * dx + xmin;
      std::vector<double> xvec;
      xvec.push_back(x);
      double y = FitFunction(xvec, parameters);

      xvec.clear();

      graph->SetPoint(i, x, y);
   }

   return graph;
}

TGraph * BCModel::GetFitFunctionGraph ( const std::vector< double > &  parameters  ) 

Definition at line 429 of file BCModel.cxx.

{
   if (!fErrorBandXY)
      return 0;

   // define new graph
   int nx = fErrorBandXY->GetNbinsX();
   TGraph * graph = new TGraph(nx);

   // loop over x values
   for (int i = 0; i < nx; i++) {
      double x = fErrorBandXY->GetXaxis()->GetBinCenter(i + 1);

      std::vector<double> xvec;
      xvec.push_back(x);
      double y = FitFunction(xvec, parameters);
      xvec.clear();

      graph->SetPoint(i, x, y);
   }

   return graph;
}

bool BCModel::GetFixedDataAxis ( unsigned int  index  ) 

Definition at line 1766 of file BCModel.cxx.

{
   // check if index is within range
   if (index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::GetFixedDataAxis : Index out of range.");
      return false;
   }

   return fDataFixedValues[index];
}

bool BCModel::GetFlagBoundaries (  ) 

Checks if the boundaries have been defined

Returns:
true, if the boundaries have been set, false otherwise

Definition at line 477 of file BCModel.cxx.

{
   if (!fDataPointLowerBoundaries)
      return false;

   if (!fDataPointUpperBoundaries)
      return false;

   if (fDataPointLowerBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
      return false;

   if (fDataPointUpperBoundaries->GetNValues() != fDataSet->GetDataPoint(0)->GetNValues())
      return false;

   return true;
}

int BCModel::GetIndex (  )  [inline]
Returns:
The index of the model.

Definition at line 96 of file BCModel.h.

         { return fIndex; }

BCH1D * BCModel::GetMarginalized ( BCParameter *  parameter  ) 

If MarginalizeAll method was used, the individual marginalized distributions with respect to one parameter can be retrieved using this method.

Parameters:
parameter Model parameter
Returns:
1D marginalized probability

Definition at line 950 of file BCModel.cxx.

{
   if (!parameter) {
// don't print any error message, should be done upstream
//      BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist.");
      return 0;
   }

   int index = parameter->GetIndex();
   if (fMCMCFlagsFillHistograms[index] == false) {
// don't print any error message, should be done upstream
//      BCLog::OutError(Form(
//            "BCModel::GetMarginalized : Distribuion for '%s' not filled.",
//            parameter->GetName().data()));
      return 0;
   }

   // get histogram
   TH1D * hist = MCMCGetH1Marginalized(index);
   if (!hist)
      return 0;

   // set axis labels
   hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data()));
   hist->SetXTitle(parameter->GetName().data());
   hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data()));
   hist->SetStats(kFALSE);

   // set histogram
   BCH1D * hprob = new BCH1D();
   hprob->SetHistogram(hist);

   // set best fit parameter
   double bestfit = hprob->GetMode();

   if (fBestFitParametersMarginalized.empty())
       fBestFitParametersMarginalized.assign(GetNParameters(), 0.0);

   fBestFitParametersMarginalized[index] = bestfit;

   hprob->SetGlobalMode(fBestFitParameters.at(index));

   return hprob;
}

BCH1D* BCModel::GetMarginalized ( const char *  name  )  [inline]

Definition at line 613 of file BCModel.h.

         { return this -> GetMarginalized(this -> GetParameter(name)); }

BCH2D * BCModel::GetMarginalized ( BCParameter *  parameter1,
BCParameter *  parameter2 
)

If MarginalizeAll method was used, the individual marginalized distributions with respect to two parameters can be retrieved using this method.

Parameters:
parameter1 First parameter
parameter2 Second parameter
Returns:
2D marginalized probability

Definition at line 1314 of file BCModel.cxx.

{
   int index1 = par1->GetIndex();
   int index2 = par2->GetIndex();

   if (fMCMCFlagsFillHistograms[index1] == false || fMCMCFlagsFillHistograms[index2] == false) {
// don't print any error message, should be done upstream
//      BCLog::OutError(
//            Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.",
//                  par1->GetName().data(), par2->GetName().data()));
      return 0;
   }

   if (index1 == index2) {
// don't print any error message, should be done upstream
//      BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
      return 0;
   }

   if (par1->GetRangeWidth() == 0 || par2->GetRangeWidth() == 0){
       return 0;
   }

   BCParameter * npar1 = par1;
   BCParameter * npar2 = par2;

   if (index1 > index2) {
      npar1 = par2;
      npar2 = par1;

      int itmp = index1;
      index1 = index2;
      index2 = itmp;
   }

   // get histogram
   TH2D * hist = MCMCGetH2Marginalized(index1, index2);

   if (hist == 0)
      return 0;

   BCH2D * hprob = new BCH2D();

   // set axis labels
   hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data()));
   hist->SetXTitle(Form("%s", npar1->GetName().data()));
   hist->SetYTitle(Form("%s", npar2->GetName().data()));
   hist->SetStats(kFALSE);

   double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
   hprob->SetGlobalMode(gmode);

   // set histogram
   hprob->SetHistogram(hist);

   return hprob;
}

BCH2D* BCModel::GetMarginalized ( const char *  name1,
const char *  name2 
) [inline]

Definition at line 624 of file BCModel.h.

         { return this -> GetMarginalized(this -> GetParameter(name1), this -> GetParameter(name2)); }

double BCModel::GetModelAPosterioriProbability (  )  [inline]
Returns:
The a posteriori probability.

Definition at line 106 of file BCModel.h.

         { return fModelAPosteriori; }

double BCModel::GetModelAPrioriProbability (  )  [inline]
Returns:
The a priori probability.

Definition at line 101 of file BCModel.h.

         { return fModelAPriori; }

std::string BCModel::GetName (  )  [inline]
Returns:
The name of the model.

Definition at line 91 of file BCModel.h.

         { return fName; }

int BCModel::GetNDataPoints (  ) 
Returns:
The number of data points in the current data set.

Definition at line 215 of file BCModel.cxx.

{
   int npoints = 0;
   if (fDataSet)
      npoints = fDataSet->GetNDataPoints();
   else {
      BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined.");
      return ERROR_NOEVENTS;
   }

   return npoints;
}

unsigned int BCModel::GetNDataPointsMaximum (  )  [inline]
Returns:
The maximum number of data points.

Definition at line 162 of file BCModel.h.

         { return fNDataPointsMaximum; }

unsigned int BCModel::GetNDataPointsMinimum (  )  [inline]
Returns:
The minimum number of data points.

Definition at line 157 of file BCModel.h.

         { return fNDataPointsMinimum; }

double BCModel::GetNormalization (  )  [inline]
Returns:
The normalization factor of the probability

Definition at line 111 of file BCModel.h.

         { return fNormalization; }

unsigned int BCModel::GetNParameters (  )  [inline]
Returns:
The number of parameters of the model.

Definition at line 167 of file BCModel.h.

         { return fParameterSet ? fParameterSet -> size() : 0; }

BCParameter * BCModel::GetParameter ( int  index  ) 
Parameters:
index The index of the parameter in the parameter set.
Returns:
The parameter.

Definition at line 239 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

   if (index < 0 || index >= (int)GetNParameters()) {
      BCLog::OutWarning(
            Form("BCModel::GetParameter : Parameter index %d not within range.",index));
      return 0;
   }

   return fParameterSet->at(index);
}

BCParameter * BCModel::GetParameter ( const char *  name  ) 
Parameters:
name The name of the parameter in the parameter set.
Returns:
The parameter.

Definition at line 254 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
         index = i;

   if (index < 0) {
      BCLog::OutWarning(
            Form("BCModel::GetParameter : Model %s has no parameter named '%s'",
                 GetName().data(), name));
      return 0;
   }

   return GetParameter(index);
}

BCParameterSet* BCModel::GetParameterSet (  )  [inline]
Returns:
parameter set

Definition at line 182 of file BCModel.h.

         { return fParameterSet; }

double BCModel::GetPValue (  )  [inline]
Returns:
The p-value

Definition at line 691 of file BCModel.h.

         { return fPValue; }

double BCModel::GetPvalueFromChi2 ( const std::vector< double > &  par,
int  sigma_index 
)

Calculate p-value from Chi2 distribution for Gaussian problems

Parameters:
par Parameter set for the calculation of the likelihood
sigma_index Index of the sigma/uncertainty for the data points (for data in format "x y erry" the index would be 2)

Definition at line 1373 of file BCModel.cxx.

{
   double ll = LogLikelihood(par);
   int n = GetNDataPoints();

   double sum_sigma = 0;
   for (int i = 0; i < n; i++)
      sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));

   double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);

   fPValue = TMath::Prob(chi2, n);

   return fPValue;
}

double BCModel::GetPvalueFromChi2Johnson ( std::vector< double >  par  ) 

Calculate p-value from asymptotic Chi2 distribution for arbitrary problems using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).

Parameters:
par Parameter set for the calculation of the likelihood

Definition at line 1397 of file BCModel.cxx.

{
   double chi2 = GetChi2Johnson(par);
   // look up corresponding p value
   fPValue = TMath::Prob(chi2, NumberBins() - 1);
   return fPValue;
}

double BCModel::GetPvalueFromChi2NDoF ( std::vector< double >  par,
int  sigma_index 
)

Definition at line 1572 of file BCModel.cxx.

{
   double ll = LogLikelihood(par);
   int n = GetNDataPoints();
   int npar = GetNParameters();

   double sum_sigma = 0;
   for (int i = 0; i < n; i++)
      sum_sigma += log(GetDataPoint(i)->GetValue(sigma_index));

   double chi2 = -2. * (ll + (double) n / 2. * log(2. * M_PI) + sum_sigma);

   fChi2NDoF = chi2 / double(n - npar);
   fPValueNDoF = TMath::Prob(chi2, n - npar);

   return fPValueNDoF;
}

double BCModel::GetPvalueFromKolmogorov ( const std::vector< double > &  par,
int  index 
)

Calculate p-value from Kolmogorov-Smirnov test statistic for 1D - datasets.

Parameters:
par Parameter set for the calculation of the likelihood
index Index of the data point in the BCDataSet (for data in format "x y erry" the index would be 1)

Definition at line 1591 of file BCModel.cxx.

{
   if (flag_discrete) {
      BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
         "test defined only for continuous distributions."));
      return -1.;
   }

   //calculate the ECDF from the 1D data
   std::vector<double> yData = fDataSet->GetDataComponents(index);
   TH1D * ECDF = BCMath::ECDF(yData);

   int N = GetNDataPoints();

   // calculated expected CDF for unique points only
   std::set<double> uniqueObservations;
   for (int i = 0; i < N; i++)
       uniqueObservations.insert(CDF(par, i, false));

   int nUnique = uniqueObservations.size();
   if (nUnique != ECDF->GetNbinsX() + 1) {
      BCLog::OutError(Form("BCModel::GetPvalueFromKolmogorov : "
         "Number of unique data points doesn't match (%d vs %d)", nUnique,
            ECDF->GetNbinsX() + 1));
      return -1.;
   }

   // find maximum distance
   double distMax = 0.;

   // current distance
   double dist = 0.;

   std::set<double>::const_iterator iter = uniqueObservations.begin();
   for (int iBin = 0; iBin < nUnique; ++iBin) {
      // distance between current points
      dist = TMath::Abs(*iter - ECDF->GetBinContent(iBin + 1));
      // update maximum if necessary
      distMax = TMath::Max(dist, distMax);

//      BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
//         "expected vs empirical (%f vs %f)", *iter, ECDF->GetBinContent(iBin
//            + 1)));

      // advance to next entry in the set
      ++iter;
   }

   // correct for total #points, not unique #points.
   // would need sqrt(n1*n2/(n1+n2)) if we had two experimental datasets
   double z = distMax * sqrt(N);

   fPValue = TMath::KolmogorovProb(z);

//   BCLog::OutDebug(Form("BCModel::GetPvalueFromKolmogorov : "
//      "max distance vs corrected (%f vs %f)", distMax, z));

   // clean up
   delete ECDF;

   return fPValue;
}

double BCModel::GetPValueNDoF (  )  [inline]

Definition at line 694 of file BCModel.h.

         { return fPValueNDoF; }

double BCModel::HessianMatrixElement ( BCParameter *  parameter1,
BCParameter *  parameter2,
std::vector< double >  point 
)

Calculates the matrix element of the Hessian matrix

Parameters:
parameter1 The parameter for the first derivative
parameter2 The parameter for the first derivative
Returns:
The matrix element of the Hessian matrix

Definition at line 1705 of file BCModel.cxx.

{
   // check number of entries in vector
   if (point.size() != GetNParameters()) {
      BCLog::OutError("BCModel::HessianMatrixElement : Invalid number of entries in the vector.");
      return -1;
   }

   // define steps
   double nsteps = 1e5;
   double dx1 = par1->GetRangeWidth() / nsteps;
   double dx2 = par2->GetRangeWidth() / nsteps;

   // define points at which to evaluate
   std::vector<double> xpp = point;
   std::vector<double> xpm = point;
   std::vector<double> xmp = point;
   std::vector<double> xmm = point;

   int idx1 = par1->GetIndex();
   int idx2 = par2->GetIndex();

   xpp[idx1] += dx1;
   xpp[idx2] += dx2;

   xpm[idx1] += dx1;
   xpm[idx2] -= dx2;

   xmp[idx1] -= dx1;
   xmp[idx2] += dx2;

   xmm[idx1] -= dx1;
   xmm[idx2] -= dx2;

   // calculate probability at these points
   double ppp = Likelihood(xpp);
   double ppm = Likelihood(xpm);
   double pmp = Likelihood(xmp);
   double pmm = Likelihood(xmm);

   // return derivative
   return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
}

double BCModel::Likelihood ( const std::vector< double > &  params  )  [inline]

Returns the likelihood

Parameters:
params A set of parameter values
Returns:
The likelihood

Definition at line 489 of file BCModel.h.

         { return exp( this->LogLikelihood(params) ); }

double BCModel::LogAPrioriProbability ( const std::vector< double > &  parameters  )  [virtual]

Returns natural logarithm of the prior probability. Method needs to be overloaded by the user.

Parameters:
parameters A set of parameter values
Returns:
The prior probability p(parameters)
See also:
GetPrior(std::vector<double> parameters)

Reimplemented in BCGoFTest, BCSummaryPriorModel, and BCRooInterface.

Definition at line 648 of file BCModel.cxx.

{
   double logprob = 0;

   if(!fPriorConstantAll) {

      // get number of parameters
      int npar = GetNParameters();

      // loop over all 1-d priors
      for (int i = 0; i < npar; ++i) {
         if (fPriorContainer[i]) {
            // check what type of object is stored
            TF1 * f = dynamic_cast<TF1*>(fPriorContainer[i]);
            TH1 * h = dynamic_cast<TH1*>(fPriorContainer[i]);

            if (f) // TF1
               logprob += log(f->Eval(parameters[i]));
            else if (h) { // TH1
               if(fPriorContainerInterpolate[i])
                  logprob += log(h->Interpolate(parameters[i]));
               else
                  logprob += log(h->GetBinContent(h->FindBin(parameters[i])));
            }
            else
                 BCLog::OutError(Form(
                  "BCModel::LogAPrioriProbability : Prior for parameter %s "
                  "is defined but not recodnized.",
                  GetParameter(i)->GetName().c_str())); // this should never happen
         }
         // use constant only if user has defined it
         else if (!fPriorContainerConstant[i]) {
            BCLog::OutWarning(Form(
                  "BCModel::LogAPrioriProbability: Prior for parameter %s "
                  "is undefined. Using constant prior to proceed.",
                  GetParameter(i)->GetName().c_str()));
            logprob -= log(GetParameter(i)->GetRangeWidth());
         }
      }
   }

   // add the contribution from constant priors in one step
   logprob += fPriorConstantValue;

   return logprob;
}

double BCModel::LogEval ( const std::vector< double > &  parameters  )  [virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 696 of file BCModel.cxx.

{
   return LogProbabilityNN(parameters);
}

virtual double BCModel::LogLikelihood ( const std::vector< double > &  params  )  [pure virtual]

Calculates natural logarithm of the likelihood. Method needs to be overloaded by the user.

Parameters:
params A set of parameter values
Returns:
Natural logarithm of the likelihood

Implemented in BCGoFTest, BCSummaryPriorModel, BCEfficiencyFitter, BCGraphFitter, BCHistogramFitter, BCRooInterface, BCTemplateFitter, and BCMTF.

double BCModel::LogProbability ( const std::vector< double > &  parameter  ) 

Returns natural logarithm of the a posteriori probability given a set of parameter values

Parameters:
parameters A set of parameter values
Returns:
The a posteriori probability

Definition at line 636 of file BCModel.cxx.

{
   // check if normalized
   if (fNormalization < 0. || fNormalization == 0.) {
      BCLog::OutError("BCModel::LogProbability. Normalization not available or zero.");
      return 0.;
   }

   return LogProbabilityNN(parameters) - log(fNormalization);
}

double BCModel::LogProbabilityNN ( const std::vector< double > &  parameter  ) 

Returns the natural logarithm of likelihood times prior probability given a set of parameter values

Parameters:
parameters A set of parameter values
Returns:
The likelihood times prior probability

Definition at line 624 of file BCModel.cxx.

{
   // add log of likelihood
   double logprob = LogLikelihood(parameters);

   // add log of prior probability
   logprob += LogAPrioriProbability(parameters);

   return logprob;
}

int BCModel::MarginalizeAll (  ) 

Marginalize all probabilities wrt. single parameters and all combinations of two parameters. The individual distributions can be retrieved using the GetMarginalized method.

Returns:
Total number of marginalized distributions

Definition at line 903 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("MarginalizeAll : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return 0;
   }

   BCLog::OutSummary(Form("Running MCMC for model \'%s\'",GetName().data()));

   // prepare function fitting
   double dx = 0.;
   double dy = 0.;

   if (fFitFunctionIndexX >= 0) {
      dx = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX)
            - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX))
            / (double) fErrorBandNbinsX;

      dy = (fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY)
            - fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY))
            / (double) fErrorBandNbinsY;

      fErrorBandXY
            = new TH2D(TString::Format("errorbandxy_%d", BCLog::GetHIndex()), "",
                  fErrorBandNbinsX,
                  fDataPointLowerBoundaries->GetValue(fFitFunctionIndexX) - .5 * dx,
                  fDataPointUpperBoundaries->GetValue(fFitFunctionIndexX) + .5 * dx,
                  fErrorBandNbinsY,
                  fDataPointLowerBoundaries->GetValue(fFitFunctionIndexY) - .5 * dy,
                  fDataPointUpperBoundaries->GetValue(fFitFunctionIndexY) + .5 * dy);
      fErrorBandXY->SetStats(kFALSE);

      for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
         for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
            fErrorBandXY->SetBinContent(ix, iy, 0.);
   }

   // run the Markov chains
   MCMCMetropolis();

   // store the mode found by the chains
   StoreMode();

   return 1;
}

double BCModel::Normalize (  ) 

Integrates over the un-normalized probability and updates fNormalization.

Definition at line 717 of file BCModel.cxx.

{
   if(fParameterSet->size()<1) {
      BCLog::OutError(Form("Normalize : No parameters defined in model \'%s\'. Aborting.",GetName().data()));
      return -1.;
   }

   BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",GetName().data()));

   unsigned int n = GetNvar();

   // initialize BCIntegrate if not done already
   if (n == 0) {
      SetParameters(fParameterSet);
      n = GetNvar();
   }

   // integrate and get best fit parameters
   // maybe we have to remove the mode finding from here in the future
   fNormalization = Integrate();

   BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization));

   return fNormalization;
}

int BCModel::NumberBins (  )  [inline, private]

rule of thumb for good number of bins (Wald1942, Johnson2004) to group observations updated so minimum is three bins (for 1-5 observations)!

Parameters:
 

Definition at line 868 of file BCModel.h.

         { return (int) exp( 0.4 * log( (float)GetNDataPoints() ) ) + 2; }

BCModel & BCModel::operator= ( const BCModel &  bcmodel  ) 

Defaut assignment operator

Definition at line 164 of file BCModel.cxx.

{
   BCIntegrate::operator=(bcmodel);
   fIndex                           = bcmodel.fIndex;
   fName                            = bcmodel.fName;
   fModelAPriori                    = bcmodel.fModelAPriori;
   fModelAPosteriori                = bcmodel.fModelAPosteriori;
   for (int i = 0; i < int(bcmodel.fParameterSet->size()); ++i) {
      if (bcmodel.fParameterSet->at(i)) {
         fParameterSet->push_back(new BCParameter(*(bcmodel.fParameterSet->at(i))));

      }
      else
         fParameterSet->push_back(0);
   }
   if (fDataSet)
      fDataSet = bcmodel.fDataSet;
   else
      fDataSet = 0;
   if (bcmodel.fNDataPointsMinimum)
      fNDataPointsMinimum = bcmodel.fNDataPointsMinimum;
   else
      fNDataPointsMinimum = 0;
   if (bcmodel.fNDataPointsMaximum)
      fNDataPointsMaximum = bcmodel.fNDataPointsMaximum;
   else
      fNDataPointsMaximum = 0;
   fPValue                          = bcmodel.fPValue;
   fChi2NDoF                        = bcmodel.fChi2NDoF;
   fPValueNDoF                      = bcmodel.fPValueNDoF;
   flag_discrete                    = bcmodel.flag_discrete;
   fGoFNIterationsMax               = bcmodel.fGoFNIterationsMax;
   fGoFNIterationsRun               = bcmodel.fGoFNIterationsRun;
   fGoFNChains                      = bcmodel.fGoFNChains;
   for (int i = 0; i < int(bcmodel.fPriorContainer.size()); ++i) {
      if (bcmodel.fPriorContainer.at(i))
         fPriorContainer.push_back(new TNamed(*bcmodel.fPriorContainer.at(i)));
      else
         fPriorContainer.push_back(0);
   }
   fPriorConstantAll                = bcmodel.fPriorConstantAll;
   fPriorConstantValue              = bcmodel.fPriorConstantValue;
   fPriorContainerConstant          = bcmodel.fPriorContainerConstant;
   fPriorContainerInterpolate       = bcmodel.fPriorContainerInterpolate;
   fNormalization                   = bcmodel.fNormalization;

   // return this
   return *this;
}

int BCModel::PrintAllMarginalized ( const char *  file,
unsigned int  hdiv = 1,
unsigned int  ndiv = 1 
)

Definition at line 1134 of file BCModel.cxx.

{
   if (!fMCMCFlagFillHistograms) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled.");
      return 0;
   }

   int npar = GetNParameters();

   if (fMCMCH1Marginalized.size() == 0 || (fMCMCH2Marginalized.size() == 0
         && npar > 1)) {
      BCLog::OutError(
            "BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   // if there's only one parameter, we just want to call Print()
   if (fMCMCH1Marginalized.size() == 1 && fMCMCH2Marginalized.size() == 0) {
      BCParameter * a = GetParameter(0);
      if (GetMarginalized(a))
         GetMarginalized(a)->Print(file);
      return 1;
   }

   int c_width = 600; // default canvas width
   int c_height = 600; // default canvas height

   int type = 112; // landscape

   if (hdiv > vdiv) {
      if (hdiv > 3) {
         c_width = 1000;
         c_height = 700;
      }
      else {
         c_width = 800;
         c_height = 600;
      }
   }
   else if (hdiv < vdiv) {
      if (hdiv > 3) {
         c_height = 1000;
         c_width = 700;
      }
      else {
         c_height = 800;
         c_width = 600;
      }
      type = 111;
   }

   // calculate number of plots
   int nplots2d = npar * (npar - 1) / 2;
   int nplots = npar + nplots2d;

   // give out warning if too many plots
   BCLog::OutSummary(
         Form(
               "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
               npar, nplots2d, nplots, file));
   if (nplots > 100)
      BCLog::OutDetail("This can take a while...");

   // setup the canvas and postscript file
   TCanvas * c = new TCanvas("c", "canvas", c_width, c_height);

   TPostScript * ps = new TPostScript(file, type);

   if (type == 112)
      ps->Range(24, 16);
   else
      ps->Range(16, 24);

   // draw all 1D distributions
   ps->NewPage();
   c->cd();
   c->Clear();
   c->Divide(hdiv, vdiv);

   int n = 0;
   for (int i = 0; i < npar; i++) {
      // get corresponding parameter
      BCParameter * a = GetParameter(i);

      // check if histogram exists
      if (!GetMarginalized(a))
         continue;

      // check if histogram is filled
      if (GetMarginalized(a)->GetHistogram()->Integral() <= 0)
         continue;

      // if current page is full, switch to new page
      if (i != 0 && i % (hdiv * vdiv) == 0) {
         c->Update();
         ps->NewPage();
         c->cd();
         c->Clear();
         c->Divide(hdiv, vdiv);
      }

      // go to next pad
      c->cd(i % (hdiv * vdiv) + 1);

      // just draw a line for a delta prior
      if (a->GetRangeWidth() == 0)
          GetMarginalized(a)->Draw(4, a->GetLowerLimit());
      else
          GetMarginalized(a)->Draw();

      n++;

      if (n % 100 == 0)
         BCLog::OutDetail(Form(" --> %d plots done", n));
   }

   if (n > 0) {
      c->Update();
   }

   // check how many 2D plots are actually drawn, despite no histogram filling or delta prior
   int k = 0;
   for (int i = 0; i < npar - 1; i++) {
      for (int j = i + 1; j < npar; j++) {
         // get corresponding parameters
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);

         // check if histogram exists, or skip if one par has a delta prior
         if (!GetMarginalized(a, b))
            continue;

         // check if histogram is filled
         if (GetMarginalized(a, b)->GetHistogram()->Integral() <= 0)
            continue;

         // if current page is full, switch to new page, but only if there is data to plot
         if ((k != 0 && k % (hdiv * vdiv) == 0) || k == 0) {
            c->Update();
            ps->NewPage();
            c->cd();
            c->Clear();
            c->Divide(hdiv, vdiv);
         }

         // go to next pad
         c->cd(k % (hdiv * vdiv) + 1);

         double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
         double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
         if (deltaa <= 1e-7 * meana)
            continue;

         double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
         double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
         if (deltab <= 1e-7 * meanb)
            continue;

         GetMarginalized(a, b)->Draw(52);
         k++;

         if ((n + k) % 100 == 0)
            BCLog::OutDetail(Form(" --> %d plots done", n + k));
      }
   }

   if ((n + k) > 100 && (n + k) % 100 != 0)
      BCLog::OutDetail(Form(" --> %d plots done", n + k));

   c->Update();
   ps->Close();

   delete c;
   delete ps;

   // return total number of drawn histograms
   return n + k;
}

int BCModel::PrintAllMarginalized1D ( const char *  filebase  ) 

Definition at line 1081 of file BCModel.cxx.

{
   if (fMCMCH1Marginalized.size() == 0) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   int n = GetNParameters();
   for (int i = 0; i < n; i++) {
      BCParameter * a = GetParameter(i);
      if (GetMarginalized(a))
         GetMarginalized(a)->Print(Form("%s_1D_%s.ps", filebase, a->GetName().data()));
   }

   return n;
}

int BCModel::PrintAllMarginalized2D ( const char *  filebase  ) 

Definition at line 1099 of file BCModel.cxx.

{
   if (fMCMCH2Marginalized.size() == 0) {
      BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
      return 0;
   }

   int k = 0;
   int n = GetNParameters();
   for (int i = 0; i < n - 1; i++) {
      for (int j = i + 1; j < n; j++) {
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);

         double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
         double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
         if (deltaa <= 1e-7 * meana)
            continue;

         double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
         double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
         if (deltab <= 1e-7 * meanb)
            continue;

         if (GetMarginalized(a, b))
            GetMarginalized(a, b)->Print(
                  Form("%s_2D_%s_%s.ps",filebase, a->GetName().data(), b->GetName().data()));
         k++;
      }
   }

   return k;
}

void BCModel::PrintHessianMatrix ( std::vector< double >  parameters  ) 

Prints matrix elements of the Hessian matrix

Parameters:
parameters The parameter values at which point to evaluate the matrix

Definition at line 2352 of file BCModel.cxx.

{
   // check number of entries in vector
   if (parameters.size() != GetNParameters()) {
      BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
      return;
   }

   // print to screen
   std::cout << std::endl << " Hessian matrix elements: " << std::endl << " Point: ";

   for (int i = 0; i < int(parameters.size()); i++)
      std::cout << parameters.at(i) << " ";
   std::cout << std::endl;

   // loop over all parameter pairs
   for (unsigned int i = 0; i < GetNParameters(); i++)
      for (unsigned int j = 0; j < i; j++) {
         // calculate Hessian matrix element
         double hessianmatrixelement = HessianMatrixElement(
               fParameterSet->at(i), fParameterSet->at(j), parameters);

         // print to screen
         std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl;
      }
}

void BCModel::PrintResults ( const char *  file  ) 

Prints a summary of the Markov Chain Monte Carlo to a file.

Definition at line 2147 of file BCModel.cxx.

{
   // print summary of Markov Chain Monte Carlo

   // open file
   ofstream ofi(file);

   // check if file is open
   if (!ofi.is_open()) {
      std::cerr << "Couldn't open file " << file << std::endl;
      return;
   }

   // number of parameters and chains
   int npar = MCMCGetNParameters();
   int nchains = MCMCGetNChains();

   // check convergence
   bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0;

   ofi << std::endl 
          << " -----------------------------------------------------" << std::endl
          << " Summary" << std::endl
       << " -----------------------------------------------------" << std::endl
       << std::endl;

   ofi << " Model summary" << std::endl << " =============" << std::endl
       << " Model: " << fName.data() << std::endl
       << " Number of parameters: " << npar << std::endl
       << " List of Parameters and ranges:" << std::endl;
   for (int i = 0; i < npar; ++i)
      ofi << "  (" << i << ") Parameter \""
          << fParameterSet->at(i)->GetName().data() << "\"" << ": "
          << "(" << fParameterSet->at(i)->GetLowerLimit() << ", "
          << fParameterSet->at(i)->GetUpperLimit() << ")" << std::endl;
   ofi << std::endl;

   ofi << " Results of the optimization" << std::endl
       << " ===========================" << std::endl
       << " Optimization algorithm used: "
       << DumpUsedOptimizationMethod()<< std::endl;

   if (int(fBestFitParameters.size()) > 0) {
      ofi << " List of parameters and global mode:" << std::endl;
      for (int i = 0; i < npar; ++i) {
         ofi << "  (" << i << ") Parameter \""
             << fParameterSet->at(i)->GetName().data() << "\": "
             << fBestFitParameters[i];
             if (int(fBestFitParameterErrors.size()) == npar)
                if(fBestFitParameterErrors[i]>=0.)
                   ofi << " +- " << fBestFitParameterErrors[i];
         ofi << std::endl;
      }
      ofi << std::endl;
   }
   else {
      ofi << " No best fit information available." << std::endl;
      ofi << std::endl;
   }

   if (fPValue >= 0.) {
      ofi << " Results of the model test" << std::endl
          << " =========================" << std::endl
          << " p-value at global mode: " << fPValue << std::endl << std::endl;
   }

   if (fNormalization >= 0.) {
       ofi << " Results of the normalization" << std::endl
             << " ============================" << std::endl
             << " Integration method used:"
             << DumpIntegrationMethod() << std::endl;
       ofi << " Normalization factor: " << fNormalization << std::endl << std::endl;
   }

   // give warning if MCMC did not converge
   if (!flag_conv && fMCMCFlagRun)
      ofi << " WARNING: the Markov Chain did not converge!" << std::endl
          << " Be cautious using the following results!" << std::endl
          << std::endl;

   // print results of marginalization (if MCMC was run)
   if (fMCMCFlagRun && fMCMCFlagFillHistograms) {
      ofi << " Results of the marginalization" << std::endl
          << " ==============================" << std::endl
          << " List of parameters and properties of the marginalized"
          << std::endl << " distributions:" << std::endl;
      for (int i = 0; i < npar; ++i) {
         if (!fMCMCFlagsFillHistograms[i])
            continue;

         BCH1D * bch1d = GetMarginalized(fParameterSet->at(i));

         ofi << "  (" << i << ") Parameter \""
             << fParameterSet->at(i)->GetName().data() << "\":" << std::endl

             << "      Mean +- sqrt(V):                " << std::setprecision(4)
             << bch1d->GetMean() << " +- " << std::setprecision(4)
             << bch1d->GetRMS() << std::endl

                   << "      Median +- central 68% interval: "
             << std::setprecision(4) << bch1d->GetMedian() << " +  "
             << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian()
             << " - " << std::setprecision(4)
             << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl

             << "      (Marginalized) mode:            " << bch1d->GetMode() << std::endl;

         ofi << "       5% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.05) << std::endl
             << "      10% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.10) << std::endl
             << "      16% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.16) << std::endl
             << "      84% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.85) << std::endl
             << "      90% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.90) << std::endl
             << "      95% quantile:                   " << std::setprecision(4)
             << bch1d->GetQuantile(0.95) << std::endl;

         std::vector<double> v;
         v = bch1d->GetSmallestIntervals(0.68);
         int ninter = int(v.size());
         ofi << "      Smallest interval(s) containing 68% and local modes:"
             << std::endl;
         for (int j = 0; j < ninter; j += 5)
            ofi << "       (" << v[j] << ", " << v[j + 1]
                << ") (local mode at " << v[j + 3] << " with rel. height "
                << v[j + 2] << "; rel. area " << v[j + 4] << ")"
                << std::endl;
            ofi << std::endl;
      }
   }
   if (fMCMCFlagRun) {
      ofi << " Status of the MCMC" << std::endl << " ==================" << std::endl
          << " Convergence reached:                    " << (flag_conv ? "yes" : "no")
          << std::endl;

      if (flag_conv)
         ofi << " Number of iterations until convergence: "
             << MCMCGetNIterationsConvergenceGlobal() << std::endl;
      else
         ofi << " WARNING: the Markov Chain did not converge! Be\n"
             << " cautious using the following results!" << std::endl
             << std::endl;
      ofi << " Number of chains:                       " << MCMCGetNChains() << std::endl
          << " Number of iterations per chain:         " << MCMCGetNIterationsRun() << std::endl
          << " Average efficiencies:" << std::endl;

      std::vector<double> efficiencies;
      efficiencies.assign(npar, 0.);

      for (int ipar = 0; ipar < npar; ++ipar)
         for (int ichain = 0; ichain < nchains; ++ichain) {
            int index = ichain * npar + ipar;
            efficiencies[ipar] +=
                  double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index)
                  + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.;
         }

      for (int ipar = 0; ipar < npar; ++ipar)
         ofi << "  (" << ipar << ") Parameter \""
             << fParameterSet->at(ipar)->GetName().data() << "\": "
             << efficiencies.at(ipar) << "%" << std::endl;
      ofi << std::endl;
   }

    ofi << " -----------------------------------------------------" << std::endl
          << " Notation:" << std::endl
          << " Mean        : mean value of the marg. pdf" << std::endl
          << " Median      : maximum of the marg. pdf" << std::endl
          << " Marg. mode  : most probable value of the marg. pdf" << std::endl
          << " V           : Variance of the marg. pdf" << std::endl
          << " Quantiles   : most commonly used quantiles" <<std::endl
          << " -----------------------------------------------------" << std::endl
          << std::endl;

   // close file
   //   ofi.close;
}

void BCModel::PrintShortFitSummary ( int  chi2flag = 0  ) 

Prints a short summary of the fit results on the screen.

Definition at line 2329 of file BCModel.cxx.

{
   BCLog::OutSummary("---------------------------------------------------");
   BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data()));
   BCLog::OutSummary(Form("   Number of parameters:  Npar  = %i", GetNParameters()));
   BCLog::OutSummary(Form("   Number of data points: Ndata = %i", GetNDataPoints()));
   BCLog::OutSummary("   Number of degrees of freedom:");
   BCLog::OutSummary(Form("      NDoF = Ndata - Npar = %i", GetNDataPoints() - GetNParameters()));

   BCLog::OutSummary("   Best fit parameters (global):");
   for (unsigned int i = 0; i < GetNParameters(); ++i)
      BCLog::OutSummary(Form("      %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i)));

   BCLog::OutSummary("   Goodness-of-fit test:");
   BCLog::OutSummary(Form("      p-value = %.3g", GetPValue()));
   if (chi2flag) {
      BCLog::OutSummary(Form("      p-value corrected for NDoF = %.3g", GetPValueNDoF()));
      BCLog::OutSummary(Form("      chi2 / NDoF = %.3g", GetChi2NDoF()));
   }
   BCLog::OutSummary("---------------------------------------------------");
}

void BCModel::PrintSummary (  ) 

Prints a summary on the screen.

Definition at line 2099 of file BCModel.cxx.

{
   int nparameters = GetNParameters();

   // model summary
   std::cout << std::endl
             << "   ---------------------------------" << std::endl
             << "    Model : " << fName.data() << std::endl
             << "   ---------------------------------" << std::endl
             << "     Index                : " << fIndex << std::endl
             << "     Number of parameters : " << nparameters << std::endl << std::endl
             << "     - Parameters : " << std::endl << std::endl;

   // parameter summary
   for (int i = 0; i < nparameters; i++)
      fParameterSet->at(i)->PrintSummary();

   // best fit parameters
   if (GetBestFitParameters().size() > 0) {
      std::cout << std::endl << "     - Best fit parameters :" << std::endl << std::endl;

      for (int i = 0; i < nparameters; i++) {
         std::cout << "       " << fParameterSet->at(i)->GetName().data()
                   << " = " << GetBestFitParameter(i) << " (overall)" << std::endl;
         if ((int) fBestFitParametersMarginalized.size() == nparameters)
            std::cout << "       "
                      << fParameterSet->at(i)->GetName().data() << " = "
                      << GetBestFitParameterMarginalized(i)
                      << " (marginalized)" << std::endl;
      }
   }

   std::cout << std::endl;

   // model testing
   if (fPValue >= 0) {
      double likelihood = Likelihood(GetBestFitParameters());
      std::cout << "   - Model testing:" << std::endl << std::endl
                << "       p(data|lambda*) = " << likelihood << std::endl
                << "       p-value         = " << fPValue << std::endl << std::endl;
   }

   // normalization
   if (fNormalization > 0)
      std::cout << "     Normalization        : " << fNormalization << std::endl;
}

double BCModel::Probability ( const std::vector< double > &  parameter  )  [inline]

Returns the a posteriori probability given a set of parameter values

Parameters:
parameters A set of parameter values
Returns:
The a posteriori probability

Definition at line 517 of file BCModel.h.

         { return exp( this->LogProbability(parameter) ); }

double BCModel::ProbabilityNN ( const std::vector< double > &  params  )  [inline]

Returns the likelihood times prior probability given a set of parameter values

Parameters:
params A set of parameter values
Returns:
The likelihood times prior probability

Definition at line 503 of file BCModel.h.

         { return exp( this->LogProbabilityNN(params) ); }

int BCModel::ReadErrorBandFromFile ( const char *  file  ) 

Read

Definition at line 1054 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if (!froot->IsOpen()) {
      BCLog::OutError(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.", file));
      return 0;
   }

   int r = 0;

   TH2D * h2 = (TH2D*) froot->Get("errorbandxy");
   if (h2) {
      h2->SetDirectory(0);
      h2->SetName(TString::Format("errorbandxy_%d", BCLog::GetHIndex()));
      SetErrorBandHisto(h2);
      r = 1;
   }
   else
      BCLog::OutWarning(
            Form("BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",file));

   froot->Close();

   return r;
}

int BCModel::ReadMarginalizedFromFile ( const char *  file  ) 

Read

Definition at line 996 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if (!froot->IsOpen()) {
      BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.", file));
      return 0;
   }

   // We reset the MCMCEngine here for the moment.
   // In the future maybe we only want to do this if the engine
   // wans't initialized at all or when there were some changes
   // in the model.
   // But maybe we want reset everything since we're overwriting
   // the marginalized distributions anyway.
   MCMCInitialize();

   int k = 0;
   int n = GetNParameters();
   for (int i = 0; i < n; i++) {
      BCParameter * a = GetParameter(i);
      TKey * key = froot->GetKey(TString::Format("hist_%s_%s",GetName().data(), a->GetName().data()));
      if (key) {
         TH1D * h1 = (TH1D*) key->ReadObjectAny(TH1D::Class());
         h1->SetDirectory(0);
         if (SetMarginalized(i, h1))
            k++;
      }
      else
         BCLog::OutWarning(
               Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.",
                     GetName().data(), a->GetName().data(), file));
   }

   for (int i = 0; i < n - 1; i++) {
      for (int j = i + 1; j < n; j++) {
         BCParameter * a = GetParameter(i);
         BCParameter * b = GetParameter(j);
         TKey * key = froot->GetKey(
               TString::Format("hist_%s_%s_%s",GetName().data(), a->GetName().data(),b->GetName().data()));
         if (key) {
            TH2D * h2 = (TH2D*) key->ReadObjectAny(TH2D::Class());
            h2->SetDirectory(0);
            if (SetMarginalized(i, j, h2))
               k++;
         }
         else
            BCLog::OutWarning(
                  Form("BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.",
                        GetName().data(), a->GetName().data(), b->GetName().data(), file));
      }
   }

   froot->Close();

   return k;
}

int BCModel::ReadMode ( const char *  file  ) 

Read mode from file created by WriteMode() call

Definition at line 864 of file BCModel.cxx.

{
   ifstream ifi(file);
   if (!ifi.is_open()) {
      BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.", file));
      return 0;
   }

   int npar = fParameterSet->size();
   std::vector<double> mode;

   int i = 0;
   while (i < npar && !ifi.eof()) {
      double a;
      ifi >> a;
      mode.push_back(a);
      i++;
   }

   if (i < npar) {
      BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.", file));
      BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.", npar, i));
      return 0;
   }

   BCLog::OutSummary(
         Form("#  Read in best fit parameters (mode) for model \'%s\' from file %s:",
               fName.data(), file));
   SetMode(mode);
   for (int j = 0; j < npar; j++)
      BCLog::OutSummary(Form("#   ->Parameter %d : %s = %e", j,
            fParameterSet->at(j)->GetName().data(), fBestFitParameters[j]));

   BCLog::OutWarning("#  ! Best fit values obtained before this call will be overwritten !");

   return npar;
}

void BCModel::RecalculatePriorConstant (  )  [private]

Update the constant part of the prior.

Definition at line 2041 of file BCModel.cxx.

{
   fPriorConstantValue = 0.;

   // get number of parameters
   int npar = GetNParameters();

   int nconstant = 0;

   for (int i=0; i<npar; ++i)
      if (fPriorContainerConstant[i]) {
         // default case
         if (GetParameter(i)->GetRangeWidth() > 0)
             fPriorConstantValue -= log(GetParameter(i)->GetRangeWidth());
         // do not add infinity due to zero width for delta prior
         ++nconstant;
      }

   if (nconstant == npar)
      fPriorConstantAll = true;
   else
      fPriorConstantAll = false;
}

int BCModel::ResetResults (  ) 

Reset all results.

Returns:
An error code.

Definition at line 2088 of file BCModel.cxx.

{
   BCIntegrate::IntegrateResetResults();

   BCEngineMCMC::MCMCResetResults();

   // no error
   return 1;
}

double BCModel::SamplingFunction ( const std::vector< double > &  parameters  )  [virtual]

Sampling function used for importance sampling. Method needs to be overloaded by the user.

Parameters:
parameters A set of parameter values
Returns:
The probability density at the parameter values

Definition at line 708 of file BCModel.cxx.

{
   double probability = 1.;
   for (std::vector<BCParameter *>::const_iterator it = fParameterSet->begin(); it != fParameterSet->end(); ++it)
      probability *= 1. / ((*it)->GetUpperLimit() - (*it)->GetLowerLimit());
   return probability;
}

void BCModel::SetDataBoundaries ( unsigned int  index,
double  lowerboundary,
double  upperboundary,
bool  fixed = false 
)

Definition at line 517 of file BCModel.cxx.

{
   // check if data set exists
   if (!fDataSet) {
      BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
      return;
   }

   // check if index is within range
   if (index > fDataSet->GetDataPoint(0)->GetNValues()) {
      BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
      return;
   }

   // check if boundary data points exist
   if (!fDataPointLowerBoundaries)
      fDataPointLowerBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());

   if (!fDataPointUpperBoundaries)
      fDataPointUpperBoundaries = new BCDataPoint(fDataSet->GetDataPoint(0)->GetNValues());

   if (fDataFixedValues.size() == 0)
      fDataFixedValues.assign(fDataSet->GetDataPoint(0)->GetNValues(), false);

   // set boundaries
   fDataPointLowerBoundaries->SetValue(index, lowerboundary);
   fDataPointUpperBoundaries->SetValue(index, upperboundary);
   fDataFixedValues[index] = fixed;
}

void BCModel::SetDataSet ( BCDataSet *  dataset  )  [inline]

Sets the data set.

Parameters:
dataset A data set

Definition at line 302 of file BCModel.h.

         { fDataSet = dataset; fNormalization = -1.0; }

void BCModel::SetErrorBandContinuous ( bool  flag  ) 

Sets the error band flag to continuous function

Definition at line 548 of file BCModel.cxx.

{
   fErrorBandContinuous = flag;

   if (flag)
      return;

   // clear x-values
   fErrorBandX.clear();

   // copy data x-values
   for (unsigned int i = 0; i < fDataSet->GetNDataPoints(); ++i)
      fErrorBandX.push_back(fDataSet->GetDataPoint(i)->GetValue(fFitFunctionIndexX));
}

void BCModel::SetGoFNChains ( int  n  )  [inline]

Set number of chains in the MCMC of the p-value evaluation using MCMC

Definition at line 722 of file BCModel.h.

         { fGoFNChains=n; }

void BCModel::SetGoFNIterationsMax ( int  n  )  [inline]

Set maximum number of iterations in the MCMC pre-run of the p-value evaluation using MCMC

Definition at line 710 of file BCModel.h.

         { fGoFNIterationsMax=n; }

void BCModel::SetGoFNIterationsRun ( int  n  )  [inline]

Set number of iterations in the MCMC normal run of the p-value evaluation using MCMC

Definition at line 716 of file BCModel.h.

         { fGoFNIterationsRun=n; }

void BCModel::SetIndex ( int  index  )  [inline]

Sets the index of the model within the BCModelManager.

Parameters:
index The index of the model

Definition at line 261 of file BCModel.h.

         { fIndex = index; }

void BCModel::SetModelAPosterioriProbability ( double  probability  )  [inline]

Sets the a posteriori probability for a model.

Parameters:
model The model
probability The a posteriori probability

Definition at line 289 of file BCModel.h.

         { fModelAPosteriori = probability; }

void BCModel::SetModelAPrioriProbability ( double  probability  )  [inline]

Sets the a priori probability for a model.

Parameters:
model The model
probability The a priori probability

Definition at line 282 of file BCModel.h.

         { fModelAPriori = probability; }

void BCModel::SetName ( const char *  name  )  [inline]

Sets the name of the model.

Parameters:
name Name of the model

Definition at line 255 of file BCModel.h.

         { fName = name; }

void BCModel::SetNbins ( const char *  parname,
int  nbins 
)

Set the number of bins for the marginalized distribution of a parameter.

Parameters:
parname The name of the parameter in the parameter set
nbins Number of bins (default = 100)

Definition at line 326 of file BCModel.cxx.

{
   BCParameter * p = GetParameter(parname);
   if (!p) {
      BCLog::OutWarning(
            Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname));
      return;
   }

   BCIntegrate::SetNbins(nbins, p->GetIndex());
}

void BCModel::SetNDataPointsMaximum ( unsigned int  maximum  )  [inline]

Sets the maximum number of data points.

Definition at line 319 of file BCModel.h.

         { fNDataPointsMaximum = maximum; }

void BCModel::SetNDataPointsMinimum ( unsigned int  minimum  )  [inline]

Sets the minimum number of data points.

Definition at line 314 of file BCModel.h.

         { fNDataPointsMinimum = minimum; }

void BCModel::SetNormalization ( double  norm  )  [inline]

Sets the normalization of the likelihood. The normalization is the integral of likelihood over all parameters.

Parameters:
norm The normalization of the likelihood

Definition at line 296 of file BCModel.h.

         { fNormalization = norm; }

int BCModel::SetParameterRange ( int  index,
double  parmin,
double  parmax 
)

Set the range of a parameter

Parameters:
index The parameter index
parmin The parameter minimum
parmax The parameter maximum
Returns:
An error code.

Definition at line 2066 of file BCModel.cxx.

{
   // check index
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetParameterRange : Index out of range.");
      return 0;
   }

   // set parameter ranges in BAT
   GetParameter(index)->SetLowerLimit(parmin);
   GetParameter(index)->SetUpperLimit(parmax);
   fMCMCBoundaryMin[index] = parmin;
   fMCMCBoundaryMax[index] = parmax;

   // reset results
   ResetResults();

   // no error
   return 1;
}

void BCModel::SetParameterSet ( BCParameterSet *  parset  )  [inline]

Set all parameters of the model using a BCParameterSet container.

pointer to parameter set

Definition at line 267 of file BCModel.h.

         { fParameterSet = parset; }

int BCModel::SetPrior ( int  index,
TF1 *  f 
)

Set prior for a parameter.

Parameters:
index The parameter index
f A pointer to a function describing the prior
Returns:
An error code.

Definition at line 1778 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPrior : Index out of range.");
      return 0;
   }

   if (fPriorContainer[index])
      delete fPriorContainer[index];

   // copy function
   fPriorContainer[index] = new TF1(*f);

   fPriorContainerConstant[index] = false;

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPrior ( const char *  name,
TH1 *  h,
bool  flag = false 
)

Set prior for a parameter.

Parameters:
name parameter name
h pointer to a histogram describing the prior
flag whether or not to use linear interpolation
Returns:
An error code.

Definition at line 1954 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPrior(index, h, interpolate);
}

int BCModel::SetPrior ( const char *  name,
TF1 *  f 
)

Set prior for a parameter.

Parameters:
name The parameter name
f A pointer to a function describing the prior
Returns:
An error code.

Definition at line 1804 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPrior(index, f);
}

int BCModel::SetPrior ( int  index,
TH1 *  h,
bool  flag = false 
)

Set prior for a parameter.

Parameters:
index parameter index
h pointer to a histogram describing the prior
flag whether or not to use linear interpolation
Returns:
An error code.

Definition at line 1912 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPrior : Index out of range.");
      return 0;
   }

   // if the histogram exists
   if(h) {

      // check if histogram is 1d
      if (h->GetDimension() != 1) {
         BCLog::OutError(Form("BCModel::SetPrior : Histogram given for parameter %d is not 1D.",index)); 
         return 0;
      }

      // normalize the histogram
      h->Scale(1./h->Integral("width"));

      if(fPriorContainer[index])
         delete fPriorContainer[index];

      // set function
      fPriorContainer[index] = new TH1(*h);

      if (interpolate)
         fPriorContainerInterpolate[index] = true;

      fPriorContainerConstant[index] = false;
   }

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPriorConstant ( const char *  name  ) 

Set constant prior for this parameter

Parameters:
name the name of the parameter
Returns:
An error code

Definition at line 1993 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++) {
      if (name == GetParameter(i)->GetName()) {
         index = i;
         break;
      }
   }

   if (index == -1) {
      BCLog::OutError(Form(
            "BCModel::SetPriorConstant : parameter '%s' doesn't exist.", name));
      return 0;
   }

   return SetPriorConstant(index);
}

int BCModel::SetPriorConstant ( int  index  ) 

Set constant prior for this parameter

Parameters:
index the index of the parameter
Returns:
An error code

Definition at line 1967 of file BCModel.cxx.

{
   // check index range
   if (index < 0 && index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorConstant : Index out of range.");
      return 0;
   }

   if(fPriorContainer[index]) {
      delete fPriorContainer[index];
      fPriorContainer[index] = 0;
   }

   // set prior to a constant
   fPriorContainerConstant[index] = true;

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPriorConstantAll (  ) 

Enable caching the constant value of the prior, so LogAPrioriProbability is called only once. Note that the prior for ALL parameters is assumed to be constant. The value is computed from the parameter ranges, so make sure these are defined before this method is called.

Returns:
An error code

Definition at line 2014 of file BCModel.cxx.

{
   // get number of parameters
   int nPar = GetNParameters();

   if (nPar == 0)
      BCLog::OutWarning("BCModel::SetPriorConstantAll : No parameters defined.");

   // loop over all 1-d priors
   for (int i = 0; i < nPar; ++i) {
      if (fPriorContainer[i]) {
         delete fPriorContainer[i];
         fPriorContainer[i]=0;
      }
      fPriorContainerConstant[i] = true;
   }

   RecalculatePriorConstant();

   // reset all results
   ResetResults();

   // no error
   return 1;
}

int BCModel::SetPriorDelta ( int  index,
double  value 
)

Set delta-function prior for a parameter. Note: this sets the parameter range to the specified value. The old parameter range is lost.

Parameters:
index The parameter index
value The position of the delta function.
Returns:
An error code.

Definition at line 1817 of file BCModel.cxx.

{
   // set range to value
   SetParameterRange(index, value, value);

   // set prior
   return SetPriorConstant(index);
}

int BCModel::SetPriorDelta ( const char *  name,
double  value 
)

Set delta-function prior for a parameter. Note: this sets the parameter range to the specified value. The old parameter range is lost.

Parameters:
name The parameter name
value The position of the delta function.
Returns:
An error code.

Definition at line 1827 of file BCModel.cxx.

{
  // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
    return SetPriorDelta(index, value);
}

int BCModel::SetPriorGauss ( int  index,
double  mean,
double  sigma 
)

Set Gaussian prior for a parameter.

Parameters:
index The parameter index
mean The mean of the Gaussian
sigma The sigma of the Gaussian
Returns:
An error code.

Definition at line 1840 of file BCModel.cxx.

{
   // check index range
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorGauss : Index out of range.");
      return 0;
   }

   // create new function
   TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
                     "1./sqrt(2.*TMath::Pi())/[1] * exp(- (x-[0])*(x-[0])/2./[1]/[1])",
                     GetParameter(index)->GetLowerLimit(),
                     GetParameter(index)->GetUpperLimit());
   f->SetParameter(0, mean);
   f->SetParameter(1, sigma);

   // set prior
   return SetPrior(index, f);
}

int BCModel::SetPriorGauss ( const char *  name,
double  mean,
double  sigmadown,
double  sigmaup 
)

Set Gaussian prior for a parameter with two different widths.

Parameters:
name The parameter name
mean The mean of the Gaussian
sigmadown The sigma (down) of the Gaussian
sigmaup The sigma (up)of the Gaussian
Returns:
An error code.

Definition at line 1899 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
        index = i;

   // set prior
   return SetPriorGauss(index, mean, sigmadown, sigmaup);
}

int BCModel::SetPriorGauss ( int  index,
double  mean,
double  sigmadown,
double  sigmaup 
)

Set Gaussian prior for a parameter with two different widths.

Parameters:
index The parameter index
mean The mean of the Gaussian
sigmadown The sigma (down) of the Gaussian
sigmaup The sigma (up)of the Gaussian
Returns:
An error code.

Definition at line 1874 of file BCModel.cxx.

{
   // check index range
   if (index < 0 || index >= int(GetNParameters())) {
      BCLog::OutError("BCModel::SetPriorGauss : Index out of range.");
      return 0;
   }

   // create new function
   TF1 * f = new TF1(Form("prior_%s", GetParameter(index)->GetName().c_str()),
                     BCMath::SplitGaussian,
                     GetParameter(index)->GetLowerLimit(),
                     GetParameter(index)->GetUpperLimit(),
                     3);
   f->SetParameter(0, mean);
   f->SetParameter(1, sigmadown);
   f->SetParameter(2, sigmaup);

   // set prior
   return SetPrior(index, f);

   return 0;
}

int BCModel::SetPriorGauss ( const char *  name,
double  mean,
double  sigma 
)

Set Gaussian prior for a parameter.

Parameters:
name The parameter name
mean The mean of the Gaussian
sigma The sigma of the Gaussian
Returns:
An error code.

Definition at line 1861 of file BCModel.cxx.

{
   // find index
   int index = -1;
   for (unsigned int i = 0; i < GetNParameters(); i++)
      if (name == GetParameter(i)->GetName())
         index = i;

   // set prior
   return SetPriorGauss(index, mean, sigma);
}

void BCModel::SetSingleDataPoint ( BCDataPoint *  datapoint  ) 

Sets a single data point as data set.

Parameters:
datapoint A data point

Definition at line 495 of file BCModel.cxx.

{
   // create new data set consisting of a single data point
   BCDataSet * dataset = new BCDataSet();

   // add the data point
   dataset->AddDataPoint(datapoint);

   // set this new data set
   SetDataSet(dataset);
}

void BCModel::SetSingleDataPoint ( BCDataSet *  dataset,
unsigned int  index 
)

Definition at line 508 of file BCModel.cxx.

{
   if (index > dataset->GetNDataPoints())
      return;

   SetSingleDataPoint(dataset->GetDataPoint(index));
}

void BCModel::StoreMode (  )  [private]

save parameters at the mode after marginalization

Definition at line 2403 of file BCModel.cxx.

{
   //start with max. of posterior at first chain
   double probmax = fMCMCprobMax.at(0);
   int probmaxindex = 0;

   // loop over all chains and find the maximum point
   for (int i = 1; i < fMCMCNChains; ++i) {
      if (fMCMCprobMax.at(i) > probmax) {
         probmax = fMCMCprobMax.at(i);
         probmaxindex = i;
      }
   }

   // save if improved the log posterior
   if (fBestFitParameters.empty() || probmax > LogEval(fBestFitParameters)) {
      fBestFitParameters.assign(fMCMCNParameters, 0.0);
      for (int i = 0; i < fMCMCNParameters; ++i)
         fBestFitParameters[i] = fMCMCxMax[probmaxindex * fMCMCNParameters + i];

      SetOptimizationMethodMode(BCIntegrate::kOptMetropolis);
   }
}

BCDataPoint * BCModel::VectorToDataPoint ( const std::vector< double > &  data  )  [private]

Converts a vector of doubles into a BCDataPoint

Definition at line 2380 of file BCModel.cxx.

{
   int sizeofvector = int(data.size());
   BCDataPoint * datapoint = new BCDataPoint(sizeofvector);
   datapoint->SetValues(data);
   return datapoint;
}

void BCModel::WriteMode ( const char *  file  ) 

Write mode into file

Definition at line 822 of file BCModel.cxx.

{
   ofstream ofi(file);
   if (!ofi.is_open()) {
      std::cerr << "Couldn't open file " << file << std::endl;
      return;
   }

   int npar = fParameterSet->size();
   for (int i = 0; i < npar; i++)
      ofi << fBestFitParameters.at(i) << std::endl;

   ofi << std::endl;
   ofi << "#######################################################################"
       << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  This file was created automatically by BCModel::WriteMode() call."
       << std::endl;
   ofi << "#  It can be read in by call to BCModel::ReadMode()." << std::endl;
   ofi << "#  Do not modify it unless you know what you're doing." << std::endl;
   ofi << "#" << std::endl;
   ofi << "#######################################################################"
       << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  Best fit parameters (mode) for model:" << std::endl;
   ofi << "#  \'" << fName.data() << "\'" << std::endl;
   ofi << "#" << std::endl;
   ofi << "#  Number of parameters: " << npar << std::endl;
   ofi << "#  Parameters ordered as above:" << std::endl;

   for (int i = 0; i < npar; i++) {
      ofi << "#     " << i << ": ";
      ofi << fParameterSet->at(i)->GetName().data() << " = ";
      ofi << fBestFitParameters.at(i) << std::endl;
   }

   ofi << "#" << std::endl;
   ofi << "########################################################################"
       << std::endl;
}

Member Data Documentation

double BCModel::fChi2NDoF [protected]

Definition at line 810 of file BCModel.h.

BCDataSet* BCModel::fDataSet [protected]

A data set

Definition at line 796 of file BCModel.h.

int BCModel::fGoFNChains [protected]

Number of chains in the MCMC of the p-value evaluation using MCMC

Definition at line 830 of file BCModel.h.

int BCModel::fGoFNIterationsMax [protected]

Maximum number of iterations in the MCMC pre-run of the p-value evaluation using MCMC

Definition at line 820 of file BCModel.h.

int BCModel::fGoFNIterationsRun [protected]

Number of iterations in the MCMC normal run of the p-value evaluation using MCMC

Definition at line 825 of file BCModel.h.

int BCModel::fIndex [protected]

Index of the model.

Definition at line 776 of file BCModel.h.

bool BCModel::flag_discrete [protected]

true for a discrete probability, false for continuous pdf

Definition at line 815 of file BCModel.h.

double BCModel::fModelAPosteriori [protected]

The model a posteriori probability.

Definition at line 788 of file BCModel.h.

double BCModel::fModelAPriori [protected]

The model prior probability.

Definition at line 784 of file BCModel.h.

std::string BCModel::fName [protected]

Name of the model.

Definition at line 780 of file BCModel.h.

unsigned int BCModel::fNDataPointsMaximum [protected]

Maximum number of data points

Definition at line 804 of file BCModel.h.

unsigned int BCModel::fNDataPointsMinimum [protected]

Minimum number of data points

Definition at line 800 of file BCModel.h.

double BCModel::fNormalization [private]

The posterior normalization (evidence).

Definition at line 858 of file BCModel.h.

BCParameterSet* BCModel::fParameterSet [protected]

A model parameter container.

Definition at line 792 of file BCModel.h.

bool BCModel::fPriorConstantAll [protected]

Flag to indicate that all parameters have constant prior.

Definition at line 839 of file BCModel.h.

double BCModel::fPriorConstantValue [protected]

The value of the product of constant priors of individual parameters.

Definition at line 844 of file BCModel.h.

std::vector<TNamed*> BCModel::fPriorContainer [protected]

A vector of prior functions/histograms/graphs.

Definition at line 835 of file BCModel.h.

std::vector<bool> BCModel::fPriorContainerConstant [protected]

List the parameters whose prior is constant

Definition at line 848 of file BCModel.h.

std::vector<bool> BCModel::fPriorContainerInterpolate [protected]

List the parameters for which the histogram prior should be interpolated

Definition at line 852 of file BCModel.h.

double BCModel::fPValue [protected]

The p-value

Definition at line 808 of file BCModel.h.

double BCModel::fPValueNDoF [protected]

Definition at line 811 of file BCModel.h.

The documentation for this class was generated from the following files:
Generated by  doxygen 1.7.1