BCModel Class Reference

#include <BCModel.h>

Inheritance diagram for BCModel:

[legend]

Collaboration diagram for BCModel:

[legend]

List of all members.

---

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 conditional probabilities given those parameters. The methods which implement the prior and conditional probabilities have to be overloaded by the user in the user defined model class which will inherit from this class.

Definition at line 50 of file BCModel.h.

Public Member Functions
Member functions (miscellaneous methods)
int	AddParameter (BCParameter *parameter)
int	AddParameter (const char *name, double lowerlimit, double upperlimit)
double	APrioriProbability (std::vector< double > parameters)
BCH1D *	CalculatePValue (std::vector< double > par, bool flag_histogram=false)
int	CheckParameters (std::vector< double > parameters)
double	ConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters)
virtual void	CorrelateDataPointValues (std::vector< double > &x)
virtual void	DefineParameters ()
double	Eval (std::vector< double > parameters)
double	EvalSampling (std::vector< double > parameters)
void	FindMode ()
void	FixDataAxis (unsigned int index, bool fixed)
BCH2D *	GetMarginalized (const char *name1, const char *name2)
BCH2D *	GetMarginalized (BCParameter *parameter1, BCParameter *parameter2)
BCH1D *	GetMarginalized (const char *name)
BCH1D *	GetMarginalized (BCParameter *parameter)
double	GetPValue ()
double	GetPvalueFromChi2 (std::vector< double > par, int sigma_index)
double	HessianMatrixElement (BCParameter *parameter1, BCParameter *parameter2, std::vector< double > point)
double	Likelihood (std::vector< double > parameter)
virtual double	LogAPrioriProbability (std::vector< double > parameters)
virtual double	LogConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters)
double	LogEval (std::vector< double > parameters)
virtual double	LogLikelihood (std::vector< double > parameter)
double	LogProbability (std::vector< double > parameter)
double	LogProbabilityNN (std::vector< double > parameter)
int	MarginalizeAll ()
double	Normalize ()
int	PrintAllMarginalized (const char *file, unsigned int hdiv=1, unsigned int ndiv=1)
int	PrintAllMarginalized1D (const char *filebase)
int	PrintAllMarginalized2D (const char *filebase)
void	PrintHessianMatrix (std::vector< double > parameters)
void	PrintResults (const char *file)
void	PrintShortFitSummary ()
void	PrintSummary ()
double	Probability (std::vector< double > parameter)
double	ProbabilityNN (std::vector< double > parameter)
int	ReadErrorBandFromFile (const char *file)
int	ReadMarginalizedFromFile (const char *file)
int	ReadMode (const char *file)
virtual double	SamplingFunction (std::vector< double > parameters)
void	SetGoFNChains (int n)
void	SetGoFNIterationsMax (int n)
void	SetGoFNIterationsRun (int n)
void	WriteMode (const char *file)
Constructors and destructors
	BCModel (const char *name)
	BCModel ()
virtual	~BCModel ()
Member functions (get)
double	GetBestFitParameter (unsigned int index)
double	GetBestFitParameterMarginalized (unsigned int index)
std::vector< double >	GetBestFitParameters ()
std::vector< double >	GetBestFitParametersMarginalized ()
BCDataPoint *	GetDataPoint (unsigned int index)
BCDataPoint *	GetDataPointLowerBoundaries ()
double	GetDataPointLowerBoundary (unsigned int index)
BCDataPoint *	GetDataPointUpperBoundaries ()
double	GetDataPointUpperBoundary (unsigned int index)
BCDataSet *	GetDataSet ()
std::vector< double >	GetErrorBand (double level)
TGraph *	GetErrorBandGraph (double level1, double level2)
TH2D *	GetErrorBandXY ()
TH2D *	GetErrorBandXY_yellow (double level=.68, int nsmooth=0)
TGraph *	GetFitFunctionGraph (std::vector< double > parameters, double xmin, double xmax, int n=1000)
TGraph *	GetFitFunctionGraph ()
TGraph *	GetFitFunctionGraph (std::vector< double > parameters)
bool	GetFixedDataAxis (unsigned int index)
bool	GetFlagBoundaries ()
int	GetIndex ()
double	GetModelAPosterioriProbability ()
double	GetModelAPrioriProbability ()
std::string	GetName ()
int	GetNDataPoints ()
unsigned int	GetNDataPointsMaximum ()
unsigned int	GetNDataPointsMinimum ()
double	GetNormalization ()
unsigned int	GetNParameters ()
BCParameter *	GetParameter (const char *name)
BCParameter *	GetParameter (int index)
Member functions (set)
void	SetDataBoundaries (unsigned int index, double lowerboundary, double upperboundary, bool fixed=false)
void	SetDataSet (BCDataSet *dataset)
void	SetErrorBandContinuous (bool flag)
void	SetIndex (int index)
void	SetModelAPosterioriProbability (double probability)
void	SetModelAPrioriProbability (double probability)
void	SetName (const char *name)
void	SetNDataPointsMaximum (unsigned int maximum)
void	SetNDataPointsMinimum (unsigned int minimum)
void	SetNormalization (double norm)
void	SetSingleDataPoint (BCDataSet *dataset, unsigned int index)
void	SetSingleDataPoint (BCDataPoint *datapoint)
Protected Attributes
BCDataSet *	fDataSet
int	fGoFNChains
int	fGoFNIterationsMax
int	fGoFNIterationsRun
int	fIndex
bool	flag_ConditionalProbabilityEntry
double	fModelAPosteriori
double	fModelAPriori
std::string	fName
unsigned int	fNDataPointsMaximum
unsigned int	fNDataPointsMinimum
BCParameterSet *	fParameterSet
double	fPValue
Private Member Functions
int	CompareStrings (const char *string1, const char *string2)
BCDataPoint *	VectorToDataPoint (std::vector< double > data)
Private Attributes
double	fNormalization

---

Constructor & Destructor Documentation

BCModel::BCModel

(

)

The default constructor.

Definition at line 61 of file BCModel.cxx.

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

   fIndex = -1;
   fPValue = -1;

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

   flag_ConditionalProbabilityEntry = true;

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

BCModel::BCModel

(

const char *

name

)

A constructor.

Parameters:

name

The name of the model

Definition at line 37 of file BCModel.cxx.

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

   fIndex = -1;
   fPValue = -1;

   fName = (char *) name;
   flag_ConditionalProbabilityEntry = true;

   fDataPointUpperBoundaries = 0;
   fDataPointLowerBoundaries = 0;

   fErrorBandXY = 0;

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

BCModel::~BCModel

(

)

[virtual]

The default destructor.

Definition at line 83 of file BCModel.cxx.

{
   delete fParameterSet;

   if (fDataPointLowerBoundaries)
      delete fDataPointLowerBoundaries;

   if (fDataPointUpperBoundaries)
      delete fDataPointUpperBoundaries;
}

---

Member Function Documentation

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 420 of file BCModel.cxx.

{
   // check if parameter set exists
   if (!fParameterSet)
   {
      BCLog::Out(BCLog::error, BCLog::error,
            "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 < this -> GetNParameters(); i++)
      if (this -> CompareStrings(parameter -> GetName().data(), this -> GetParameter(i) -> GetName().data()) == 0)
         flag_exists = -1;

   if (flag_exists < 0)
   {
      BCLog::Out(BCLog::error, BCLog::error,
            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 parameters to integation methods
   this -> SetParameters(fParameterSet);

   // add parameter to MCMC
   // this -> MCMCAddParameter(parameter -> GetLowerLimit(), parameter -> GetUpperLimit());

   // re-initialize
// this -> MCMCInitialize();

   return 0;
}

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 406 of file BCModel.cxx.

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

   int flag_ok = this -> AddParameter(parameter);
   if (flag_ok)
      delete parameter;

   return flag_ok;
}

double BCModel::APrioriProbability

(

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 318 of file BCModel.h.

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

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

Definition at line 1166 of file BCModel.cxx.

{
   BCH1D * hist = 0;

   // print log
   BCLog::Out(BCLog::summary, BCLog::summary, "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;

   // 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;
}

int BCModel::CheckParameters

(

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 556 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::Out(BCLog::error, BCLog::error,
                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 1571 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;
}

double BCModel::ConditionalProbabilityEntry	(	BCDataPoint *	datapoint,
		std::vector< double >	parameters
	)			[inline]

Returns a conditional probability. Method needs to be overloaded by the user.

Parameters:

	datapoint	A data point
	parameters	A set of parameter values

Returns:

The conditional probability p(datapoint|parameters)

See also:

GetConditionalEntry(BCDataPoint* datapoint, std::vector <double> parameters)

Definition at line 378 of file BCModel.h.

         { return exp( this->LogConditionalProbabilityEntry(datapoint, parameters) ); };

void BCModel::CorrelateDataPointValues

(

std::vector< double > &

x

)

[virtual]

Constrains a data point

Parameters:

x

A vector of double

Definition at line 1208 of file BCModel.cxx.

{
   // ...
}

virtual void BCModel::DefineParameters

(

)

[inline, virtual]

Defines the parameters of the model Method needs to be overloaded by the user.

Definition at line 310 of file BCModel.h.

         { ;};

double BCModel::Eval

(

std::vector< double >

parameters

)

[inline, virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 400 of file BCModel.h.

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

double BCModel::EvalSampling

(

std::vector< double >

parameters

)

[virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 515 of file BCModel.cxx.

{
   return this -> SamplingFunction(parameters);
}

void BCModel::FindMode

(

)

Does the mode finding

Definition at line 581 of file BCModel.cxx.

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

   BCLog::Out(BCLog::summary, BCLog::summary,
      Form("Model \'%s\': Finding mode", this -> GetName().data()));

   // synchronize parameters in BCIntegrate
   this -> SetParameters(fParameterSet);

   switch(this -> GetOptimizationMethod())
   {
      case BCIntegrate::kOptSimulatedAnnealing:

         BCLog::Out(BCLog::warning, BCLog::warning, "BCModel::FindMode. Simulated annaeling not yet implemented.");
         return;
      case BCIntegrate::kOptMinuit:
         this -> FindModeMinuit();
         return;
      case BCIntegrate::kOptMetropolis:
         this -> MarginalizeAll();
         return;
      }

   BCLog::Out(BCLog::warning, BCLog::warning, Form("BCModel::FindMode. Invalid mode finding method: %d. Return.",
         this->GetOptimizationMethod()));

   return;
}

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

Definition at line 1261 of file BCModel.cxx.

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

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

   fDataFixedValues[index] = fixed;
}

double BCModel::GetBestFitParameter

(

unsigned int

index

)

[inline]

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

Parameters:

index

The index of the parameter.

Returns:

The best fit parameter.

Definition at line 172 of file BCModel.h.

         { return fBestFitParameters[index]; };

double BCModel::GetBestFitParameterMarginalized

(

unsigned int

index

)

[inline]

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

Parameters:

index

The index of the parameter.

Returns:

The best fit parameter

Definition at line 187 of file BCModel.h.

         { 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 179 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:

The best fit parameters.

Definition at line 194 of file BCModel.h.

         { return fBestFitParametersMarginalized; };

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 112 of file BCModel.cxx.

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

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

BCDataPoint* BCModel::GetDataPointLowerBoundaries

(

)

[inline]

Returns:

The lower boundaries of possible data values.

Definition at line 108 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 119 of file BCModel.h.

         { return fDataPointLowerBoundaries -> GetValue(index); };

BCDataPoint* BCModel::GetDataPointUpperBoundaries

(

)

[inline]

Returns:

The upper boundaries of possible data values.

Definition at line 113 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 125 of file BCModel.h.

         { return fDataPointUpperBoundaries -> GetValue(index); };

BCDataSet* BCModel::GetDataSet

(

)

[inline]

Returns:

The data set.

Definition at line 103 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:

The vector of y-values

Definition at line 162 of file BCModel.cxx.

{
   std::vector <double> errorband;

   if (!fErrorBandXY)
      return errorband;

   int nx = fErrorBandXY -> GetNbinsX();
   errorband.assign(nx, 0.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 192 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 = this -> GetErrorBand(level1);
   std::vector <double> ymax = this -> GetErrorBand(level2);

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

   return graph;
}

TH2D* BCModel::GetErrorBandXY

(

)

[inline]

Returns:

The 2-d histogram of the error band.

Definition at line 199 of file BCModel.h.

         { return fErrorBandXY; };

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

Definition at line 219 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	(	std::vector< double >	parameters,
		double	xmin,
		double	xmax,
		int	n = 1000
	)

Definition at line 284 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 = this -> FitFunction(xvec, parameters);

      xvec.clear();

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

   return graph;
}

TGraph* BCModel::GetFitFunctionGraph

(

)

[inline]

Definition at line 214 of file BCModel.h.

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

TGraph * BCModel::GetFitFunctionGraph

(

std::vector< double >

parameters

)

Definition at line 257 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 = this -> FitFunction(xvec, parameters);
      xvec.clear();

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

   return graph;
}

bool BCModel::GetFixedDataAxis

(

unsigned int

index

)

Definition at line 1278 of file BCModel.cxx.

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

   return fDataFixedValues.at(index);
}

bool BCModel::GetFlagBoundaries

(

)

Definition at line 309 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 83 of file BCModel.h.

         { return fIndex; };

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

Definition at line 467 of file BCModel.h.

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

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

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

Parameters:

	parameter1	First parameter
	parameter2	Second parameter

Returns:

2D marginalized probability

Definition at line 1104 of file BCModel.cxx.

{
// if(fMCMCH2Marginalized.size()==0)
// {
//    BCLog::Out(BCLog::warning, BCLog::warning,
//           "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this.");
//    return 0;
// }

   int index1 = parameter1 -> GetIndex();
   int index2 = parameter2 -> GetIndex();

   if (index1 > index2)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
             "BCModel::GetMarginalized. Index 1 has to be smaller than index 2.");
      return 0;
   }

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

   if(hist==0)
      return 0;

   BCH2D * hprob = new BCH2D();

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

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

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

   return hprob;
}

BCH1D* BCModel::GetMarginalized

(

const char *

name

)

[inline]

Definition at line 456 of file BCModel.h.

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

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 743 of file BCModel.cxx.

{
// if(fMCMCH1Marginalized.size()==0)
// {
//    BCLog::Out(BCLog::warning, BCLog::warning,
//          "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this.");
//    return 0;
// }

   int index = parameter -> GetIndex();

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

   BCH1D * hprob = new BCH1D();

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

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

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

   if (fBestFitParametersMarginalized.size() == 0)
      for (unsigned int i = 0; i < this -> GetNParameters(); i++)
         fBestFitParametersMarginalized.push_back(0.0);

   fBestFitParametersMarginalized[index] = bestfit;

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

   return hprob;
}

double BCModel::GetModelAPosterioriProbability

(

)

[inline]

Returns:

The a posteriori probability.

Definition at line 93 of file BCModel.h.

         { return fModelAPosteriori; };

double BCModel::GetModelAPrioriProbability

(

)

[inline]

Returns:

The a priori probability.

Definition at line 88 of file BCModel.h.

         { return fModelAPriori; };

std::string BCModel::GetName

(

)

[inline]

Returns:

The name of the model.

Definition at line 78 of file BCModel.h.

         { return fName; };

int BCModel::GetNDataPoints

(

)

Returns:

The number of data points in the current data set.

Definition at line 96 of file BCModel.cxx.

{
   int npoints = 0;
   if (fDataSet)
      npoints = fDataSet -> GetNDataPoints();
   else
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"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 149 of file BCModel.h.

         { return fNDataPointsMaximum; };

unsigned int BCModel::GetNDataPointsMinimum

(

)

[inline]

Returns:

The minimum number of data points.

Definition at line 144 of file BCModel.h.

         { return fNDataPointsMinimum; };

double BCModel::GetNormalization

(

)

[inline]

Returns:

The normalization factor of the probability

Definition at line 98 of file BCModel.h.

         { return fNormalization; };

unsigned int BCModel::GetNParameters

(

)

[inline]

Returns:

The number of parameters of the model.

Definition at line 154 of file BCModel.h.

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

BCParameter * BCModel::GetParameter

(

const char *

name

)

Parameters:

name

The name of the parameter in the parameter set.

Returns:

The parameter.

Definition at line 140 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

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

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

   return this->GetParameter(index);
}

BCParameter * BCModel::GetParameter

(

int

index

)

Parameters:

index

The index of the parameter in the parameter set.

Returns:

The parameter.

Definition at line 123 of file BCModel.cxx.

{
   if (!fParameterSet)
      return 0;

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

   return fParameterSet -> at(index);
}

double BCModel::GetPValue

(

)

[inline]

Definition at line 492 of file BCModel.h.

         { return fPValue; };

double BCModel::GetPvalueFromChi2	(	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 1148 of file BCModel.cxx.

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

   double sum_sigma=0;
   for (int i=0;i<n;i++)
      sum_sigma += log(this -> 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::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 1215 of file BCModel.cxx.

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

   // define steps
   double nsteps = 1e5;
   double dx1 = (parameter1 -> GetUpperLimit() - parameter1 -> GetLowerLimit()) / nsteps;
   double dx2 = (parameter2 -> GetUpperLimit() - parameter2 -> GetLowerLimit()) / 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;

   xpp.at(parameter1 -> GetIndex()) = xpp.at(parameter1 -> GetIndex()) + dx1;
   xpp.at(parameter2 -> GetIndex()) = xpp.at(parameter2 -> GetIndex()) + dx2;

   xpm.at(parameter1 -> GetIndex()) = xpm.at(parameter1 -> GetIndex()) + dx1;
   xpm.at(parameter2 -> GetIndex()) = xpm.at(parameter2 -> GetIndex()) - dx2;

   xmp.at(parameter1 -> GetIndex()) = xmp.at(parameter1 -> GetIndex()) - dx1;
   xmp.at(parameter2 -> GetIndex()) = xmp.at(parameter2 -> GetIndex()) + dx2;

   xmm.at(parameter1 -> GetIndex()) = xmm.at(parameter1 -> GetIndex()) - dx1;
   xmm.at(parameter2 -> GetIndex()) = xmm.at(parameter2 -> GetIndex()) - dx2;

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

   // calculate derivative
   double derivative = (ppp + pmm - ppm - pmp) / (4.0 * dx1 * dx2);

   return derivative;
}

double BCModel::Likelihood

(

std::vector< double >

parameter

)

[inline]

Returns the likelihood

Parameters:

parameters

A set of parameter values

Returns:

The likelihood

Definition at line 334 of file BCModel.h.

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

virtual double BCModel::LogAPrioriProbability

(

std::vector< double >

parameters

)

[inline, 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 BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.

Definition at line 327 of file BCModel.h.

         { return 1.; };

virtual double BCModel::LogConditionalProbabilityEntry	(	BCDataPoint *	datapoint,
		std::vector< double >	parameters
	)			[inline, virtual]

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

Parameters:

	datapoint	A data point
	parameters	A set of parameter values

Returns:

The conditional probability p(datapoint|parameters)

See also:

GetConditionalEntry(BCDataPoint* datapoint, std::vector <double> parameters)

Definition at line 388 of file BCModel.h.

         { flag_ConditionalProbabilityEntry = false; return 0.; };

double BCModel::LogEval

(

std::vector< double >

parameters

)

[virtual]

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 508 of file BCModel.cxx.

{
   return this -> LogProbabilityNN(parameters);
}

double BCModel::LogLikelihood

(

std::vector< double >

parameter

)

[virtual]

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

Parameters:

parameters

A set of parameter values

Returns:

Natural logarithm of the likelihood

Reimplemented in BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.

Definition at line 492 of file BCModel.cxx.

{
   double logprob = 0.;

   // add log of conditional probabilities event-by-event
   for (unsigned int i=0;i<fDataSet -> GetNDataPoints();i++)
   {
      BCDataPoint * datapoint = this -> GetDataPoint(i);
      logprob += this -> LogConditionalProbabilityEntry(datapoint, parameters);
   }

   return logprob;
}

double BCModel::LogProbability

(

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 477 of file BCModel.cxx.

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

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

double BCModel::LogProbabilityNN

(

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 464 of file BCModel.cxx.

{
   // add log of conditional probability
   double logprob = this -> LogLikelihood(parameters);

   // add log of prior probability
   logprob += this -> 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 700 of file BCModel.cxx.

{
   BCLog::Out(BCLog::summary,BCLog::summary,
      Form("Running MCMC for model \'%s\'",this->GetName().data()));

   // prepare function fitting
   double dx = 0.0;
   double dy = 0.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.0);
   }

   this -> MCMCMetropolis();
   this -> FindModeMCMC();

   // this -> PrintResults(Form("%s.txt", this -> GetName().data()));

   return 1;
}

double BCModel::Normalize

(

)

Integrates over the un-normalized probability and updates fNormalization.

Definition at line 532 of file BCModel.cxx.

{
   BCLog::Out(BCLog::summary, BCLog::summary, Form("Model \'%s\': Normalizing probability",this->GetName().data()));

   unsigned int n = this -> GetNvar();

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

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

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

   return fNormalization;
}

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

Definition at line 940 of file BCModel.cxx.

{
   if(fMCMCH1Marginalized.size()==0 || (fMCMCH2Marginalized.size()==0 && this -> GetNParameters() > 1))
   {
      BCLog::Out(BCLog::error, BCLog::error,
            "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 = this->GetParameter(0);
      this -> 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;
   }

   // get number of parameters of the model and calculate number of plots
   int npar = this -> GetNParameters();
   int nplots2d = npar * (npar-1)/2;
   int nplots = npar + nplots2d;

   // give out warning if too many plots
   BCLog::Out(BCLog::summary,BCLog::summary,
         Form("Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
               npar,nplots2d,nplots,file));
   if(nplots>100)
      BCLog::Out(BCLog::detail,BCLog::detail,"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++)
   {
      // if current page is full, swith 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);

      BCParameter * a = this->GetParameter(i);
      this -> GetMarginalized(a) -> Draw();
      n++;

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

   c->Update();

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

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

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

         BCParameter * a = this->GetParameter(i);
         BCParameter * b = this->GetParameter(j);

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

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

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

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

   if( (n+k)>100 && (n+k)%100 != 0 )
      BCLog::Out(BCLog::detail,BCLog::detail,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 881 of file BCModel.cxx.

{
   if(fMCMCH1Marginalized.size()==0)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
            "BCModel::PrintAllMarginalized1D : MarginalizeAll() has to be run prior to this to fill the distributions.");
      return 0;
   }

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

   return n;
}

int BCModel::PrintAllMarginalized2D

(

const char *

filebase

)

Definition at line 902 of file BCModel.cxx.

{
   if(fMCMCH2Marginalized.size()==0)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
            "BCModel::PrintAllMarginalized2D : MarginalizeAll() has to be run prior to this to fill the distributions.");
      return 0;
   }

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

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

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

         this -> 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 1526 of file BCModel.cxx.

{
   // check number of entries in vector
   if (parameters.size() != this -> GetNParameters())
   {
      BCLog::Out(BCLog::warning, BCLog::warning,
            Form("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 < this -> GetNParameters(); i++)
      for (unsigned int j = 0; j < i; j++)
      {
         // calculate Hessian matrix element
         double hessianmatrixelement = this -> 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 1355 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
   int npar = fParameterSet -> size();

   // check convergence
   bool flag_conv = ((this -> MCMCGetNIterationsConvergenceGlobal() > 0)?1:0);

   ofi
      << std::endl
      << " -----------------------------------------------------" << std::endl
      << " Summary of the Markov Chain Monte Carlo run" << std::endl
      << " -----------------------------------------------------" << std::endl
      << std::endl;

   if (!flag_conv)
   {
      ofi
         << " WARNING: the Markov Chain did not converge! Be" << std::endl
         << " cautious using the following results!" << 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;

   if (flag_conv)
   {
      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)
      {
         BCH1D * bch1d = this -> GetMarginalized(fParameterSet -> at(i));

         ofi
            << "  (" << i << ") Parameter \""
               << fParameterSet -> at(i) -> GetName().data() << "\"" << std::endl
            << "      Mean +- RMS:         "
               << std::setprecision(4) << bch1d -> GetMean()
               << " +- "
               << std::setprecision(4) << bch1d -> GetRMS() << std::endl
            << "      Median +- sigma:     "
               << 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
            << "      Smallest interval(s) containing 68% and local modes:" << std::endl;

         std::vector <double> v;
         v = bch1d -> GetSmallestIntervals(0.68);
         int ninter = int(v.size());

         for (int j = 0; j < ninter; j+=5)
            ofi << "       " << v.at(j) << " - " << v.at(j+1) << " (local mode at " << v.at(j+3) << " with rel. height " << v.at(j+2) << "; rel. area " << v.at(j+4) << ")" << std::endl;
      }
      ofi << std::endl;
   }

   ofi
      << " Results of the optimization" << std::endl
      << " ===========================" << std::endl
      << " Optimization algorithm used: ";
   switch(this -> GetOptimizationMethod())
   {
      case BCIntegrate::kOptSimulatedAnnealing:
         ofi << " Simulated Annealing" << std::endl;
         break;
      case BCIntegrate::kOptMinuit:
         ofi << " Minuit" << std::endl;
         break;
      case BCIntegrate::kOptMetropolis:
         ofi << " MCMC " << std::endl;
         break;
   }

   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.at(i) << 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;
   }

   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: " << this -> 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:                       " << this -> MCMCGetNChains() << std::endl
      << " Number of iterations of each chain:     " << this -> MCMCGetNIterationsMax() << std::endl
      << std::endl;

   ofi
      << " -----------------------------------------------------" << std::endl
      << std::endl
      << " Notes" << std::endl
      << " =====" << std::endl
      << " (i) Median +- sigma denotes the median, m, of the" << std::endl
      << "     marginalized distribution and the intervals from" << std::endl
      << "     m to the 16% and 84% quantiles." << std::endl
      << " -----------------------------------------------------" << std::endl;

   // close file
// ofi.close;
}

void BCModel::PrintShortFitSummary

(

)

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

Definition at line 1509 of file BCModel.cxx.

{
   BCLog::Out(BCLog::summary, BCLog::summary, "---------------------------------------------------");
   BCLog::Out(BCLog::summary, BCLog::summary, Form("Fit summary for model \'%s\':",this -> GetName().data()));
   BCLog::Out(BCLog::summary, BCLog::summary, Form("   Number of parameters = %i", this -> GetNParameters()));

   BCLog::Out(BCLog::summary, BCLog::summary, "   Best fit parameters (global):");
   for (unsigned int i = 0; i < this -> GetNParameters(); ++i)
      BCLog::Out(BCLog::summary, BCLog::summary, Form("      %s = %.2lf", this -> GetParameter(i) -> GetName().data(), this -> GetBestFitParameter(i)));

   BCLog::Out(BCLog::summary, BCLog::summary, "   Goodness-of-fit test:");
   BCLog::Out(BCLog::summary, BCLog::summary, Form("      p-value = %.2lf", this -> GetPValue()));
   BCLog::Out(BCLog::summary, BCLog::summary, "---------------------------------------------------");
}

void BCModel::PrintSummary

(

)

Prints a summary on the screen.

Definition at line 1292 of file BCModel.cxx.

{
   int nparameters = this -> 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 (this -> 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()
            << " = " << this -> GetBestFitParameter(i)
            << " (overall)" << std::endl;
         if ((int)fBestFitParametersMarginalized.size() == nparameters)
            std::cout
               << "       " << fParameterSet -> at(i) -> GetName().data()
               << " = " << this -> GetBestFitParameterMarginalized(i)
               << " (marginalized)" << std::endl;
      }
   }

   std::cout << std::endl;

   // model testing
   if (fPValue >= 0)
   {
      double likelihood = this -> Likelihood(this -> 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

(

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 362 of file BCModel.h.

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

double BCModel::ProbabilityNN

(

std::vector< double >

parameter

)

[inline]

Returns the 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 348 of file BCModel.h.

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

int BCModel::ReadErrorBandFromFile

(

const char *

file

)

Read

Definition at line 851 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if(!froot->IsOpen())
   {
      BCLog::Out(BCLog::warning, BCLog::warning, 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()));
      this->SetErrorBandHisto(h2);
      r=1;
   }
   else
      BCLog::Out(BCLog::warning, BCLog::warning,
         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 786 of file BCModel.cxx.

{
   TFile * froot = new TFile(file);
   if(!froot->IsOpen())
   {
      BCLog::Out(BCLog::error, BCLog::error,
            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.
   this -> MCMCInitialize();

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

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

{
   ifstream ifi(file);
   if(!ifi.is_open())
   {
      BCLog::Out(BCLog::error,BCLog::error,
            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::Out(BCLog::error,BCLog::error,
            Form("BCModel::ReadMode : Couldn't read mode from file %s.",file));
      BCLog::Out(BCLog::error,BCLog::error,
            Form("BCModel::ReadMode : Expected %d parameters, found %d.",npar,i));
      return 0;
   }

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

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

   return npar;
}

double BCModel::SamplingFunction

(

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 522 of file BCModel.cxx.

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

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

Definition at line 353 of file BCModel.cxx.

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

   // check if index is within range
   if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
   {
      BCLog::Out(BCLog::error, BCLog::error,
            "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 262 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 389 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]

Definition at line 510 of file BCModel.h.

         { fGoFNChains=n; };

void BCModel::SetGoFNIterationsMax

(

int

n

)

[inline]

Definition at line 498 of file BCModel.h.

         { fGoFNIterationsMax=n; };

void BCModel::SetGoFNIterationsRun

(

int

n

)

[inline]

Definition at line 504 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 235 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 249 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 242 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 229 of file BCModel.h.

         { fName = name; };

void BCModel::SetNDataPointsMaximum

(

unsigned int

maximum

)

[inline]

Sets the maximum number of data points.

Definition at line 279 of file BCModel.h.

         { fNDataPointsMaximum = maximum; };

void BCModel::SetNDataPointsMinimum

(

unsigned int

minimum

)

[inline]

Sets the minimum number of data points.

Definition at line 274 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 256 of file BCModel.h.

         { fNormalization = norm; };

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

Definition at line 343 of file BCModel.cxx.

{
   if (index < 0 || index > dataset -> GetNDataPoints())
      return;

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

void BCModel::SetSingleDataPoint

(

BCDataPoint *

datapoint

)

Sets a single data point as data set.

Parameters:

datapoint

A data point

Definition at line 328 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
   this -> SetDataSet(dataset);

}

BCDataPoint * BCModel::VectorToDataPoint

(

std::vector< double >

data

)

[private]

Converts a vector of doubles into a BCDataPoint

Definition at line 1561 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 613 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

BCDataSet* BCModel::fDataSet [protected]

A data set

Definition at line 565 of file BCModel.h.

int BCModel::fGoFNChains [protected]

Definition at line 596 of file BCModel.h.

int BCModel::fGoFNIterationsMax [protected]

Definition at line 586 of file BCModel.h.

int BCModel::fGoFNIterationsRun [protected]

Definition at line 591 of file BCModel.h.

int BCModel::fIndex [protected]

Index of the model.

Definition at line 545 of file BCModel.h.

bool BCModel::flag_ConditionalProbabilityEntry [protected]

A flag for overloading ConditionalProbabilityEntry

Definition at line 577 of file BCModel.h.

double BCModel::fModelAPosteriori [protected]

The model a posteriori probability.

Definition at line 557 of file BCModel.h.

double BCModel::fModelAPriori [protected]

The model prior probability.

Definition at line 553 of file BCModel.h.

std::string BCModel::fName [protected]

Name of the model.

Definition at line 549 of file BCModel.h.

unsigned int BCModel::fNDataPointsMaximum [protected]

Maximum number of data points

Definition at line 573 of file BCModel.h.

unsigned int BCModel::fNDataPointsMinimum [protected]

Minimum number of data points

Definition at line 569 of file BCModel.h.

double BCModel::fNormalization [private]

The Likelihood normalization.

Definition at line 610 of file BCModel.h.

BCParameterSet* BCModel::fParameterSet [protected]

A model parameter container.

Definition at line 561 of file BCModel.h.

double BCModel::fPValue [protected]

The p-value

Definition at line 581 of file BCModel.h.