BCModel Class Reference

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

#include <BCModel.h>

Inherits BCIntegrate.

Inherited by BCEfficiencyFitter, BCGoFTest, BCGraphFitter, BCHistogramFitter, and BCRooInterface.

Collaboration diagram for BCModel:
Collaboration graph
[legend]

List of all members.

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)
double Eval (std::vector< double > parameters)
double EvalSampling (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 FixDataAxis (unsigned int index, bool fixed)
double GetChi2NDoF ()
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 GetPvalueFromChi2NDoF (std::vector< double > par, int sigma_index)
double GetPValueNDoF ()
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 (int chi2flag=0)
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 GetBestFitParameterError (unsigned int index)
std::vector< double > GetBestFitParameterErrors ()
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)
BCParameterSet * GetParameterSet ()
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 SetNbins (const char *parname, int nbins)
void SetNDataPointsMaximum (unsigned int maximum)
void SetNDataPointsMinimum (unsigned int minimum)
void SetNormalization (double norm)
void SetParameterSet (BCParameterSet *parset)
void SetSingleDataPoint (BCDataSet *dataset, unsigned int index)
void SetSingleDataPoint (BCDataPoint *datapoint)

Protected Attributes

double fChi2NDoF
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
double fPValueNDoF

Private Member Functions

int CompareStrings (const char *string1, const char *string2)
BCDataPoint * VectorToDataPoint (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 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.

Constructor & Destructor Documentation

BCModel::BCModel (  ) 

The default constructor.

Definition at line 63 of file BCModel.cxx.

00063                  : BCIntegrate()
00064 {
00065    fNormalization = -1.0;
00066    fDataSet = 0;
00067    fParameterSet = new BCParameterSet();
00068 
00069    fIndex = -1;
00070    fPValue = -1;
00071    fPValueNDoF = -1;
00072    fChi2NDoF   = -1;
00073 
00074    fName = "model";
00075    fDataPointUpperBoundaries = 0;
00076    fDataPointLowerBoundaries = 0;
00077 
00078    flag_ConditionalProbabilityEntry = true;
00079 
00080    fGoFNChains = 5;
00081    fGoFNIterationsMax = 100000;
00082    fGoFNIterationsRun = 2000;
00083 }

BCModel::BCModel ( const char *  name  ) 

A constructor.

Parameters:
name The name of the model

Definition at line 37 of file BCModel.cxx.

00037                                   : BCIntegrate()
00038 {
00039    fNormalization = -1.0;
00040    fDataSet = 0;
00041    fParameterSet = new BCParameterSet;
00042 
00043    fIndex = -1;
00044    fPValue = -1;
00045    fPValueNDoF = -1;
00046    fChi2NDoF   = -1;
00047 
00048    fName = (char *) name;
00049    flag_ConditionalProbabilityEntry = true;
00050 
00051    fDataPointUpperBoundaries = 0;
00052    fDataPointLowerBoundaries = 0;
00053 
00054    fErrorBandXY = 0;
00055 
00056    fGoFNChains = 5;
00057    fGoFNIterationsMax = 100000;
00058    fGoFNIterationsRun = 2000;
00059 }

BCModel::~BCModel (  )  [virtual]

The default destructor.

Definition at line 87 of file BCModel.cxx.

00088 {
00089    delete fParameterSet;
00090 
00091    if (fDataPointLowerBoundaries)
00092       delete fDataPointLowerBoundaries;
00093 
00094    if (fDataPointUpperBoundaries)
00095       delete fDataPointUpperBoundaries;
00096 }

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

00440 {
00441    // check if parameter set exists
00442    if (!fParameterSet)
00443    {
00444       BCLog::OutError("BCModel::AddParameter : Parameter set does not exist");
00445       return ERROR_PARAMETERSETDOESNOTEXIST;
00446    }
00447 
00448    // check if parameter with same name exists
00449    int flag_exists = 0;
00450    for (unsigned int i = 0; i < this -> GetNParameters(); i++)
00451       if (this -> CompareStrings(parameter -> GetName().data(), this -> GetParameter(i) -> GetName().data()) == 0)
00452          flag_exists = -1;
00453 
00454    if (flag_exists < 0)
00455    {
00456       BCLog::OutError(
00457             Form("BCModel::AddParameter : Parameter with name %s exists already. ", parameter -> GetName().data()));
00458       return ERROR_PARAMETEREXISTSALREADY;
00459    }
00460 
00461    // define index of new parameter
00462    unsigned int index = fParameterSet -> size();
00463    parameter -> SetIndex(index);
00464 
00465    // add parameter to parameter container
00466    fParameterSet -> push_back(parameter);
00467 
00468    // add parameters to integation methods
00469    this -> SetParameters(fParameterSet);
00470 
00471    return 0;
00472 }

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

00426 {
00427    // create new parameter
00428    BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit);
00429 
00430    int flag_ok = this -> AddParameter(parameter);
00431    if (flag_ok)
00432       delete parameter;
00433 
00434    return flag_ok;
00435 }

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

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

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

Definition at line 1258 of file BCModel.cxx.

01259 {
01260    BCH1D * hist = 0;
01261 
01262    // print log
01263    BCLog::OutSummary("Do goodness-of-fit-test");
01264 
01265    // create model test
01266    BCGoFTest * goftest = new BCGoFTest("modeltest");
01267 
01268    // set this model as the model to be tested
01269    goftest -> SetTestModel(this);
01270 
01271    // set the point in parameter space which is tested an initialize
01272    // the model testing
01273    if (!goftest -> SetTestPoint(par))
01274       return 0;
01275    
01276    // disable the creation of histograms to save _a lot_ of memory
01277    // (histograms are not needed for p-value calculation)
01278    goftest->MCMCSetFlagFillHistograms(false);
01279 
01280    // set parameters of the MCMC for the GoFTest
01281    goftest -> MCMCSetNChains(fGoFNChains);
01282    goftest -> MCMCSetNIterationsMax(fGoFNIterationsMax);
01283    goftest -> MCMCSetNIterationsRun(fGoFNIterationsRun);
01284 
01285    // get p-value
01286    fPValue = goftest -> GetCalculatedPValue(flag_histogram);
01287 
01288    // get histogram
01289    if (flag_histogram)
01290    {
01291       hist = new BCH1D();
01292       hist -> SetHistogram(goftest -> GetHistogramLogProb());
01293    }
01294 
01295    // delete model test
01296    delete goftest;
01297 
01298    // return histogram
01299    return hist;
01300 }

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

00569 {
00570    // check if vectors are of equal size
00571    if (!parameters.size() == fParameterSet -> size())
00572       return  ERROR_INVALIDNUMBEROFPARAMETERS;
00573 
00574    // check if parameters are within limits
00575    for (unsigned int i = 0; i < fParameterSet -> size(); i++)
00576    {
00577       BCParameter * modelparameter = fParameterSet -> at(i);
00578 
00579       if (modelparameter -> GetLowerLimit() > parameters.at(i) ||
00580             modelparameter -> GetUpperLimit() < parameters.at(i))
00581       {
00582          BCLog::OutError(
00583                 Form("BCModel::CheckParameters : Parameter %s not within limits.", fParameterSet -> at(i) -> GetName().data()));
00584          return ERROR_PARAMETERNOTWITHINRANGE;
00585       }
00586    }
00587 
00588    return 0;
00589 }

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

Compares to strings

Definition at line 1719 of file BCModel.cxx.

01720 {
01721    int flag_same = 0;
01722 
01723    if (strlen(string1) != strlen(string2))
01724       return -1;
01725 
01726    for (int i = 0; i < int(strlen(string1)); i++)
01727       if (string1[i] != string2[i])
01728          flag_same = -1;
01729 
01730    return flag_same;
01731 }

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

00396          { 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 1304 of file BCModel.cxx.

01305 {
01306    // ...
01307 }

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

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 417 of file BCModel.h.

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

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

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 527 of file BCModel.cxx.

00528 {
00529    return this -> SamplingFunction(parameters);
00530 }

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

00594 {
00595 // this implementation is CLEARLY not good we have to work on this.
00596 
00597    BCLog::OutSummary(Form("Model \'%s\': Finding mode", this -> GetName().data()));
00598 
00599    // synchronize parameters in BCIntegrate
00600    this -> SetParameters(fParameterSet);
00601 
00602    switch(this -> GetOptimizationMethod())
00603    {
00604       case BCIntegrate::kOptSA:
00605          this -> FindModeSA(start);
00606          return;
00607 
00608       case BCIntegrate::kOptMinuit:
00609          {
00610             int printlevel = -1; 
00611             if (BCLog::GetLogLevelScreen() <= BCLog::detail) 
00612                printlevel = 0; 
00613             if (BCLog::GetLogLevelScreen() <= BCLog::debug)
00614                printlevel = 1; 
00615             this -> BCIntegrate::FindModeMinuit(start, printlevel);
00616             return;
00617          }
00618 
00619       case BCIntegrate::kOptMetropolis:
00620          this -> MarginalizeAll();
00621          return;
00622       }
00623 
00624    BCLog::OutError(
00625       Form("BCModel::FindMode : Invalid mode finding method: %d",
00626          this->GetOptimizationMethod()));
00627 
00628    return;
00629 }

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

00634 {
00635    // synchronize parameters in BCIntegrate
00636    this -> SetParameters(fParameterSet);
00637 
00638    this -> BCIntegrate::FindModeMinuit(start,printlevel);
00639 }

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

Definition at line 1358 of file BCModel.cxx.

01359 {
01360    // check if index is within range
01361    if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
01362    {
01363       BCLog::OutWarning("BCModel::FixDataAxis. Index out of range.");
01364       return;
01365    }
01366 
01367    if (fDataFixedValues.size() == 0)
01368       fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false);
01369 
01370    fDataFixedValues[index] = fixed;
01371 }

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

00178          { return fBestFitParameters[index]; };

double BCModel::GetBestFitParameterError ( unsigned int  index  )  [inline]

Definition at line 180 of file BCModel.h.

00181       { return fBestFitParameterErrors[index]; };

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

Definition at line 190 of file BCModel.h.

00191          { return fBestFitParameterErrors; };

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

00199          { 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 187 of file BCModel.h.

00188          { 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 205 of file BCModel.h.

00206          { return fBestFitParametersMarginalized; };

double BCModel::GetChi2NDoF (  )  [inline]

Definition at line 539 of file BCModel.h.

00540          { return fChi2NDoF; };

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

00117 {
00118    if (fDataSet)
00119       return fDataSet -> GetDataPoint(index);
00120 
00121    BCLog::OutWarning("BCModel::GetDataPoint. No data set defined.");
00122    return 0;
00123 }

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

Definition at line 108 of file BCModel.h.

00109          { 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.

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

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

Definition at line 113 of file BCModel.h.

00114          { 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.

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

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

Definition at line 103 of file BCModel.h.

00104          { 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 183 of file BCModel.cxx.

00184 {
00185    std::vector <double> errorband;
00186 
00187    if (!fErrorBandXY)
00188       return errorband;
00189 
00190    int nx = fErrorBandXY -> GetNbinsX();
00191    errorband.assign(nx, 0.0);
00192 
00193    // loop over x and y bins
00194    for (int ix = 1; ix <= nx; ix++)
00195    {
00196       TH1D * temphist = fErrorBandXY -> ProjectionY("temphist", ix, ix);
00197 
00198       int nprobSum = 1;
00199       double q[1];
00200       double probSum[1];
00201       probSum[0] = level;
00202 
00203       temphist -> GetQuantiles(nprobSum, q, probSum);
00204 
00205       errorband[ix-1] = q[0];
00206    }
00207 
00208    return errorband;
00209 }

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

Definition at line 213 of file BCModel.cxx.

00214 {
00215    if (!fErrorBandXY)
00216       return 0;
00217 
00218    // define new graph
00219    int nx = fErrorBandXY -> GetNbinsX();
00220 
00221    TGraph * graph = new TGraph(2 * nx);
00222    graph -> SetFillStyle(1001);
00223    graph -> SetFillColor(kYellow);
00224 
00225    // get error bands
00226    std::vector <double> ymin = this -> GetErrorBand(level1);
00227    std::vector <double> ymax = this -> GetErrorBand(level2);
00228 
00229    for (int i = 0; i < nx; i++)
00230    {
00231       graph -> SetPoint(i,      fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1), ymin.at(i));
00232       graph -> SetPoint(nx + i, fErrorBandXY -> GetXaxis() -> GetBinCenter(nx - i), ymax.at(nx - i - 1));
00233    }
00234 
00235    return graph;
00236 }

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

Definition at line 210 of file BCModel.h.

00211          { return fErrorBandXY; };

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

Definition at line 240 of file BCModel.cxx.

00241 {
00242    if (!fErrorBandXY)
00243       return 0;
00244 
00245    int nx = fErrorBandXY -> GetNbinsX();
00246    int ny = fErrorBandXY -> GetNbinsY();
00247 
00248    // copy existing histogram
00249    TH2D * hist_tempxy = (TH2D*) fErrorBandXY -> Clone(TString::Format("%s_sub_%f.2",fErrorBandXY->GetName(),level));
00250    hist_tempxy -> Reset();
00251    hist_tempxy -> SetFillColor(kYellow);
00252 
00253    // loop over x bins
00254    for (int ix = 1; ix < nx; ix++)
00255    {
00256       BCH1D * hist_temp = new BCH1D();
00257 
00258       TH1D * hproj = fErrorBandXY -> ProjectionY("temphist", ix, ix);
00259       if(nsmooth>0)
00260          hproj->Smooth(nsmooth);
00261 
00262       hist_temp -> SetHistogram(hproj);
00263 
00264       TH1D * hist_temp_yellow = hist_temp -> GetSmallestIntervalHistogram(level);
00265 
00266       for (int iy = 1; iy <= ny; ++iy)
00267          hist_tempxy -> SetBinContent(ix, iy, hist_temp_yellow -> GetBinContent(iy));
00268 
00269       delete hist_temp_yellow;
00270       delete hist_temp;
00271    }
00272 
00273    return hist_tempxy;
00274 }

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

Definition at line 305 of file BCModel.cxx.

00306 {
00307    // define new graph
00308    TGraph * graph = new TGraph(n+1);
00309 
00310    double dx = (xmax-xmin)/(double)n;
00311 
00312    // loop over x values
00313    for (int i = 0; i <= n; i++)
00314    {
00315       double x = (double)i*dx+xmin;
00316       std::vector <double> xvec;
00317       xvec.push_back(x);
00318       double y = this -> FitFunction(xvec, parameters);
00319 
00320       xvec.clear();
00321 
00322       graph -> SetPoint(i, x, y);
00323    }
00324 
00325    return graph;
00326 }

TGraph* BCModel::GetFitFunctionGraph (  )  [inline]

Definition at line 225 of file BCModel.h.

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

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

Definition at line 278 of file BCModel.cxx.

00279 {
00280    if (!fErrorBandXY)
00281       return 0;
00282 
00283    // define new graph
00284    int nx = fErrorBandXY -> GetNbinsX();
00285    TGraph * graph = new TGraph(nx);
00286 
00287    // loop over x values
00288    for (int i = 0; i < nx; i++)
00289    {
00290       double x = fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1);
00291 
00292       std::vector <double> xvec;
00293       xvec.push_back(x);
00294       double y = this -> FitFunction(xvec, parameters);
00295       xvec.clear();
00296 
00297       graph -> SetPoint(i, x, y);
00298    }
00299 
00300    return graph;
00301 }

bool BCModel::GetFixedDataAxis ( unsigned int  index  ) 

Definition at line 1375 of file BCModel.cxx.

01376 {
01377    // check if index is within range
01378    if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
01379    {
01380       BCLog::OutWarning("BCModel::GetFixedDataAxis. Index out of range.");
01381       return false;
01382    }
01383 
01384    return fDataFixedValues.at(index);
01385 }

bool BCModel::GetFlagBoundaries (  ) 

Definition at line 330 of file BCModel.cxx.

00331 {
00332    if (!fDataPointLowerBoundaries)
00333       return false;
00334 
00335    if (!fDataPointUpperBoundaries)
00336       return false;
00337 
00338    if (fDataPointLowerBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues())
00339       return false;
00340 
00341    if (fDataPointUpperBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues())
00342       return false;
00343 
00344    return true;
00345 }

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

Definition at line 83 of file BCModel.h.

00084          { return fIndex; };

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

Definition at line 506 of file BCModel.h.

00507          { 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 1164 of file BCModel.cxx.

01165 {
01166    int index1 = par1->GetIndex();
01167    int index2 = par2->GetIndex();
01168 
01169    if (fMCMCFlagsFillHistograms[index1]==false || fMCMCFlagsFillHistograms[index2]==false )
01170    {
01171       BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.",
01172             par1->GetName().data(),par2->GetName().data()));
01173       return 0;
01174    }
01175 
01176    if (index1 == index2)
01177    {
01178       BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
01179       return 0;
01180    }
01181 
01182    BCParameter * npar1 = par1;
01183    BCParameter * npar2 = par2;
01184 
01185    if (index1>index2)
01186    {
01187       npar1 = par2;
01188       npar2 = par1;
01189 
01190       int itmp = index1;
01191       index1 = index2;
01192       index2 = itmp;
01193    }
01194 
01195    // get histogram
01196    TH2D * hist = this -> MCMCGetH2Marginalized(index1, index2);
01197 
01198    if(hist==0)
01199       return 0;
01200 
01201    BCH2D * hprob = new BCH2D();
01202 
01203    // set axis labels
01204    hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data()));
01205    hist->SetXTitle(Form("%s", npar1->GetName().data()));
01206    hist->SetYTitle(Form("%s", npar2->GetName().data()));
01207    hist->SetStats(kFALSE);
01208 
01209    double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
01210    hprob->SetGlobalMode(gmode);
01211 
01212    // set histogram
01213    hprob->SetHistogram(hist);
01214 
01215    return hprob;
01216 }

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

Definition at line 495 of file BCModel.h.

00496          { 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 766 of file BCModel.cxx.

00767 {
00768    if (!parameter)
00769    {
00770       BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist.");
00771       return 0;
00772    }
00773 
00774    int index = parameter -> GetIndex();
00775    if (fMCMCFlagsFillHistograms[index] == false)
00776    {
00777       BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' not filled.",parameter->GetName().data()));
00778       return 0;
00779    }
00780 
00781    // get histogram
00782    TH1D * hist = this -> MCMCGetH1Marginalized(index);
00783    if(!hist)
00784       return 0;
00785 
00786    // set axis labels
00787    hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data()));
00788    hist->SetXTitle(parameter->GetName().data());
00789    hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data()));
00790    hist->SetStats(kFALSE);
00791 
00792    // set histogram
00793    BCH1D * hprob = new BCH1D();
00794    hprob->SetHistogram(hist);
00795 
00796    // set best fit parameter
00797    double bestfit = hprob->GetMode();
00798 
00799    if (fBestFitParametersMarginalized.size() == 0)
00800       for (unsigned int i = 0; i < GetNParameters(); i++)
00801          fBestFitParametersMarginalized.push_back(0.);
00802 
00803    fBestFitParametersMarginalized[index] = bestfit;
00804 
00805    hprob->SetGlobalMode(fBestFitParameters.at(index));
00806 
00807    return hprob;
00808 }

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

Definition at line 93 of file BCModel.h.

00094          { return fModelAPosteriori; };

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

Definition at line 88 of file BCModel.h.

00089          { return fModelAPriori; };

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

Definition at line 78 of file BCModel.h.

00079          { return fName; };

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

Definition at line 100 of file BCModel.cxx.

00101 {
00102    int npoints = 0;
00103    if (fDataSet)
00104       npoints = fDataSet -> GetNDataPoints();
00105    else
00106    {
00107       BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined.");
00108       return ERROR_NOEVENTS;
00109    }
00110 
00111    return npoints;
00112 }

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

Definition at line 149 of file BCModel.h.

00150          { return fNDataPointsMaximum; };

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

Definition at line 144 of file BCModel.h.

00145          { return fNDataPointsMinimum; };

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

Definition at line 98 of file BCModel.h.

00099          { return fNormalization; };

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

Definition at line 154 of file BCModel.h.

00155          { 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 144 of file BCModel.cxx.

00145 {
00146    if (!fParameterSet)
00147       return 0;
00148 
00149    int index = -1;
00150    for (unsigned int i = 0; i < this->GetNParameters(); i++)
00151       if (name == this -> GetParameter(i) -> GetName())
00152          index = i;
00153 
00154    if (index<0)
00155    {
00156       BCLog::OutWarning(Form(
00157                                         "BCModel::GetParameter : Model %s has no parameter named '%s'",
00158                                         (this -> GetName()).data(), name
00159                                         )
00160                                  );
00161       return 0;
00162    }
00163 
00164    return this->GetParameter(index);
00165 }

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

Definition at line 127 of file BCModel.cxx.

00128 {
00129    if (!fParameterSet)
00130       return 0;
00131 
00132    if (index < 0 || index >= (int)this -> GetNParameters())
00133    {
00134       BCLog::OutWarning(
00135             Form("BCModel::GetParameter. Parameter index %d not within range.", index));
00136       return 0;
00137    }
00138 
00139    return fParameterSet -> at(index);
00140 }

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

Definition at line 169 of file BCModel.h.

00170          { return fParameterSet; };

double BCModel::GetPValue (  )  [inline]

Definition at line 533 of file BCModel.h.

00534          { 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 1220 of file BCModel.cxx.

01221 {
01222    double ll = this -> LogLikelihood(par);
01223    int n = this -> GetNDataPoints();
01224 
01225    double sum_sigma=0;
01226    for (int i=0;i<n;i++)
01227       sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index));
01228 
01229    double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma);
01230 
01231    fPValue = TMath::Prob(chi2,n);
01232 
01233    return fPValue;
01234 }

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

Definition at line 1238 of file BCModel.cxx.

01239 {
01240    double ll = this -> LogLikelihood(par);
01241    int n = this -> GetNDataPoints();
01242    int npar = this -> GetNParameters();
01243 
01244    double sum_sigma=0;
01245    for (int i=0;i<n;i++)
01246       sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index));
01247 
01248    double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma);
01249 
01250    fChi2NDoF = chi2/double(n-npar);
01251    fPValueNDoF = TMath::Prob(chi2,n-npar);
01252 
01253    return fPValueNDoF;
01254 }

double BCModel::GetPValueNDoF (  )  [inline]

Definition at line 536 of file BCModel.h.

00537          { 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 1311 of file BCModel.cxx.

01312 {
01313    // check number of entries in vector
01314    if (point.size() != this -> GetNParameters())
01315    {
01316       BCLog::OutWarning("BCModel::HessianMatrixElement. Invalid number of entries in the vector.");
01317       return -1;
01318    }
01319 
01320    // define steps
01321    double nsteps = 1e5;
01322    double dx1 = par1 -> GetRangeWidth() / nsteps;
01323    double dx2 = par2 -> GetRangeWidth() / nsteps;
01324 
01325    // define points at which to evaluate
01326    std::vector<double> xpp = point;
01327    std::vector<double> xpm = point;
01328    std::vector<double> xmp = point;
01329    std::vector<double> xmm = point;
01330 
01331    int idx1 = par1 -> GetIndex();
01332    int idx2 = par2 -> GetIndex();
01333 
01334    xpp[idx1] += dx1;
01335    xpp[idx2] += dx2;
01336 
01337    xpm[idx1] += dx1;
01338    xpm[idx2] -= dx2;
01339 
01340    xmp[idx1] -= dx1;
01341    xmp[idx2] += dx2;
01342 
01343    xmm[idx1] -= dx1;
01344    xmm[idx2] -= dx2;
01345 
01346    // calculate probability at these points
01347    double ppp = this -> Likelihood(xpp);
01348    double ppm = this -> Likelihood(xpm);
01349    double pmp = this -> Likelihood(xmp);
01350    double pmm = this -> Likelihood(xmm);
01351 
01352    // return derivative
01353    return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2);
01354 }

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

Returns the likelihood

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

Definition at line 351 of file BCModel.h.

00352          { 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, BCHistogramFitter, and BCRooInterface.

Definition at line 344 of file BCModel.h.

00345          { return 0.; };

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

00406          { flag_ConditionalProbabilityEntry = false; return 0.; };

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

Overloaded function to evaluate integral.

Reimplemented from BCIntegrate.

Definition at line 520 of file BCModel.cxx.

00521 {
00522    return this -> LogProbabilityNN(parameters);
00523 }

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, BCHistogramFitter, and BCRooInterface.

Definition at line 504 of file BCModel.cxx.

00505 {
00506    double logprob = 0.;
00507 
00508    // add log of conditional probabilities event-by-event
00509    for (unsigned int i=0;i<fDataSet -> GetNDataPoints();i++)
00510    {
00511       BCDataPoint * datapoint = this -> GetDataPoint(i);
00512       logprob += this -> LogConditionalProbabilityEntry(datapoint, parameters);
00513    }
00514 
00515    return logprob;
00516 }

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

00490 {
00491    // check if normalized
00492    if (fNormalization<0. || fNormalization==0.)
00493    {
00494       BCLog::Out(BCLog::warning, BCLog::warning,
00495              "BCModel::LogProbability. Normalization not available or zero.");
00496       return 0.;
00497    }
00498 
00499    return this -> LogProbabilityNN(parameters) - log(fNormalization);
00500 }

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

00477 {
00478    // add log of conditional probability
00479    double logprob = this -> LogLikelihood(parameters);
00480 
00481    // add log of prior probability
00482    logprob += this -> LogAPrioriProbability(parameters);
00483 
00484    return logprob;
00485 }

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

00725 {
00726    BCLog::OutSummary(Form("Running MCMC for model \'%s\'",this->GetName().data()));
00727 
00728    // prepare function fitting
00729    double dx = 0.0;
00730    double dy = 0.0;
00731 
00732    if (fFitFunctionIndexX >= 0)
00733    {
00734       dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX))
00735             / (double)fErrorBandNbinsX;
00736 
00737       dy = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY))
00738             / (double)fErrorBandNbinsY;
00739 
00740       fErrorBandXY = new TH2D(
00741             TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "",
00742             fErrorBandNbinsX,
00743             fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - .5 * dx,
00744             fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + .5 * dx,
00745             fErrorBandNbinsY,
00746             fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - .5 * dy,
00747             fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + .5 * dy);
00748       fErrorBandXY -> SetStats(kFALSE);
00749 
00750       for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
00751          for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
00752             fErrorBandXY -> SetBinContent(ix, iy, 0.0);
00753    }
00754 
00755    this -> MCMCMetropolis();
00756    this -> FindModeMCMC();
00757 
00758 // PrintResults(Form("%s.txt", this -> GetName().data()));
00759 
00760    return 1;
00761 }

double BCModel::Normalize (  ) 

Integrates over the un-normalized probability and updates fNormalization.

Definition at line 544 of file BCModel.cxx.

00545 {
00546    BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",this->GetName().data()));
00547 
00548    unsigned int n = this -> GetNvar();
00549 
00550    // initialize BCIntegrate if not done already
00551    if (n == 0)
00552    {
00553       this->SetParameters(fParameterSet);
00554       n = this->GetNvar();
00555    }
00556 
00557    // integrate and get best fit parameters
00558    // maybe we have to remove the mode finding from here in the future
00559    fNormalization = this -> Integrate();
00560 
00561    BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization));
00562 
00563    return fNormalization;
00564 }

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

Definition at line 966 of file BCModel.cxx.

00967 {
00968    if(!fMCMCFlagFillHistograms)
00969    {
00970       BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled.");
00971       return 0;
00972    }
00973 
00974    int npar = GetNParameters();
00975 
00976    if(fMCMCH1Marginalized.size()==0 || (fMCMCH2Marginalized.size()==0 && npar>1))
00977    {
00978       BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00979       return 0;
00980    }
00981 
00982    // if there's only one parameter, we just want to call Print()
00983    if (fMCMCH1Marginalized.size()==1 && fMCMCH2Marginalized.size()==0)
00984    {
00985       BCParameter * a = GetParameter(0);
00986       GetMarginalized(a)->Print(file);
00987       return 1;
00988    }
00989 
00990    int c_width=600; // default canvas width
00991    int c_height=600; // default canvas height
00992 
00993    int type = 112; // landscape
00994 
00995    if (hdiv > vdiv)
00996    {
00997       if(hdiv>3)
00998       {
00999          c_width=1000;
01000          c_height=700;
01001       }
01002       else
01003       {
01004          c_width=800;
01005          c_height=600;
01006       }
01007    }
01008    else if(hdiv < vdiv)
01009    {
01010       if(hdiv>3)
01011       {
01012          c_height=1000;
01013          c_width=700;
01014       }
01015       else
01016       {
01017          c_height=800;
01018          c_width=600;
01019       }
01020       type=111;
01021    }
01022 
01023    // calculate number of plots
01024    int nplots2d = npar * (npar-1)/2;
01025    int nplots = npar + nplots2d;
01026 
01027    // give out warning if too many plots
01028    BCLog::OutSummary(Form(
01029          "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s",
01030          npar,nplots2d,nplots,file));
01031    if(nplots>100)
01032       BCLog::OutDetail("This can take a while...");
01033 
01034    // setup the canvas and postscript file
01035    TCanvas * c = new TCanvas( "c","canvas",c_width,c_height);
01036 
01037    TPostScript * ps = new TPostScript(file,type);
01038 
01039    if(type==112)
01040       ps->Range(24,16);
01041    else
01042       ps->Range(16,24);
01043 
01044    // draw all 1D distributions
01045    ps->NewPage();
01046    c->cd();
01047    c->Clear();
01048    c->Divide(hdiv,vdiv);
01049 
01050    int n=0;
01051    for(int i=0;i<npar;i++)
01052    {
01053       // get corresponding parameter
01054       BCParameter * a = GetParameter(i);
01055 
01056       // check if histogram exists
01057       if ( !GetMarginalized(a) )
01058          continue;
01059 
01060       // check if histogram is filled
01061       if ( GetMarginalized(a)->GetHistogram()->Integral() <= 0)
01062          continue;
01063 
01064       // if current page is full, switch to new page
01065       if(i!=0 && i%(hdiv*vdiv)==0)
01066       {
01067          c->Update();
01068          ps->NewPage();
01069          c->cd();
01070          c->Clear();
01071          c->Divide(hdiv,vdiv);
01072       }
01073 
01074       // go to next pad
01075       c->cd(i%(hdiv*vdiv)+1);
01076 
01077       this -> GetMarginalized(a) -> Draw();
01078       n++;
01079 
01080       if(n%100==0)
01081          BCLog::OutDetail(Form(" --> %d plots done",n));
01082    }
01083 
01084    if (n > 0)
01085    {
01086       c->Update();
01087 
01088       // draw all the 2D distributions
01089       ps->NewPage();
01090       c->cd();
01091       c->Clear();
01092    }
01093 
01094    c->Divide(hdiv,vdiv);
01095 
01096    int k=0;
01097    for(int i=0;i<npar-1;i++)
01098    {
01099       if(!fMCMCFlagsFillHistograms[i])
01100          continue;
01101       for(int j=i+1;j<npar;j++)
01102       {
01103          if(!fMCMCFlagsFillHistograms[j])
01104             continue;
01105 
01106          // get corresponding parameters
01107          BCParameter * a = GetParameter(i);
01108          BCParameter * b = GetParameter(j);
01109 
01110          // check if histogram exists
01111          if ( !GetMarginalized(a,b) )
01112             continue;
01113 
01114          // check if histogram is filled
01115          if ( GetMarginalized(a,b)->GetHistogram()->Integral() <= 0 )
01116             continue;
01117 
01118          // if current page is full, switch to new page
01119          if(k!=0 && k%(hdiv*vdiv)==0)
01120          {
01121             c->Update();
01122             ps->NewPage();
01123             c->cd();
01124             c->Clear();
01125             c->Divide(hdiv,vdiv);
01126          }
01127 
01128          // go to next pad
01129          c->cd(k%(hdiv*vdiv)+1);
01130 
01131          double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.;
01132          double deltaa = (a->GetUpperLimit() - a->GetLowerLimit());
01133          if (deltaa <= 1e-7 * meana)
01134             continue;
01135 
01136          double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.;
01137          double deltab = (b->GetUpperLimit() - b->GetLowerLimit());
01138          if (deltab <= 1e-7 * meanb)
01139             continue;
01140 
01141          GetMarginalized(a,b)->Draw(52);
01142          k++;
01143 
01144          if((n+k)%100==0)
01145             BCLog::OutDetail(Form(" --> %d plots done",n+k));
01146       }
01147    }
01148 
01149    if( (n+k)>100 && (n+k)%100 != 0 )
01150       BCLog::OutDetail(Form(" --> %d plots done",n+k));
01151 
01152    c->Update();
01153    ps->Close();
01154 
01155    delete c;
01156    delete ps;
01157 
01158    // return total number of drawn histograms
01159    return n+k;
01160 }

int BCModel::PrintAllMarginalized1D ( const char *  filebase  ) 

Definition at line 907 of file BCModel.cxx.

00908 {
00909    if(fMCMCH1Marginalized.size()==0)
00910    {
00911       BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00912       return 0;
00913    }
00914 
00915    int n=this->GetNParameters();
00916    for(int i=0;i<n;i++)
00917    {
00918       BCParameter * a = this->GetParameter(i);
00919       if (this -> GetMarginalized(a))
00920          this -> GetMarginalized(a) -> Print(Form("%s_1D_%s.ps",filebase,a->GetName().data()));
00921    }
00922 
00923    return n;
00924 }

int BCModel::PrintAllMarginalized2D ( const char *  filebase  ) 

Definition at line 928 of file BCModel.cxx.

00929 {
00930    if(fMCMCH2Marginalized.size()==0)
00931    {
00932       BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00933       return 0;
00934    }
00935 
00936    int k=0;
00937    int n=this->GetNParameters();
00938    for(int i=0;i<n-1;i++)
00939    {
00940       for(int j=i+1;j<n;j++)
00941       {
00942          BCParameter * a = this->GetParameter(i);
00943          BCParameter * b = this->GetParameter(j);
00944 
00945          double meana = (a -> GetLowerLimit() + a -> GetUpperLimit()) / 2.0;
00946          double deltaa = (a -> GetUpperLimit() - a -> GetLowerLimit());
00947          if (deltaa <= 1e-7 * meana)
00948             continue;
00949 
00950          double meanb = (b -> GetLowerLimit() + b -> GetUpperLimit()) / 2.0;
00951          double deltab = (b -> GetUpperLimit() - b -> GetLowerLimit());
00952          if (deltab <= 1e-7 * meanb)
00953             continue;
00954 
00955          if (this -> GetMarginalized(a,b))
00956          this -> GetMarginalized(a,b) -> Print(Form("%s_2D_%s_%s.ps",filebase,a->GetName().data(),b->GetName().data()));
00957          k++;
00958       }
00959    }
00960 
00961    return k;
00962 }

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

01676 {
01677    // check number of entries in vector
01678    if (parameters.size() != GetNParameters())
01679    {
01680       BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector");
01681       return;
01682    }
01683 
01684    // print to screen
01685    std::cout
01686       << std::endl
01687       << " Hessian matrix elements: " << std::endl
01688       << " Point: ";
01689 
01690    for (int i = 0; i < int(parameters.size()); i++)
01691       std::cout << parameters.at(i) << " ";
01692    std::cout << std::endl;
01693 
01694    // loop over all parameter pairs
01695    for (unsigned int i = 0; i < GetNParameters(); i++)
01696       for (unsigned int j = 0; j < i; j++)
01697       {
01698          // calculate Hessian matrix element
01699          double hessianmatrixelement = HessianMatrixElement(fParameterSet->at(i),
01700                fParameterSet->at(j), parameters);
01701 
01702          // print to screen
01703          std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl;
01704       }
01705 }

void BCModel::PrintResults ( const char *  file  ) 

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

Definition at line 1452 of file BCModel.cxx.

01453 {
01454    // print summary of Markov Chain Monte Carlo
01455 
01456    // open file
01457    ofstream ofi(file);
01458 
01459    // check if file is open
01460    if(!ofi.is_open())
01461    {
01462       std::cerr << "Couldn't open file " << file <<std::endl;
01463       return;
01464    }
01465 
01466    // number of parameters and chains
01467    int npar = MCMCGetNParameters();
01468    int nchains = MCMCGetNChains();
01469 
01470    // check convergence
01471    bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0;
01472 
01473    ofi
01474       << std::endl
01475       << " -----------------------------------------------------" << std::endl
01476       << " Summary" << std::endl
01477       << " -----------------------------------------------------" << std::endl
01478       << std::endl;
01479 
01480    ofi
01481       << " Model summary" << std::endl
01482       << " =============" << std::endl
01483       << " Model: " << fName.data() << std::endl
01484       << " Number of parameters: " << npar << std::endl
01485       << " List of Parameters and ranges:" << std::endl;
01486    for (int i = 0; i < npar; ++i)
01487       ofi
01488          << "  (" << i << ") Parameter \""
01489          << fParameterSet -> at(i) -> GetName().data() << "\""
01490          << ": " << fParameterSet -> at(i) -> GetLowerLimit()
01491          << " - "
01492          << fParameterSet -> at(i) -> GetUpperLimit() << std::endl;
01493    ofi << std::endl;
01494 
01495    // give warning if MCMC did not converge
01496    if (!flag_conv && fMCMCFlagRun)
01497       ofi
01498          << " WARNING: the Markov Chain did not converge!" << std::endl
01499          << " Be cautious using the following results!" << std::endl
01500          << std::endl;
01501 
01502    // print results of marginalization (if MCMC was run)
01503    if ( fMCMCFlagRun && fMCMCFlagFillHistograms)
01504    {
01505       ofi
01506          << " Results of the marginalization" << std::endl
01507          << " ==============================" << std::endl
01508          << " List of parameters and properties of the marginalized" << std::endl
01509          << " distributions:" << std::endl;
01510       for (int i = 0; i < npar; ++i)
01511       {
01512          if (!fMCMCFlagsFillHistograms.at(i))
01513             continue;
01514 
01515          BCH1D * bch1d = GetMarginalized(fParameterSet -> at(i));
01516 
01517          ofi
01518             << "  (" << i << ") Parameter \""
01519             << fParameterSet->at(i)->GetName().data() << "\"" << std::endl
01520             << "      Mean +- RMS:         "
01521             << std::setprecision(4) << bch1d->GetMean()
01522             << " +- "
01523             << std::setprecision(4) << bch1d->GetRMS() << std::endl
01524             << "      Median +- sigma:     "
01525             << std::setprecision(4) << bch1d->GetMedian()
01526             << " +  " << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian()
01527             << " - " << std::setprecision(4) << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl
01528             << "      (Marginalized) mode: " << bch1d->GetMode() << std::endl
01529             << "      Smallest interval(s) containing 68% and local modes:" << std::endl;
01530 
01531          std::vector <double> v;
01532          v = bch1d->GetSmallestIntervals(0.68);
01533          int ninter = int(v.size());
01534 
01535          for (int j = 0; j < ninter; j+=5)
01536             ofi
01537                << "       " << v.at(j) << " - " << v.at(j+1)
01538                << " (local mode at " << v.at(j+3)
01539                << " with rel. height " << v.at(j+2)
01540                << "; rel. area " << v.at(j+4) << ")"
01541                << std::endl;
01542 
01543          ofi
01544             << "       5% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.05) << std::endl
01545             << "      10% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.10) << std::endl
01546             << "      16% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.16) << std::endl
01547             << "      84% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.85) << std::endl
01548             << "      90% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.90) << std::endl
01549             << "      95% quantile:        " << std::setprecision(4) << bch1d->GetQuantile(0.95) << std::endl
01550             << std::endl;
01551       }
01552       ofi << std::endl;
01553    }
01554 
01555    ofi
01556       << " Results of the optimization" << std::endl
01557       << " ===========================" << std::endl
01558       << " Optimization algorithm used: ";
01559    switch(GetOptimizationMethodMode())
01560    {
01561       case BCIntegrate::kOptSA:
01562          ofi << " Simulated Annealing" << std::endl;
01563          break;
01564       case BCIntegrate::kOptMinuit:
01565          ofi << " Minuit" << std::endl;
01566          break;
01567       case BCIntegrate::kOptMetropolis:
01568          ofi << " MCMC" << std::endl;
01569          break;
01570    }
01571 
01572    if (int(fBestFitParameters.size()) > 0)
01573    {
01574       ofi << " List of parameters and global mode:" << std::endl;
01575       for (int i = 0; i < npar; ++i)
01576          ofi
01577             << "  (" << i << ") Parameter \""
01578             << fParameterSet->at(i)->GetName().data() << "\": "
01579             << fBestFitParameters.at(i) << " +- " << fBestFitParameterErrors.at(i) << std::endl;
01580       ofi << std::endl;
01581    }
01582    else
01583    {
01584       ofi << " No best fit information available." << std::endl;
01585       ofi << std::endl;
01586    }
01587 
01588    if (fPValue >= 0.)
01589    {
01590       ofi
01591          << " Results of the model test" << std::endl
01592          << " =========================" << std::endl
01593          << " p-value at global mode: " << fPValue << std::endl
01594          << std::endl;
01595    }
01596 
01597    if ( fMCMCFlagRun )
01598    {
01599       ofi
01600          << " Status of the MCMC" << std::endl
01601          << " ==================" << std::endl
01602          << " Convergence reached: " << (flag_conv ? "yes" : "no") << std::endl;
01603 
01604       if (flag_conv)
01605          ofi << " Number of iterations until convergence: " << MCMCGetNIterationsConvergenceGlobal() << std::endl;
01606       else
01607          ofi
01608             << " WARNING: the Markov Chain did not converge! Be\n"
01609             << " cautious using the following results!" << std::endl
01610             << std::endl;
01611       ofi
01612          << " Number of chains:                       " << MCMCGetNChains() << std::endl
01613          << " Number of iterations per chain:         " << MCMCGetNIterationsRun() << std::endl
01614          << " Average efficiencies:" << std::endl;
01615 
01616       std::vector <double> efficiencies;
01617       efficiencies.assign(npar, 0.);
01618 
01619       for (int ipar = 0; ipar < npar; ++ipar)
01620          for (int ichain = 0; ichain < nchains; ++ichain)
01621          {
01622             int index = ichain * npar + ipar;
01623             efficiencies[ipar] += double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index) + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.;
01624          }
01625 
01626       for (int ipar = 0; ipar < npar; ++ipar)
01627          ofi
01628             << "  (" << ipar << ") Parameter \""
01629             << fParameterSet->at(ipar)->GetName().data() << "\": "
01630             << efficiencies.at(ipar) << "%" << std::endl;
01631       ofi << std::endl;
01632    }
01633 
01634    ofi
01635       << " -----------------------------------------------------" << std::endl
01636       << std::endl
01637       << " Notes" << std::endl
01638       << " =====" << std::endl
01639       << " (i) Median +- sigma denotes the median, m, of the" << std::endl
01640       << "     marginalized distribution and the intervals from" << std::endl
01641       << "     m to the 16% and 84% quantiles." << std::endl
01642       << " -----------------------------------------------------" << std::endl;
01643    
01644    // close file
01645 // ofi.close;
01646 }

void BCModel::PrintShortFitSummary ( int  chi2flag = 0  ) 

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

Definition at line 1650 of file BCModel.cxx.

01651 {
01652    BCLog::OutSummary("---------------------------------------------------");
01653    BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data()));
01654    BCLog::OutSummary(Form("   Number of parameters:  Npar  = %i", GetNParameters()));
01655    BCLog::OutSummary(Form("   Number of data points: Ndata = %i", GetNDataPoints()));
01656    BCLog::OutSummary("   Number of degrees of freedom:");
01657    BCLog::OutSummary(Form("      NDoF = Ndata - Npar = %i", GetNDataPoints()-GetNParameters()));
01658 
01659    BCLog::OutSummary("   Best fit parameters (global):");
01660    for (unsigned int i = 0; i < GetNParameters(); ++i)
01661       BCLog::OutSummary(Form("      %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i)));
01662 
01663    BCLog::OutSummary("   Goodness-of-fit test:");
01664    BCLog::OutSummary(Form("      p-value = %.3g", GetPValue()));
01665    if(chi2flag)
01666    {
01667       BCLog::OutSummary(Form("      p-value corrected for NDoF = %.3g", GetPValueNDoF()));
01668       BCLog::OutSummary(Form("      chi2 / NDoF = %.3g", GetChi2NDoF()));
01669    }
01670    BCLog::OutSummary("---------------------------------------------------");
01671 }

void BCModel::PrintSummary (  ) 

Prints a summary on the screen.

Definition at line 1389 of file BCModel.cxx.

01390 {
01391    int nparameters = this -> GetNParameters();
01392 
01393    // model summary
01394    std::cout
01395       << std::endl
01396       << "   ---------------------------------" << std::endl
01397       << "    Model : " << fName.data() << std::endl
01398       << "   ---------------------------------"<< std::endl
01399       << "     Index                : " << fIndex << std::endl
01400       << "     Number of parameters : " << nparameters << std::endl
01401       << std::endl
01402       << "     - Parameters : " << std::endl
01403       << std::endl;
01404 
01405    // parameter summary
01406    for (int i=0; i<nparameters; i++)
01407       fParameterSet -> at(i) -> PrintSummary();
01408 
01409    // best fit parameters
01410    if (this -> GetBestFitParameters().size() > 0)
01411    {
01412       std::cout
01413          << std::endl
01414          << "     - Best fit parameters :" << std::endl
01415          << std::endl;
01416 
01417       for (int i=0; i<nparameters; i++)
01418       {
01419          std::cout
01420             << "       " << fParameterSet -> at(i) -> GetName().data()
01421             << " = " << this -> GetBestFitParameter(i)
01422             << " (overall)" << std::endl;
01423          if ((int)fBestFitParametersMarginalized.size() == nparameters)
01424             std::cout
01425                << "       " << fParameterSet -> at(i) -> GetName().data()
01426                << " = " << this -> GetBestFitParameterMarginalized(i)
01427                << " (marginalized)" << std::endl;
01428       }
01429    }
01430 
01431    std::cout << std::endl;
01432 
01433    // model testing
01434    if (fPValue >= 0)
01435    {
01436       double likelihood = this -> Likelihood(this -> GetBestFitParameters());
01437       std::cout
01438          << "   - Model testing:" << std::endl
01439          << std::endl
01440          << "       p(data|lambda*) = " << likelihood << std::endl
01441          << "       p-value         = " << fPValue << std::endl
01442          << std::endl;
01443    }
01444 
01445    // normalization
01446    if (fNormalization > 0)
01447       std::cout << "     Normalization        : " << fNormalization << std::endl;
01448 }

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

00380          { 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 365 of file BCModel.h.

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

int BCModel::ReadErrorBandFromFile ( const char *  file  ) 

Read

Definition at line 876 of file BCModel.cxx.

00877 {
00878    TFile * froot = new TFile(file);
00879    if(!froot->IsOpen())
00880    {
00881       BCLog::OutWarning(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.",file));
00882       return 0;
00883    }
00884 
00885    int r=0;
00886 
00887    TH2D * h2 = (TH2D*)froot->Get("errorbandxy");
00888    if(h2)
00889    {
00890       h2->SetDirectory(0);
00891       h2->SetName(TString::Format("errorbandxy_%d",BCLog::GetHIndex()));
00892       this->SetErrorBandHisto(h2);
00893       r=1;
00894    }
00895    else
00896       BCLog::OutWarning(Form(
00897             "BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",
00898             file));
00899 
00900    froot->Close();
00901 
00902    return r;
00903 }

int BCModel::ReadMarginalizedFromFile ( const char *  file  ) 

Read

Definition at line 812 of file BCModel.cxx.

00813 {
00814    TFile * froot = new TFile(file);
00815    if(!froot->IsOpen())
00816    {
00817       BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.",file));
00818       return 0;
00819    }
00820 
00821    // We reset the MCMCEngine here for the moment.
00822    // In the future maybe we only want to do this if the engine
00823    // wans't initialized at all or when there were some changes
00824    // in the model.
00825    // But maybe we want reset everything since we're overwriting
00826    // the marginalized distributions anyway.
00827    this -> MCMCInitialize();
00828 
00829    int k=0;
00830    int n=this->GetNParameters();
00831    for(int i=0;i<n;i++)
00832    {
00833       BCParameter * a = this -> GetParameter(i);
00834       TKey * key = froot -> GetKey(TString::Format("hist_%s_%s", this -> GetName().data(), a -> GetName().data()));
00835       if(key)
00836       {
00837          TH1D * h1 = (TH1D*) key -> ReadObjectAny(TH1D::Class());
00838          h1->SetDirectory(0);
00839          if(this->SetMarginalized(i,h1))
00840             k++;
00841       }
00842       else
00843          BCLog::OutWarning(Form(
00844                "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.",
00845                this -> GetName().data(), a -> GetName().data(), file));
00846    }
00847 
00848    for(int i=0;i<n-1;i++)
00849    {
00850       for(int j=i+1;j<n;j++)
00851       {
00852          BCParameter * a = this -> GetParameter(i);
00853          BCParameter * b = this -> GetParameter(j);
00854          TKey * key = froot -> GetKey(TString::Format("hist_%s_%s_%s", this -> GetName().data(), a -> GetName().data(), b -> GetName().data()));
00855          if(key)
00856          {
00857             TH2D * h2 = (TH2D*) key -> ReadObjectAny(TH2D::Class());
00858             h2->SetDirectory(0);
00859             if(this->SetMarginalized(i,j,h2))
00860                k++;
00861          }
00862          else
00863             BCLog::OutWarning(Form(
00864                   "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.",
00865                   this -> GetName().data(), a -> GetName().data(), b -> GetName().data(), file));
00866       }
00867    }
00868 
00869    froot->Close();
00870 
00871    return k;
00872 }

int BCModel::ReadMode ( const char *  file  ) 

Read mode from file created by WriteMode() call

Definition at line 684 of file BCModel.cxx.

00685 {
00686    ifstream ifi(file);
00687    if(!ifi.is_open())
00688    {
00689       BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.",file));
00690       return 0;
00691    }
00692 
00693    int npar = fParameterSet -> size();
00694    std::vector <double> mode;
00695 
00696    int i=0;
00697    while (i<npar && !ifi.eof())
00698    {
00699       double a;
00700       ifi>>a;
00701       mode.push_back(a);
00702       i++;
00703    }
00704 
00705    if(i<npar)
00706    {
00707       BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.",file));
00708       BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.",npar,i));
00709       return 0;
00710    }
00711 
00712    BCLog::OutSummary(Form("#  Read in best fit parameters (mode) for model \'%s\' from file %s:",fName.data(),file));
00713    this->SetMode(mode);
00714    for(int j=0 ; j<npar; j++)
00715       BCLog::OutSummary(Form("#    -> Parameter %d : %s = %e", j, fParameterSet->at(j)->GetName().data(), fBestFitParameters[j]));
00716 
00717    BCLog::OutWarning("#  ! Best fit values obtained before this call will be overwritten !");
00718 
00719    return npar;
00720 }

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

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

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

Definition at line 374 of file BCModel.cxx.

00375 {
00376    // check if data set exists
00377    if (!fDataSet)
00378    {
00379       BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first.");
00380       return;
00381    }
00382 
00383    // check if index is within range
00384    if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues())
00385    {
00386       BCLog::OutError("BCModel::SetDataBoundaries : Index out of range.");
00387       return;
00388    }
00389 
00390    // check if boundary data points exist
00391    if (!fDataPointLowerBoundaries)
00392       fDataPointLowerBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues());
00393 
00394    if (!fDataPointUpperBoundaries)
00395       fDataPointUpperBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues());
00396 
00397    if (fDataFixedValues.size() == 0)
00398       fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false);
00399 
00400    // set boundaries
00401    fDataPointLowerBoundaries -> SetValue(index, lowerboundary);
00402    fDataPointUpperBoundaries -> SetValue(index, upperboundary);
00403    fDataFixedValues[index] = fixed;
00404 }

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

Sets the data set.

Parameters:
dataset A data set

Definition at line 279 of file BCModel.h.

00280          { fDataSet = dataset; fNormalization = -1.0; };

void BCModel::SetErrorBandContinuous ( bool  flag  ) 

Sets the error band flag to continuous function

Definition at line 408 of file BCModel.cxx.

00409 {
00410    fErrorBandContinuous = flag;
00411 
00412    if (flag)
00413       return;
00414 
00415    // clear x-values
00416    fErrorBandX.clear();
00417 
00418    // copy data x-values
00419    for (unsigned int i = 0; i < fDataSet -> GetNDataPoints(); ++i)
00420       fErrorBandX.push_back(fDataSet -> GetDataPoint(i) -> GetValue(fFitFunctionIndexX));
00421 }

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

Definition at line 557 of file BCModel.h.

00558          { fGoFNChains=n; };

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

Definition at line 545 of file BCModel.h.

00546          { fGoFNIterationsMax=n; };

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

Definition at line 551 of file BCModel.h.

00552          { 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 246 of file BCModel.h.

00247          { 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 266 of file BCModel.h.

00267          { 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 259 of file BCModel.h.

00260          { fModelAPriori = probability; };

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

Sets the name of the model.

Parameters:
name Name of the model

Definition at line 240 of file BCModel.h.

00241          { 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 169 of file BCModel.cxx.

00170 {
00171    BCParameter * p = this -> GetParameter(parname);
00172    if(!p)
00173    {
00174       BCLog::OutWarning(Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname));
00175       return;
00176    }
00177 
00178    this -> BCIntegrate::SetNbins(nbins, p -> GetIndex());
00179 }

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

Sets the maximum number of data points.

Definition at line 296 of file BCModel.h.

00297          { fNDataPointsMaximum = maximum; };

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

Sets the minimum number of data points.

Definition at line 291 of file BCModel.h.

00292          { 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 273 of file BCModel.h.

00274          { fNormalization = norm; };

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

Set all parameters of the model using a BCParameterSet container.

pointer to parameter set

Definition at line 252 of file BCModel.h.

00253          { fParameterSet = parset; };

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

Definition at line 364 of file BCModel.cxx.

00365 {
00366    if (index < 0 || index > dataset -> GetNDataPoints())
00367       return;
00368 
00369    this -> SetSingleDataPoint(dataset -> GetDataPoint(index));
00370 }

void BCModel::SetSingleDataPoint ( BCDataPoint *  datapoint  ) 

Sets a single data point as data set.

Parameters:
datapoint A data point

Definition at line 349 of file BCModel.cxx.

00350 {
00351    // create new data set consisting of a single data point
00352    BCDataSet * dataset = new BCDataSet();
00353 
00354    // add the data point
00355    dataset -> AddDataPoint(datapoint);
00356 
00357    // set this new data set
00358    this -> SetDataSet(dataset);
00359 
00360 }

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

Converts a vector of doubles into a BCDataPoint

Definition at line 1709 of file BCModel.cxx.

01710 {
01711    int sizeofvector = int(data.size());
01712    BCDataPoint * datapoint = new BCDataPoint(sizeofvector);
01713    datapoint -> SetValues(data);
01714    return datapoint;
01715 }

void BCModel::WriteMode ( const char *  file  ) 

Write mode into file

Definition at line 643 of file BCModel.cxx.

00644 {
00645    ofstream ofi(file);
00646    if(!ofi.is_open())
00647    {
00648       std::cerr<<"Couldn't open file "<<file<<std::endl;
00649       return;
00650    }
00651 
00652    int npar = fParameterSet -> size();
00653    for (int i=0; i<npar; i++)
00654       ofi<<fBestFitParameters.at(i)<<std::endl;
00655 
00656    ofi<<std::endl;
00657    ofi<<"#######################################################################"<<std::endl;
00658    ofi<<"#"<<std::endl;
00659    ofi<<"#  This file was created automatically by BCModel::WriteMode() call."<<std::endl;
00660    ofi<<"#  It can be read in by call to BCModel::ReadMode()."<<std::endl;
00661    ofi<<"#  Do not modify it unless you know what you're doing."<<std::endl;
00662    ofi<<"#"<<std::endl;
00663    ofi<<"#######################################################################"<<std::endl;
00664    ofi<<"#"<<std::endl;
00665    ofi<<"#  Best fit parameters (mode) for model:"<<std::endl;
00666    ofi<<"#  \'"<<fName.data()<<"\'"<<std::endl;
00667    ofi<<"#"<<std::endl;
00668    ofi<<"#  Number of parameters: "<<npar<<std::endl;
00669    ofi<<"#  Parameters ordered as above:"<<std::endl;
00670 
00671    for (int i=0; i<npar; i++)
00672    {
00673       ofi<<"#     "<<i<<": ";
00674       ofi<<fParameterSet->at(i)->GetName().data()<<" = ";
00675       ofi<<fBestFitParameters.at(i)<<std::endl;
00676    }
00677 
00678    ofi<<"#"<<std::endl;
00679    ofi<<"########################################################################"<<std::endl;
00680 }

Member Data Documentation

double BCModel::fChi2NDoF [protected]

Definition at line 630 of file BCModel.h.

BCDataSet* BCModel::fDataSet [protected]

A data set

Definition at line 612 of file BCModel.h.

int BCModel::fGoFNChains [protected]

Definition at line 646 of file BCModel.h.

int BCModel::fGoFNIterationsMax [protected]

Definition at line 636 of file BCModel.h.

int BCModel::fGoFNIterationsRun [protected]

Definition at line 641 of file BCModel.h.

int BCModel::fIndex [protected]

Index of the model.

Definition at line 592 of file BCModel.h.

bool BCModel::flag_ConditionalProbabilityEntry [protected]

A flag for overloading ConditionalProbabilityEntry

Definition at line 624 of file BCModel.h.

double BCModel::fModelAPosteriori [protected]

The model a posteriori probability.

Definition at line 604 of file BCModel.h.

double BCModel::fModelAPriori [protected]

The model prior probability.

Definition at line 600 of file BCModel.h.

std::string BCModel::fName [protected]

Name of the model.

Definition at line 596 of file BCModel.h.

unsigned int BCModel::fNDataPointsMaximum [protected]

Maximum number of data points

Definition at line 620 of file BCModel.h.

unsigned int BCModel::fNDataPointsMinimum [protected]

Minimum number of data points

Definition at line 616 of file BCModel.h.

double BCModel::fNormalization [private]

The Likelihood normalization.

Definition at line 660 of file BCModel.h.

BCParameterSet* BCModel::fParameterSet [protected]

A model parameter container.

Definition at line 608 of file BCModel.h.

double BCModel::fPValue [protected]

The p-value

Definition at line 628 of file BCModel.h.

double BCModel::fPValueNDoF [protected]

Definition at line 631 of file BCModel.h.

The documentation for this class was generated from the following files:
Generated on Tue Oct 6 09:48:22 2009 for Bayesian Analysis Toolkit by  doxygen 1.6.1