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

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

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

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

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

01303 {
01304    // ...
01305 }

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

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

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

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

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

01163 {
01164    if (fMCMCFlagFillHistograms == false)
01165       {
01166          BCLog::Out(BCLog::warning, BCLog::warning,
01167                          "BCModel::GetMarginalized. Histogram filling was switched off. Marginalized distribuions are not generated."); 
01168       }
01169 
01170 // if(fMCMCH2Marginalized.size()==0)
01171 // {
01172 //    BCLog::Out(BCLog::warning, BCLog::warning,
01173 //           "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this.");
01174 //    return 0;
01175 // }
01176 
01177    int index1 = parameter1 -> GetIndex();
01178    int index2 = parameter2 -> GetIndex();
01179 
01180    if (index1 == index2)
01181    {
01182       BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available.");
01183       return 0;
01184    }
01185 
01186    if (index1 > index2)
01187    {
01188       int itmp = index1;
01189       index1 = index2;
01190       index2 = itmp;
01191    }
01192 
01193    // get histogram
01194    TH2D * hist = this -> MCMCGetH2Marginalized(index1, index2);
01195 
01196    if(hist==0)
01197       return 0;
01198 
01199    BCH2D * hprob = new BCH2D();
01200 
01201    // set axis labels
01202    hist -> SetName(Form("hist_%s_%s_%s", this -> GetName().data(), parameter1 -> GetName().data(), parameter2 -> GetName().data()));
01203    hist -> SetXTitle(Form("%s", parameter1 -> GetName().data()));
01204    hist -> SetYTitle(Form("%s", parameter2 -> GetName().data()));
01205    hist -> SetStats(kFALSE);
01206 
01207    double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) };
01208    hprob->SetGlobalMode(gmode);
01209 
01210    // set histogram
01211    hprob -> SetHistogram(hist);
01212 
01213    return hprob;
01214 }

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::Out(BCLog::warning, BCLog::warning,
00771                          "BCModel::GetMarginalized(). Parameter does not exist."); 
00772          return 0; 
00773       }
00774 
00775    if (fMCMCFlagFillHistograms == false)
00776       {
00777          BCLog::Out(BCLog::warning, BCLog::warning,
00778                          "BCModel::GetMarginalized. Histogram filling was switched off. Marginalized distribuions are not generated."); 
00779       }
00780 
00781 // if(fMCMCH1Marginalized.size()==0)
00782 // {
00783 //    BCLog::Out(BCLog::warning, BCLog::warning,
00784 //          "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this.");
00785 //    return 0;
00786 // }
00787 
00788    int index = parameter -> GetIndex();
00789 
00790    // get histogram
00791    TH1D * hist = this -> MCMCGetH1Marginalized(index);
00792    if(!hist)
00793       return 0;
00794 
00795    BCH1D * hprob = new BCH1D();
00796 
00797    // set axis labels
00798    hist -> SetName(Form("hist_%s_%s", this -> GetName().data(), parameter -> GetName().data()));
00799    hist -> SetXTitle(parameter -> GetName().data());
00800    hist -> SetYTitle(Form("p(%s|data)", parameter -> GetName().data()));
00801    hist -> SetStats(kFALSE);
00802 
00803    // set histogram
00804    hprob -> SetHistogram(hist);
00805 
00806    // set best fit parameter
00807    double bestfit = hprob -> GetMode();
00808 
00809    if (fBestFitParametersMarginalized.size() == 0)
00810       for (unsigned int i = 0; i < this -> GetNParameters(); i++)
00811          fBestFitParametersMarginalized.push_back(0.0);
00812 
00813    fBestFitParametersMarginalized[index] = bestfit;
00814 
00815    hprob->SetGlobalMode(fBestFitParameters.at(index));
00816 
00817    return hprob;
00818 }

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

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

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

Definition at line 1236 of file BCModel.cxx.

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

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

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

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

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

int BCModel::PrintAllMarginalized1D ( const char *  filebase  ) 

Definition at line 917 of file BCModel.cxx.

00918 {
00919    if(fMCMCH1Marginalized.size()==0)
00920    {
00921       BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available.");
00922       return 0;
00923    }
00924 
00925    int n=this->GetNParameters();
00926    for(int i=0;i<n;i++)
00927    {
00928       BCParameter * a = this->GetParameter(i);
00929       if (this -> GetMarginalized(a))
00930          this -> GetMarginalized(a) -> Print(Form("%s_1D_%s.ps",filebase,a->GetName().data()));
00931    }
00932 
00933    return n;
00934 }

int BCModel::PrintAllMarginalized2D ( const char *  filebase  ) 

Definition at line 938 of file BCModel.cxx.

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

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

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

void BCModel::PrintResults ( const char *  file  ) 

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

Definition at line 1450 of file BCModel.cxx.

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

void BCModel::PrintShortFitSummary ( int  chi2flag = 0  ) 

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

Definition at line 1648 of file BCModel.cxx.

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

void BCModel::PrintSummary (  ) 

Prints a summary on the screen.

Definition at line 1387 of file BCModel.cxx.

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

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

00887 {
00888    TFile * froot = new TFile(file);
00889    if(!froot->IsOpen())
00890    {
00891       BCLog::OutWarning(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.",file));
00892       return 0;
00893    }
00894 
00895    int r=0;
00896 
00897    TH2D * h2 = (TH2D*)froot->Get("errorbandxy");
00898    if(h2)
00899    {
00900       h2->SetDirectory(0);
00901       h2->SetName(TString::Format("errorbandxy_%d",BCLog::GetHIndex()));
00902       this->SetErrorBandHisto(h2);
00903       r=1;
00904    }
00905    else
00906       BCLog::OutWarning(Form(
00907             "BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.",
00908             file));
00909 
00910    froot->Close();
00911 
00912    return r;
00913 }

int BCModel::ReadMarginalizedFromFile ( const char *  file  ) 

Read

Definition at line 822 of file BCModel.cxx.

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

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

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

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 Mon Sep 21 18:06:03 2009 for Bayesian Analysis Toolkit by  doxygen 1.6.1