#include <BCModel.h>
Inherits BCIntegrate.
Inherited by BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.
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 | FixDataAxis (unsigned int index, bool fixed) |
BCH2D * | GetMarginalized (const char *name1, const char *name2) |
BCH2D * | GetMarginalized (BCParameter *parameter1, BCParameter *parameter2) |
BCH1D * | GetMarginalized (const char *name) |
BCH1D * | GetMarginalized (BCParameter *parameter) |
double | GetPValue () |
double | GetPvalueFromChi2 (std::vector< double > par, int sigma_index) |
double | HessianMatrixElement (BCParameter *parameter1, BCParameter *parameter2, std::vector< double > point) |
double | Likelihood (std::vector< double > parameter) |
virtual double | LogAPrioriProbability (std::vector< double > parameters) |
virtual double | LogConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters) |
double | LogEval (std::vector< double > parameters) |
virtual double | LogLikelihood (std::vector< double > parameter) |
double | LogProbability (std::vector< double > parameter) |
double | LogProbabilityNN (std::vector< double > parameter) |
int | MarginalizeAll () |
double | Normalize () |
int | PrintAllMarginalized (const char *file, unsigned int hdiv=1, unsigned int ndiv=1) |
int | PrintAllMarginalized1D (const char *filebase) |
int | PrintAllMarginalized2D (const char *filebase) |
void | PrintHessianMatrix (std::vector< double > parameters) |
void | PrintResults (const char *file) |
void | PrintShortFitSummary () |
void | PrintSummary () |
double | Probability (std::vector< double > parameter) |
double | ProbabilityNN (std::vector< double > parameter) |
int | ReadErrorBandFromFile (const char *file) |
int | ReadMarginalizedFromFile (const char *file) |
int | ReadMode (const char *file) |
virtual double | SamplingFunction (std::vector< double > parameters) |
void | SetGoFNChains (int n) |
void | SetGoFNIterationsMax (int n) |
void | SetGoFNIterationsRun (int n) |
void | WriteMode (const char *file) |
Constructors and destructors | |
BCModel (const char *name) | |
BCModel () | |
virtual | ~BCModel () |
Member functions (get) | |
double | GetBestFitParameter (unsigned int index) |
double | GetBestFitParameterMarginalized (unsigned int index) |
std::vector< double > | GetBestFitParameters () |
std::vector< double > | GetBestFitParametersMarginalized () |
BCDataPoint * | GetDataPoint (unsigned int index) |
BCDataPoint * | GetDataPointLowerBoundaries () |
double | GetDataPointLowerBoundary (unsigned int index) |
BCDataPoint * | GetDataPointUpperBoundaries () |
double | GetDataPointUpperBoundary (unsigned int index) |
BCDataSet * | GetDataSet () |
std::vector< double > | GetErrorBand (double level) |
TGraph * | GetErrorBandGraph (double level1, double level2) |
TH2D * | GetErrorBandXY () |
TH2D * | GetErrorBandXY_yellow (double level=.68, int nsmooth=0) |
TGraph * | GetFitFunctionGraph (std::vector< double > parameters, double xmin, double xmax, int n=1000) |
TGraph * | GetFitFunctionGraph () |
TGraph * | GetFitFunctionGraph (std::vector< double > parameters) |
bool | GetFixedDataAxis (unsigned int index) |
bool | GetFlagBoundaries () |
int | GetIndex () |
double | GetModelAPosterioriProbability () |
double | GetModelAPrioriProbability () |
std::string | GetName () |
int | GetNDataPoints () |
unsigned int | GetNDataPointsMaximum () |
unsigned int | GetNDataPointsMinimum () |
double | GetNormalization () |
unsigned int | GetNParameters () |
BCParameter * | GetParameter (const char *name) |
BCParameter * | GetParameter (int index) |
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 | |
BCDataSet * | fDataSet |
int | fGoFNChains |
int | fGoFNIterationsMax |
int | fGoFNIterationsRun |
int | fIndex |
bool | flag_ConditionalProbabilityEntry |
double | fModelAPosteriori |
double | fModelAPriori |
std::string | fName |
unsigned int | fNDataPointsMaximum |
unsigned int | fNDataPointsMinimum |
BCParameterSet * | fParameterSet |
double | fPValue |
Private Member Functions | |
int | CompareStrings (const char *string1, const char *string2) |
BCDataPoint * | VectorToDataPoint (std::vector< double > data) |
Private Attributes | |
double | fNormalization |
Definition at line 50 of file BCModel.h.
BCModel::BCModel | ( | ) |
The default constructor.
Definition at line 61 of file BCModel.cxx.
00061 : BCIntegrate() 00062 { 00063 fNormalization = -1.0; 00064 fDataSet = 0; 00065 fParameterSet = new BCParameterSet(); 00066 00067 fIndex = -1; 00068 fPValue = -1; 00069 00070 fName = "model"; 00071 fDataPointUpperBoundaries = 0; 00072 fDataPointLowerBoundaries = 0; 00073 00074 flag_ConditionalProbabilityEntry = true; 00075 00076 fGoFNChains = 5; 00077 fGoFNIterationsMax = 100000; 00078 fGoFNIterationsRun = 2000; 00079 }
BCModel::BCModel | ( | const char * | name | ) |
A constructor.
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 00046 fName = (char *) name; 00047 flag_ConditionalProbabilityEntry = true; 00048 00049 fDataPointUpperBoundaries = 0; 00050 fDataPointLowerBoundaries = 0; 00051 00052 fErrorBandXY = 0; 00053 00054 fGoFNChains = 5; 00055 fGoFNIterationsMax = 100000; 00056 fGoFNIterationsRun = 2000; 00057 }
BCModel::~BCModel | ( | ) | [virtual] |
The default destructor.
Definition at line 83 of file BCModel.cxx.
00084 { 00085 delete fParameterSet; 00086 00087 if (fDataPointLowerBoundaries) 00088 delete fDataPointLowerBoundaries; 00089 00090 if (fDataPointUpperBoundaries) 00091 delete fDataPointUpperBoundaries; 00092 }
int BCModel::AddParameter | ( | BCParameter * | parameter | ) |
Adds a parameter to the model.
parameter | A model parameter |
Definition at line 436 of file BCModel.cxx.
00437 { 00438 // check if parameter set exists 00439 if (!fParameterSet) 00440 { 00441 BCLog::OutError("BCModel::AddParameter : Parameter set does not exist"); 00442 return ERROR_PARAMETERSETDOESNOTEXIST; 00443 } 00444 00445 // check if parameter with same name exists 00446 int flag_exists = 0; 00447 for (unsigned int i = 0; i < this -> GetNParameters(); i++) 00448 if (this -> CompareStrings(parameter -> GetName().data(), this -> GetParameter(i) -> GetName().data()) == 0) 00449 flag_exists = -1; 00450 00451 if (flag_exists < 0) 00452 { 00453 BCLog::OutError( 00454 Form("BCModel::AddParameter : Parameter with name %s exists already. ", parameter -> GetName().data())); 00455 return ERROR_PARAMETEREXISTSALREADY; 00456 } 00457 00458 // define index of new parameter 00459 unsigned int index = fParameterSet -> size(); 00460 parameter -> SetIndex(index); 00461 00462 // add parameter to parameter container 00463 fParameterSet -> push_back(parameter); 00464 00465 // add parameters to integation methods 00466 this -> SetParameters(fParameterSet); 00467 00468 return 0; 00469 }
int BCModel::AddParameter | ( | const char * | name, | |
double | lowerlimit, | |||
double | upperlimit | |||
) |
Adds a parameter to the parameter set
name | The name of the parameter | |
lowerlimit | The lower limit of the parameter values | |
upperlimit | The upper limit of the parameter values |
Definition at line 422 of file BCModel.cxx.
00423 { 00424 // create new parameter 00425 BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit); 00426 00427 int flag_ok = this -> AddParameter(parameter); 00428 if (flag_ok) 00429 delete parameter; 00430 00431 return flag_ok; 00432 }
double BCModel::APrioriProbability | ( | std::vector< double > | parameters | ) | [inline] |
Returns the prior probability.
parameters | A set of parameter values |
Definition at line 329 of file BCModel.h.
00330 { return exp( this->LogAPrioriProbability(parameters) ); };
BCH1D * BCModel::CalculatePValue | ( | std::vector< double > | par, | |
bool | flag_histogram = false | |||
) |
Definition at line 1172 of file BCModel.cxx.
01173 { 01174 BCH1D * hist = 0; 01175 01176 // print log 01177 BCLog::OutSummary("Do goodness-of-fit-test"); 01178 01179 // create model test 01180 BCGoFTest * goftest = new BCGoFTest("modeltest"); 01181 01182 // set this model as the model to be tested 01183 goftest -> SetTestModel(this); 01184 01185 // set the point in parameter space which is tested an initialize 01186 // the model testing 01187 if (!goftest -> SetTestPoint(par)) 01188 return 0; 01189 01190 // set parameters of the MCMC for the GoFTest 01191 goftest -> MCMCSetNChains(fGoFNChains); 01192 goftest -> MCMCSetNIterationsMax(fGoFNIterationsMax); 01193 goftest -> MCMCSetNIterationsRun(fGoFNIterationsRun); 01194 01195 // get p-value 01196 fPValue = goftest -> GetCalculatedPValue(flag_histogram); 01197 01198 // get histogram 01199 if (flag_histogram) 01200 { 01201 hist = new BCH1D(); 01202 hist -> SetHistogram(goftest -> GetHistogramLogProb()); 01203 } 01204 01205 // delete model test 01206 delete goftest; 01207 01208 // return histogram 01209 return hist; 01210 }
int BCModel::CheckParameters | ( | std::vector< double > | parameters | ) |
Checks if a set of parameters values is within the given range.
parameters | A set of parameter values |
Definition at line 565 of file BCModel.cxx.
00566 { 00567 // check if vectors are of equal size 00568 if (!parameters.size() == fParameterSet -> size()) 00569 return ERROR_INVALIDNUMBEROFPARAMETERS; 00570 00571 // check if parameters are within limits 00572 for (unsigned int i = 0; i < fParameterSet -> size(); i++) 00573 { 00574 BCParameter * modelparameter = fParameterSet -> at(i); 00575 00576 if (modelparameter -> GetLowerLimit() > parameters.at(i) || 00577 modelparameter -> GetUpperLimit() < parameters.at(i)) 00578 { 00579 BCLog::OutError( 00580 Form("BCModel::CheckParameters : Parameter %s not within limits.", fParameterSet -> at(i) -> GetName().data())); 00581 return ERROR_PARAMETERNOTWITHINRANGE; 00582 } 00583 } 00584 00585 return 0; 00586 }
int BCModel::CompareStrings | ( | const char * | string1, | |
const char * | string2 | |||
) | [private] |
Compares to strings
Definition at line 1577 of file BCModel.cxx.
01578 { 01579 int flag_same = 0; 01580 01581 if (strlen(string1) != strlen(string2)) 01582 return -1; 01583 01584 for (int i = 0; i < int(strlen(string1)); i++) 01585 if (string1[i] != string2[i]) 01586 flag_same = -1; 01587 01588 return flag_same; 01589 }
double BCModel::ConditionalProbabilityEntry | ( | BCDataPoint * | datapoint, | |
std::vector< double > | parameters | |||
) | [inline] |
Returns a conditional probability. Method needs to be overloaded by the user.
datapoint | A data point | |
parameters | A set of parameter values |
Definition at line 389 of file BCModel.h.
00390 { return exp( this->LogConditionalProbabilityEntry(datapoint, parameters) ); };
void BCModel::CorrelateDataPointValues | ( | std::vector< double > & | x | ) | [virtual] |
Constrains a data point
x | A vector of double |
Definition at line 1214 of file BCModel.cxx.
double BCModel::Eval | ( | std::vector< double > | parameters | ) | [inline, virtual] |
Overloaded function to evaluate integral.
Reimplemented from BCIntegrate.
Definition at line 411 of file BCModel.h.
00412 { return exp( this->LogEval(parameters) ); };
double BCModel::EvalSampling | ( | std::vector< double > | parameters | ) | [virtual] |
Overloaded function to evaluate integral.
Reimplemented from BCIntegrate.
Definition at line 524 of file BCModel.cxx.
00525 { 00526 return this -> SamplingFunction(parameters); 00527 }
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 590 of file BCModel.cxx.
00591 { 00592 // this implementation is CLEARLY not good we have to work on this. 00593 00594 BCLog::OutSummary(Form("Model \'%s\': Finding mode", this -> GetName().data())); 00595 00596 // synchronize parameters in BCIntegrate 00597 this -> SetParameters(fParameterSet); 00598 00599 switch(this -> GetOptimizationMethod()) 00600 { 00601 case BCIntegrate::kOptSimulatedAnnealing: 00602 BCLog::OutError("BCModel::FindMode : Simulated annaeling not yet implemented"); 00603 return; 00604 00605 case BCIntegrate::kOptMinuit: 00606 this -> FindModeMinuit(start); 00607 return; 00608 00609 case BCIntegrate::kOptMetropolis: 00610 this -> MarginalizeAll(); 00611 return; 00612 } 00613 00614 BCLog::OutError( 00615 Form("BCModel::FindMode : Invalid mode finding method: %d", 00616 this->GetOptimizationMethod())); 00617 00618 return; 00619 }
void BCModel::FixDataAxis | ( | unsigned int | index, | |
bool | fixed | |||
) |
Definition at line 1268 of file BCModel.cxx.
01269 { 01270 // check if index is within range 01271 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 01272 { 01273 BCLog::OutWarning("BCModel::FixDataAxis. Index out of range."); 01274 return; 01275 } 01276 01277 if (fDataFixedValues.size() == 0) 01278 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false); 01279 01280 fDataFixedValues[index] = fixed; 01281 }
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.
index | The index of the parameter. |
Definition at line 177 of file BCModel.h.
00178 { return fBestFitParameters[index]; };
double BCModel::GetBestFitParameterMarginalized | ( | unsigned int | index | ) | [inline] |
Returns the value of a particular parameter (defined by index) at the modes of the marginalized posterior pdfs.
index | The index of the parameter. |
Definition at line 192 of file BCModel.h.
00193 { 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.
Definition at line 184 of file BCModel.h.
00185 { return fBestFitParameters; };
std::vector<double> BCModel::GetBestFitParametersMarginalized | ( | ) | [inline] |
Returns the set of values of the parameters at the modes of the marginalized posterior pdfs.
Definition at line 199 of file BCModel.h.
00200 { return fBestFitParametersMarginalized; };
BCDataPoint * BCModel::GetDataPoint | ( | unsigned int | index | ) |
index | The index of the data point. |
Definition at line 112 of file BCModel.cxx.
00113 { 00114 if (fDataSet) 00115 return fDataSet -> GetDataPoint(index); 00116 00117 BCLog::OutWarning("BCModel::GetDataPoint. No data set defined."); 00118 return 0; 00119 }
BCDataPoint* BCModel::GetDataPointLowerBoundaries | ( | ) | [inline] |
Definition at line 108 of file BCModel.h.
00109 { return fDataPointLowerBoundaries; };
double BCModel::GetDataPointLowerBoundary | ( | unsigned int | index | ) | [inline] |
index | The index of the variable. |
Definition at line 119 of file BCModel.h.
00120 { return fDataPointLowerBoundaries -> GetValue(index); };
BCDataPoint* BCModel::GetDataPointUpperBoundaries | ( | ) | [inline] |
Definition at line 113 of file BCModel.h.
00114 { return fDataPointUpperBoundaries; };
double BCModel::GetDataPointUpperBoundary | ( | unsigned int | index | ) | [inline] |
index | The index of the variable. |
Definition at line 125 of file BCModel.h.
00126 { return fDataPointUpperBoundaries -> GetValue(index); };
BCDataSet* BCModel::GetDataSet | ( | ) | [inline] |
std::vector< double > BCModel::GetErrorBand | ( | double | level | ) |
Returns a vector of y-values at a certain probability level.
level | The level of probability |
Definition at line 180 of file BCModel.cxx.
00181 { 00182 std::vector <double> errorband; 00183 00184 if (!fErrorBandXY) 00185 return errorband; 00186 00187 int nx = fErrorBandXY -> GetNbinsX(); 00188 errorband.assign(nx, 0.0); 00189 00190 // loop over x and y bins 00191 for (int ix = 1; ix <= nx; ix++) 00192 { 00193 TH1D * temphist = fErrorBandXY -> ProjectionY("temphist", ix, ix); 00194 00195 int nprobSum = 1; 00196 double q[1]; 00197 double probSum[1]; 00198 probSum[0] = level; 00199 00200 temphist -> GetQuantiles(nprobSum, q, probSum); 00201 00202 errorband[ix-1] = q[0]; 00203 } 00204 00205 return errorband; 00206 }
TGraph * BCModel::GetErrorBandGraph | ( | double | level1, | |
double | level2 | |||
) |
Definition at line 210 of file BCModel.cxx.
00211 { 00212 if (!fErrorBandXY) 00213 return 0; 00214 00215 // define new graph 00216 int nx = fErrorBandXY -> GetNbinsX(); 00217 00218 TGraph * graph = new TGraph(2 * nx); 00219 graph -> SetFillStyle(1001); 00220 graph -> SetFillColor(kYellow); 00221 00222 // get error bands 00223 std::vector <double> ymin = this -> GetErrorBand(level1); 00224 std::vector <double> ymax = this -> GetErrorBand(level2); 00225 00226 for (int i = 0; i < nx; i++) 00227 { 00228 graph -> SetPoint(i, fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1), ymin.at(i)); 00229 graph -> SetPoint(nx + i, fErrorBandXY -> GetXaxis() -> GetBinCenter(nx - i), ymax.at(nx - i - 1)); 00230 } 00231 00232 return graph; 00233 }
TH2D* BCModel::GetErrorBandXY | ( | ) | [inline] |
Definition at line 204 of file BCModel.h.
00205 { return fErrorBandXY; };
TH2D * BCModel::GetErrorBandXY_yellow | ( | double | level = .68 , |
|
int | nsmooth = 0 | |||
) |
Definition at line 237 of file BCModel.cxx.
00238 { 00239 if (!fErrorBandXY) 00240 return 0; 00241 00242 int nx = fErrorBandXY -> GetNbinsX(); 00243 int ny = fErrorBandXY -> GetNbinsY(); 00244 00245 // copy existing histogram 00246 TH2D * hist_tempxy = (TH2D*) fErrorBandXY -> Clone(TString::Format("%s_sub_%f.2",fErrorBandXY->GetName(),level)); 00247 hist_tempxy -> Reset(); 00248 hist_tempxy -> SetFillColor(kYellow); 00249 00250 // loop over x bins 00251 for (int ix = 1; ix < nx; ix++) 00252 { 00253 BCH1D * hist_temp = new BCH1D(); 00254 00255 TH1D * hproj = fErrorBandXY -> ProjectionY("temphist", ix, ix); 00256 if(nsmooth>0) 00257 hproj->Smooth(nsmooth); 00258 00259 hist_temp -> SetHistogram(hproj); 00260 00261 TH1D * hist_temp_yellow = hist_temp -> GetSmallestIntervalHistogram(level); 00262 00263 for (int iy = 1; iy <= ny; ++iy) 00264 hist_tempxy -> SetBinContent(ix, iy, hist_temp_yellow -> GetBinContent(iy)); 00265 00266 delete hist_temp_yellow; 00267 delete hist_temp; 00268 } 00269 00270 return hist_tempxy; 00271 }
TGraph * BCModel::GetFitFunctionGraph | ( | std::vector< double > | parameters, | |
double | xmin, | |||
double | xmax, | |||
int | n = 1000 | |||
) |
Definition at line 302 of file BCModel.cxx.
00303 { 00304 // define new graph 00305 TGraph * graph = new TGraph(n+1); 00306 00307 double dx = (xmax-xmin)/(double)n; 00308 00309 // loop over x values 00310 for (int i = 0; i <= n; i++) 00311 { 00312 double x = (double)i*dx+xmin; 00313 std::vector <double> xvec; 00314 xvec.push_back(x); 00315 double y = this -> FitFunction(xvec, parameters); 00316 00317 xvec.clear(); 00318 00319 graph -> SetPoint(i, x, y); 00320 } 00321 00322 return graph; 00323 }
TGraph* BCModel::GetFitFunctionGraph | ( | ) | [inline] |
Definition at line 219 of file BCModel.h.
00220 { return this -> GetFitFunctionGraph(this -> GetBestFitParameters()); };
TGraph * BCModel::GetFitFunctionGraph | ( | std::vector< double > | parameters | ) |
Definition at line 275 of file BCModel.cxx.
00276 { 00277 if (!fErrorBandXY) 00278 return 0; 00279 00280 // define new graph 00281 int nx = fErrorBandXY -> GetNbinsX(); 00282 TGraph * graph = new TGraph(nx); 00283 00284 // loop over x values 00285 for (int i = 0; i < nx; i++) 00286 { 00287 double x = fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1); 00288 00289 std::vector <double> xvec; 00290 xvec.push_back(x); 00291 double y = this -> FitFunction(xvec, parameters); 00292 xvec.clear(); 00293 00294 graph -> SetPoint(i, x, y); 00295 } 00296 00297 return graph; 00298 }
bool BCModel::GetFixedDataAxis | ( | unsigned int | index | ) |
Definition at line 1285 of file BCModel.cxx.
01286 { 01287 // check if index is within range 01288 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 01289 { 01290 BCLog::OutWarning("BCModel::GetFixedDataAxis. Index out of range."); 01291 return false; 01292 } 01293 01294 return fDataFixedValues.at(index); 01295 }
bool BCModel::GetFlagBoundaries | ( | ) |
Definition at line 327 of file BCModel.cxx.
00328 { 00329 if (!fDataPointLowerBoundaries) 00330 return false; 00331 00332 if (!fDataPointUpperBoundaries) 00333 return false; 00334 00335 if (fDataPointLowerBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues()) 00336 return false; 00337 00338 if (fDataPointUpperBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues()) 00339 return false; 00340 00341 return true; 00342 }
int BCModel::GetIndex | ( | ) | [inline] |
BCH2D* BCModel::GetMarginalized | ( | const char * | name1, | |
const char * | name2 | |||
) | [inline] |
Definition at line 493 of file BCModel.h.
00494 { 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.
parameter1 | First parameter | |
parameter2 | Second parameter |
Definition at line 1104 of file BCModel.cxx.
01105 { 01106 // if(fMCMCH2Marginalized.size()==0) 01107 // { 01108 // BCLog::Out(BCLog::warning, BCLog::warning, 01109 // "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this."); 01110 // return 0; 01111 // } 01112 01113 int index1 = parameter1 -> GetIndex(); 01114 int index2 = parameter2 -> GetIndex(); 01115 01116 if (index1 == index2) 01117 { 01118 BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available."); 01119 return 0; 01120 } 01121 01122 if (index1 > index2) 01123 { 01124 int itmp = index1; 01125 index1 = index2; 01126 index2 = itmp; 01127 } 01128 01129 // get histogram 01130 TH2D * hist = this -> MCMCGetH2Marginalized(index1, index2); 01131 01132 if(hist==0) 01133 return 0; 01134 01135 BCH2D * hprob = new BCH2D(); 01136 01137 // set axis labels 01138 hist -> SetName(Form("hist_%s_%s_%s", this -> GetName().data(), parameter1 -> GetName().data(), parameter2 -> GetName().data())); 01139 hist -> SetXTitle(Form("%s", parameter1 -> GetName().data())); 01140 hist -> SetYTitle(Form("%s", parameter2 -> GetName().data())); 01141 hist -> SetStats(kFALSE); 01142 01143 double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) }; 01144 hprob->SetGlobalMode(gmode); 01145 01146 // set histogram 01147 hprob -> SetHistogram(hist); 01148 01149 return hprob; 01150 }
BCH1D* BCModel::GetMarginalized | ( | const char * | name | ) | [inline] |
Definition at line 482 of file BCModel.h.
00483 { 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.
parameter | Model parameter |
Definition at line 746 of file BCModel.cxx.
00747 { 00748 // if(fMCMCH1Marginalized.size()==0) 00749 // { 00750 // BCLog::Out(BCLog::warning, BCLog::warning, 00751 // "BCModel::GetMarginalized. MarginalizeAll() has to be run prior to this."); 00752 // return 0; 00753 // } 00754 00755 int index = parameter -> GetIndex(); 00756 00757 // get histogram 00758 TH1D * hist = this -> MCMCGetH1Marginalized(index); 00759 if(!hist) 00760 return 0; 00761 00762 BCH1D * hprob = new BCH1D(); 00763 00764 // set axis labels 00765 hist -> SetName(Form("hist_%s_%s", this -> GetName().data(), parameter -> GetName().data())); 00766 hist -> SetXTitle(parameter -> GetName().data()); 00767 hist -> SetYTitle(Form("p(%s|data)", parameter -> GetName().data())); 00768 hist -> SetStats(kFALSE); 00769 00770 // set histogram 00771 hprob -> SetHistogram(hist); 00772 00773 // set best fit parameter 00774 double bestfit = hprob -> GetMode(); 00775 00776 if (fBestFitParametersMarginalized.size() == 0) 00777 for (unsigned int i = 0; i < this -> GetNParameters(); i++) 00778 fBestFitParametersMarginalized.push_back(0.0); 00779 00780 fBestFitParametersMarginalized[index] = bestfit; 00781 00782 hprob->SetGlobalMode(fBestFitParameters.at(index)); 00783 00784 return hprob; 00785 }
double BCModel::GetModelAPosterioriProbability | ( | ) | [inline] |
Definition at line 93 of file BCModel.h.
00094 { return fModelAPosteriori; };
double BCModel::GetModelAPrioriProbability | ( | ) | [inline] |
Definition at line 88 of file BCModel.h.
00089 { return fModelAPriori; };
std::string BCModel::GetName | ( | ) | [inline] |
int BCModel::GetNDataPoints | ( | ) |
Definition at line 96 of file BCModel.cxx.
00097 { 00098 int npoints = 0; 00099 if (fDataSet) 00100 npoints = fDataSet -> GetNDataPoints(); 00101 else 00102 { 00103 BCLog::OutWarning("BCModel::GetNDataPoints(). No data set defined."); 00104 return ERROR_NOEVENTS; 00105 } 00106 00107 return npoints; 00108 }
unsigned int BCModel::GetNDataPointsMaximum | ( | ) | [inline] |
Definition at line 149 of file BCModel.h.
00150 { return fNDataPointsMaximum; };
unsigned int BCModel::GetNDataPointsMinimum | ( | ) | [inline] |
Definition at line 144 of file BCModel.h.
00145 { return fNDataPointsMinimum; };
double BCModel::GetNormalization | ( | ) | [inline] |
Definition at line 98 of file BCModel.h.
00099 { return fNormalization; };
unsigned int BCModel::GetNParameters | ( | ) | [inline] |
Definition at line 154 of file BCModel.h.
00155 { return fParameterSet ? fParameterSet -> size() : 0; };
BCParameter * BCModel::GetParameter | ( | const char * | name | ) |
name | The name of the parameter in the parameter set. |
Definition at line 140 of file BCModel.cxx.
00141 { 00142 if (!fParameterSet) 00143 return 0; 00144 00145 int index = -1; 00146 for (unsigned int i = 0; i < this->GetNParameters(); i++) 00147 if (name == this -> GetParameter(i) -> GetName()) 00148 index = i; 00149 00150 if (index<0) 00151 { 00152 BCLog::OutWarning( 00153 Form( 00154 "BCModel::GetParameter : Model %s has no parameter named '%s'", 00155 (this -> GetName()).data(), name 00156 ) 00157 ); 00158 return 0; 00159 } 00160 00161 return this->GetParameter(index); 00162 }
BCParameter * BCModel::GetParameter | ( | int | index | ) |
index | The index of the parameter in the parameter set. |
Definition at line 123 of file BCModel.cxx.
00124 { 00125 if (!fParameterSet) 00126 return 0; 00127 00128 if (index < 0 || index >= (int)this -> GetNParameters()) 00129 { 00130 BCLog::OutWarning( 00131 Form("BCModel::GetParameter. Parameter index %d not within range.", index)); 00132 return 0; 00133 } 00134 00135 return fParameterSet -> at(index); 00136 }
BCParameterSet* BCModel::GetParameterSet | ( | ) | [inline] |
double BCModel::GetPValue | ( | ) | [inline] |
double BCModel::GetPvalueFromChi2 | ( | std::vector< double > | par, | |
int | sigma_index | |||
) |
Calculate p-value from Chi2 distribution for Gaussian problems
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 1154 of file BCModel.cxx.
01155 { 01156 double ll = this -> LogLikelihood(par); 01157 int n = this -> GetNDataPoints(); 01158 01159 double sum_sigma=0; 01160 for (int i=0;i<n;i++) 01161 sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index)); 01162 01163 double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma); 01164 01165 fPValue = TMath::Prob(chi2,n); 01166 01167 return fPValue; 01168 }
double BCModel::HessianMatrixElement | ( | BCParameter * | parameter1, | |
BCParameter * | parameter2, | |||
std::vector< double > | point | |||
) |
Calculates the matrix element of the Hessian matrix
parameter1 | The parameter for the first derivative | |
parameter2 | The parameter for the first derivative |
Definition at line 1221 of file BCModel.cxx.
01222 { 01223 // check number of entries in vector 01224 if (point.size() != this -> GetNParameters()) 01225 { 01226 BCLog::OutWarning("BCModel::HessianMatrixElement. Invalid number of entries in the vector."); 01227 return -1; 01228 } 01229 01230 // define steps 01231 double nsteps = 1e5; 01232 double dx1 = par1 -> GetRangeWidth() / nsteps; 01233 double dx2 = par2 -> GetRangeWidth() / nsteps; 01234 01235 // define points at which to evaluate 01236 std::vector<double> xpp = point; 01237 std::vector<double> xpm = point; 01238 std::vector<double> xmp = point; 01239 std::vector<double> xmm = point; 01240 01241 int idx1 = par1 -> GetIndex(); 01242 int idx2 = par2 -> GetIndex(); 01243 01244 xpp[idx1] += dx1; 01245 xpp[idx2] += dx2; 01246 01247 xpm[idx1] += dx1; 01248 xpm[idx2] -= dx2; 01249 01250 xmp[idx1] -= dx1; 01251 xmp[idx2] += dx2; 01252 01253 xmm[idx1] -= dx1; 01254 xmm[idx2] -= dx2; 01255 01256 // calculate probability at these points 01257 double ppp = this -> Likelihood(xpp); 01258 double ppm = this -> Likelihood(xpm); 01259 double pmp = this -> Likelihood(xmp); 01260 double pmm = this -> Likelihood(xmm); 01261 01262 // return derivative 01263 return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2); 01264 }
double BCModel::Likelihood | ( | std::vector< double > | parameter | ) | [inline] |
Returns the likelihood
parameters | A set of parameter values |
Definition at line 345 of file BCModel.h.
00346 { 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 | A set of parameter values |
Reimplemented in BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.
Definition at line 338 of file BCModel.h.
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.
datapoint | A data point | |
parameters | A set of parameter values |
Definition at line 399 of file BCModel.h.
00400 { flag_ConditionalProbabilityEntry = false; return 0.; };
double BCModel::LogEval | ( | std::vector< double > | parameters | ) | [virtual] |
Overloaded function to evaluate integral.
Reimplemented from BCIntegrate.
Definition at line 517 of file BCModel.cxx.
00518 { 00519 return this -> LogProbabilityNN(parameters); 00520 }
double BCModel::LogLikelihood | ( | std::vector< double > | parameter | ) | [virtual] |
Calculates natural logarithm of the likelihood. Method needs to be overloaded by the user.
parameters | A set of parameter values |
Reimplemented in BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.
Definition at line 501 of file BCModel.cxx.
00502 { 00503 double logprob = 0.; 00504 00505 // add log of conditional probabilities event-by-event 00506 for (unsigned int i=0;i<fDataSet -> GetNDataPoints();i++) 00507 { 00508 BCDataPoint * datapoint = this -> GetDataPoint(i); 00509 logprob += this -> LogConditionalProbabilityEntry(datapoint, parameters); 00510 } 00511 00512 return logprob; 00513 }
double BCModel::LogProbability | ( | std::vector< double > | parameter | ) |
Returns natural logarithm of the a posteriori probability given a set of parameter values
parameters | A set of parameter values |
Definition at line 486 of file BCModel.cxx.
00487 { 00488 // check if normalized 00489 if (fNormalization<0. || fNormalization==0.) 00490 { 00491 BCLog::Out(BCLog::warning, BCLog::warning, 00492 "BCModel::LogProbability. Normalization not available or zero."); 00493 return 0.; 00494 } 00495 00496 return this -> LogProbabilityNN(parameters) - log(fNormalization); 00497 }
double BCModel::LogProbabilityNN | ( | std::vector< double > | parameter | ) |
Returns the natural logarithm of likelihood times prior probability given a set of parameter values
parameters | A set of parameter values |
Definition at line 473 of file BCModel.cxx.
00474 { 00475 // add log of conditional probability 00476 double logprob = this -> LogLikelihood(parameters); 00477 00478 // add log of prior probability 00479 logprob += this -> LogAPrioriProbability(parameters); 00480 00481 return logprob; 00482 }
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.
Definition at line 704 of file BCModel.cxx.
00705 { 00706 BCLog::OutSummary(Form("Running MCMC for model \'%s\'",this->GetName().data())); 00707 00708 // prepare function fitting 00709 double dx = 0.0; 00710 double dy = 0.0; 00711 00712 if (fFitFunctionIndexX >= 0) 00713 { 00714 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX)) 00715 / (double)fErrorBandNbinsX; 00716 00717 dy = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY)) 00718 / (double)fErrorBandNbinsY; 00719 00720 fErrorBandXY = new TH2D( 00721 TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "", 00722 fErrorBandNbinsX, 00723 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - .5 * dx, 00724 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + .5 * dx, 00725 fErrorBandNbinsY, 00726 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - .5 * dy, 00727 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + .5 * dy); 00728 fErrorBandXY -> SetStats(kFALSE); 00729 00730 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix) 00731 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy) 00732 fErrorBandXY -> SetBinContent(ix, iy, 0.0); 00733 } 00734 00735 this -> MCMCMetropolis(); 00736 this -> FindModeMCMC(); 00737 00738 // this -> PrintResults(Form("%s.txt", this -> GetName().data())); 00739 00740 return 1; 00741 }
double BCModel::Normalize | ( | ) |
Integrates over the un-normalized probability and updates fNormalization.
Definition at line 541 of file BCModel.cxx.
00542 { 00543 BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",this->GetName().data())); 00544 00545 unsigned int n = this -> GetNvar(); 00546 00547 // initialize BCIntegrate if not done already 00548 if (n == 0) 00549 { 00550 this->SetParameters(fParameterSet); 00551 n = this->GetNvar(); 00552 } 00553 00554 // integrate and get best fit parameters 00555 // maybe we have to remove the mode finding from here in the future 00556 fNormalization = this -> Integrate(); 00557 00558 BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization)); 00559 00560 return fNormalization; 00561 }
int BCModel::PrintAllMarginalized | ( | const char * | file, | |
unsigned int | hdiv = 1 , |
|||
unsigned int | ndiv = 1 | |||
) |
Definition at line 941 of file BCModel.cxx.
00942 { 00943 if(fMCMCH1Marginalized.size()==0 || (fMCMCH2Marginalized.size()==0 && this -> GetNParameters() > 1)) 00944 { 00945 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available."); 00946 return 0; 00947 } 00948 00949 // if there's only one parameter, we just want to call Print() 00950 if (fMCMCH1Marginalized.size()==1 && fMCMCH2Marginalized.size()==0) 00951 { 00952 BCParameter * a = this->GetParameter(0); 00953 this -> GetMarginalized(a) -> Print(file); 00954 return 1; 00955 } 00956 00957 int c_width=600; // default canvas width 00958 int c_height=600; // default canvas height 00959 00960 int type = 112; // landscape 00961 00962 if (hdiv > vdiv) 00963 { 00964 if(hdiv>3) 00965 { 00966 c_width=1000; 00967 c_height=700; 00968 } 00969 else 00970 { 00971 c_width=800; 00972 c_height=600; 00973 } 00974 } 00975 else if(hdiv < vdiv) 00976 { 00977 if(hdiv>3) 00978 { 00979 c_height=1000; 00980 c_width=700; 00981 } 00982 else 00983 { 00984 c_height=800; 00985 c_width=600; 00986 } 00987 type=111; 00988 } 00989 00990 // get number of parameters of the model and calculate number of plots 00991 int npar = this -> GetNParameters(); 00992 int nplots2d = npar * (npar-1)/2; 00993 int nplots = npar + nplots2d; 00994 00995 // give out warning if too many plots 00996 BCLog::OutSummary(Form( 00997 "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s", 00998 npar,nplots2d,nplots,file)); 00999 if(nplots>100) 01000 BCLog::OutDetail("This can take a while..."); 01001 01002 // setup the canvas and postscript file 01003 TCanvas * c = new TCanvas( "c","canvas",c_width,c_height); 01004 01005 TPostScript * ps = new TPostScript(file,type); 01006 01007 if(type==112) 01008 ps->Range(24,16); 01009 else 01010 ps->Range(16,24); 01011 01012 // draw all 1D distributions 01013 ps->NewPage(); 01014 c->cd(); 01015 c->Clear(); 01016 c->Divide(hdiv,vdiv); 01017 01018 int n=0; 01019 for(int i=0;i<npar;i++) 01020 { 01021 // if current page is full, swith to new page 01022 if(i!=0 && i%(hdiv*vdiv)==0) 01023 { 01024 c->Update(); 01025 ps->NewPage(); 01026 c->cd(); 01027 c->Clear(); 01028 c->Divide(hdiv,vdiv); 01029 } 01030 01031 // go to next pad 01032 c->cd(i%(hdiv*vdiv)+1); 01033 01034 BCParameter * a = this->GetParameter(i); 01035 this -> GetMarginalized(a) -> Draw(); 01036 n++; 01037 01038 if(n%100==0) 01039 BCLog::OutDetail(Form(" --> %d plots done",n)); 01040 } 01041 01042 c->Update(); 01043 01044 // draw all the 2D distributions 01045 ps->NewPage(); 01046 c->cd(); 01047 c->Clear(); 01048 c->Divide(hdiv,vdiv); 01049 01050 int k=0; 01051 for(int i=0;i<npar-1;i++) 01052 { 01053 for(int j=i+1;j<npar;j++) 01054 { 01055 // if current page is full, switch to new page 01056 if(k!=0 && k%(hdiv*vdiv)==0) 01057 { 01058 c->Update(); 01059 ps->NewPage(); 01060 c->cd(); 01061 c->Clear(); 01062 c->Divide(hdiv,vdiv); 01063 } 01064 01065 // go to next pad 01066 c->cd(k%(hdiv*vdiv)+1); 01067 01068 BCParameter * a = this->GetParameter(i); 01069 BCParameter * b = this->GetParameter(j); 01070 01071 double meana = (a -> GetLowerLimit() + a -> GetUpperLimit()) / 2.0; 01072 double deltaa = (a -> GetUpperLimit() - a -> GetLowerLimit()); 01073 if (deltaa <= 1e-7 * meana) 01074 continue; 01075 01076 double meanb = (b -> GetLowerLimit() + b -> GetUpperLimit()) / 2.0; 01077 double deltab = (b -> GetUpperLimit() - b -> GetLowerLimit()); 01078 if (deltab <= 1e-7 * meanb) 01079 continue; 01080 01081 this -> GetMarginalized(a,b) -> Draw(52); 01082 k++; 01083 01084 if((n+k)%100==0) 01085 BCLog::OutDetail(Form(" --> %d plots done",n+k)); 01086 } 01087 } 01088 01089 if( (n+k)>100 && (n+k)%100 != 0 ) 01090 BCLog::OutDetail(Form(" --> %d plots done",n+k)); 01091 01092 c->Update(); 01093 ps->Close(); 01094 01095 delete c; 01096 delete ps; 01097 01098 // return total number of drawn histograms 01099 return n+k; 01100 }
int BCModel::PrintAllMarginalized1D | ( | const char * | filebase | ) |
Definition at line 884 of file BCModel.cxx.
00885 { 00886 if(fMCMCH1Marginalized.size()==0) 00887 { 00888 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available."); 00889 return 0; 00890 } 00891 00892 int n=this->GetNParameters(); 00893 for(int i=0;i<n;i++) 00894 { 00895 BCParameter * a = this->GetParameter(i); 00896 this -> GetMarginalized(a) -> Print(Form("%s_1D_%s.ps",filebase,a->GetName().data())); 00897 } 00898 00899 return n; 00900 }
int BCModel::PrintAllMarginalized2D | ( | const char * | filebase | ) |
Definition at line 904 of file BCModel.cxx.
00905 { 00906 if(fMCMCH2Marginalized.size()==0) 00907 { 00908 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available."); 00909 return 0; 00910 } 00911 00912 int k=0; 00913 int n=this->GetNParameters(); 00914 for(int i=0;i<n-1;i++) 00915 { 00916 for(int j=i+1;j<n;j++) 00917 { 00918 BCParameter * a = this->GetParameter(i); 00919 BCParameter * b = this->GetParameter(j); 00920 00921 double meana = (a -> GetLowerLimit() + a -> GetUpperLimit()) / 2.0; 00922 double deltaa = (a -> GetUpperLimit() - a -> GetLowerLimit()); 00923 if (deltaa <= 1e-7 * meana) 00924 continue; 00925 00926 double meanb = (b -> GetLowerLimit() + b -> GetUpperLimit()) / 2.0; 00927 double deltab = (b -> GetUpperLimit() - b -> GetLowerLimit()); 00928 if (deltab <= 1e-7 * meanb) 00929 continue; 00930 00931 this -> GetMarginalized(a,b) -> Print(Form("%s_2D_%s_%s.ps",filebase,a->GetName().data(),b->GetName().data())); 00932 k++; 00933 } 00934 } 00935 00936 return k; 00937 }
void BCModel::PrintHessianMatrix | ( | std::vector< double > | parameters | ) |
Prints matrix elements of the Hessian matrix
parameters | The parameter values at which point to evaluate the matrix |
Definition at line 1533 of file BCModel.cxx.
01534 { 01535 // check number of entries in vector 01536 if (parameters.size() != this -> GetNParameters()) 01537 { 01538 BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector"); 01539 return; 01540 } 01541 01542 // print to screen 01543 std::cout 01544 << std::endl 01545 << " Hessian matrix elements: " << std::endl 01546 << " Point: "; 01547 01548 for (int i = 0; i < int(parameters.size()); i++) 01549 std::cout << parameters.at(i) << " "; 01550 std::cout << std::endl; 01551 01552 // loop over all parameter pairs 01553 for (unsigned int i = 0; i < this -> GetNParameters(); i++) 01554 for (unsigned int j = 0; j < i; j++) 01555 { 01556 // calculate Hessian matrix element 01557 double hessianmatrixelement = this -> HessianMatrixElement(fParameterSet -> at(i), 01558 fParameterSet -> at(j), parameters); 01559 01560 // print to screen 01561 std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl; 01562 } 01563 }
void BCModel::PrintResults | ( | const char * | file | ) |
Prints a summary of the Markov Chain Monte Carlo to a file.
Definition at line 1362 of file BCModel.cxx.
01363 { 01364 // print summary of Markov Chain Monte Carlo 01365 01366 // open file 01367 ofstream ofi(file); 01368 01369 // check if file is open 01370 if(!ofi.is_open()) 01371 { 01372 std::cerr << "Couldn't open file " << file <<std::endl; 01373 return; 01374 } 01375 01376 // number of parameters 01377 int npar = fParameterSet -> size(); 01378 01379 // check convergence 01380 bool flag_conv = ((this -> MCMCGetNIterationsConvergenceGlobal() > 0)?1:0); 01381 01382 ofi 01383 << std::endl 01384 << " -----------------------------------------------------" << std::endl 01385 << " Summary of the Markov Chain Monte Carlo run" << std::endl 01386 << " -----------------------------------------------------" << std::endl 01387 << std::endl; 01388 01389 if (!flag_conv) 01390 { 01391 ofi 01392 << " WARNING: the Markov Chain did not converge! Be" << std::endl 01393 << " cautious using the following results!" << std::endl 01394 << std::endl; 01395 } 01396 01397 ofi 01398 << " Model summary" << std::endl 01399 << " =============" << std::endl 01400 << " Model: " << fName.data() << std::endl 01401 << " Number of parameters: " << npar << std::endl 01402 << " List of Parameters and ranges:" << std::endl; 01403 for (int i = 0; i < npar; ++i) 01404 { 01405 ofi 01406 << " (" << i << ") Parameter \"" 01407 << fParameterSet -> at(i) -> GetName().data() << "\"" 01408 << ": " << fParameterSet -> at(i) -> GetLowerLimit() 01409 << " - " 01410 << fParameterSet -> at(i) -> GetUpperLimit() << std::endl; 01411 } 01412 ofi << std::endl; 01413 01414 if (flag_conv) 01415 { 01416 ofi 01417 << " Results of the marginalization" << std::endl 01418 << " ==============================" << std::endl 01419 << " List of parameters and properties of the marginalized" << std::endl 01420 << " distributions:" << std::endl; 01421 for (int i = 0; i < npar; ++i) 01422 { 01423 BCH1D * bch1d = this -> GetMarginalized(fParameterSet -> at(i)); 01424 01425 ofi 01426 << " (" << i << ") Parameter \"" 01427 << fParameterSet -> at(i) -> GetName().data() << "\"" << std::endl 01428 << " Mean +- RMS: " 01429 << std::setprecision(4) << bch1d -> GetMean() 01430 << " +- " 01431 << std::setprecision(4) << bch1d -> GetRMS() << std::endl 01432 << " Median +- sigma: " 01433 << std::setprecision(4) << bch1d -> GetMedian() 01434 << " + " << std::setprecision(4) << bch1d -> GetQuantile(0.84) - bch1d -> GetMedian() 01435 << " - " << std::setprecision(4) << bch1d -> GetMedian() - bch1d -> GetQuantile(0.16) << std::endl 01436 << " (Marginalized) mode: " << bch1d -> GetMode() << std::endl 01437 << " Smallest interval(s) containing 68% and local modes:" << std::endl; 01438 01439 std::vector <double> v; 01440 v = bch1d -> GetSmallestIntervals(0.68); 01441 int ninter = int(v.size()); 01442 01443 for (int j = 0; j < ninter; j+=5) 01444 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; 01445 } 01446 ofi << std::endl; 01447 } 01448 01449 ofi 01450 << " Results of the optimization" << std::endl 01451 << " ===========================" << std::endl 01452 << " Optimization algorithm used: "; 01453 switch(this -> GetOptimizationMethod()) 01454 { 01455 case BCIntegrate::kOptSimulatedAnnealing: 01456 ofi << " Simulated Annealing" << std::endl; 01457 break; 01458 case BCIntegrate::kOptMinuit: 01459 ofi << " Minuit" << std::endl; 01460 break; 01461 case BCIntegrate::kOptMetropolis: 01462 ofi << " MCMC " << std::endl; 01463 break; 01464 } 01465 01466 ofi << " List of parameters and global mode:" << std::endl; 01467 for (int i = 0; i < npar; ++i) 01468 ofi 01469 << " (" << i << ") Parameter \"" 01470 << fParameterSet -> at(i) -> GetName().data() << "\": " 01471 << fBestFitParameters.at(i) << std::endl; 01472 ofi << std::endl; 01473 01474 if (fPValue >= 0.) 01475 { 01476 ofi 01477 << " Results of the model test" << std::endl 01478 << " =========================" << std::endl 01479 << " p-value at global mode: " << fPValue << std::endl 01480 << std::endl; 01481 } 01482 01483 ofi 01484 << " Status of the MCMC" << std::endl 01485 << " ==================" << std::endl 01486 << " Convergence reached: " << ((flag_conv)?"yes":"no") << std::endl; 01487 01488 if (flag_conv) 01489 ofi << " Number of iterations until convergence: " << this -> MCMCGetNIterationsConvergenceGlobal() << std::endl; 01490 else 01491 ofi 01492 << " WARNING: the Markov Chain did not converge! Be\n" 01493 << " cautious using the following results!" << std::endl 01494 << std::endl; 01495 ofi 01496 << " Number of chains: " << this -> MCMCGetNChains() << std::endl 01497 << " Number of iterations of each chain: " << this -> MCMCGetNIterationsMax() << std::endl 01498 << std::endl; 01499 01500 ofi 01501 << " -----------------------------------------------------" << std::endl 01502 << std::endl 01503 << " Notes" << std::endl 01504 << " =====" << std::endl 01505 << " (i) Median +- sigma denotes the median, m, of the" << std::endl 01506 << " marginalized distribution and the intervals from" << std::endl 01507 << " m to the 16% and 84% quantiles." << std::endl 01508 << " -----------------------------------------------------" << std::endl; 01509 01510 // close file 01511 // ofi.close; 01512 }
void BCModel::PrintShortFitSummary | ( | ) |
Prints a short summary of the fit results on the screen.
Definition at line 1516 of file BCModel.cxx.
01517 { 01518 BCLog::OutSummary("---------------------------------------------------"); 01519 BCLog::OutSummary(Form("Fit summary for model \'%s\':",this -> GetName().data())); 01520 BCLog::OutSummary(Form(" Number of parameters = %i", this -> GetNParameters())); 01521 01522 BCLog::OutSummary(" Best fit parameters (global):"); 01523 for (unsigned int i = 0; i < this -> GetNParameters(); ++i) 01524 BCLog::OutSummary(Form(" %s = %.2lf", this -> GetParameter(i) -> GetName().data(), this -> GetBestFitParameter(i))); 01525 01526 BCLog::OutSummary(" Goodness-of-fit test:"); 01527 BCLog::OutSummary(Form(" p-value = %.2lf", this -> GetPValue())); 01528 BCLog::OutSummary("---------------------------------------------------"); 01529 }
void BCModel::PrintSummary | ( | ) |
Prints a summary on the screen.
Definition at line 1299 of file BCModel.cxx.
01300 { 01301 int nparameters = this -> GetNParameters(); 01302 01303 // model summary 01304 std::cout 01305 << std::endl 01306 << " ---------------------------------" << std::endl 01307 << " Model : " << fName.data() << std::endl 01308 << " ---------------------------------"<< std::endl 01309 << " Index : " << fIndex << std::endl 01310 << " Number of parameters : " << nparameters << std::endl 01311 << std::endl 01312 << " - Parameters : " << std::endl 01313 << std::endl; 01314 01315 // parameter summary 01316 for (int i=0; i<nparameters; i++) 01317 fParameterSet -> at(i) -> PrintSummary(); 01318 01319 // best fit parameters 01320 if (this -> GetBestFitParameters().size() > 0) 01321 { 01322 std::cout 01323 << std::endl 01324 << " - Best fit parameters :" << std::endl 01325 << std::endl; 01326 01327 for (int i=0; i<nparameters; i++) 01328 { 01329 std::cout 01330 << " " << fParameterSet -> at(i) -> GetName().data() 01331 << " = " << this -> GetBestFitParameter(i) 01332 << " (overall)" << std::endl; 01333 if ((int)fBestFitParametersMarginalized.size() == nparameters) 01334 std::cout 01335 << " " << fParameterSet -> at(i) -> GetName().data() 01336 << " = " << this -> GetBestFitParameterMarginalized(i) 01337 << " (marginalized)" << std::endl; 01338 } 01339 } 01340 01341 std::cout << std::endl; 01342 01343 // model testing 01344 if (fPValue >= 0) 01345 { 01346 double likelihood = this -> Likelihood(this -> GetBestFitParameters()); 01347 std::cout 01348 << " - Model testing:" << std::endl 01349 << std::endl 01350 << " p(data|lambda*) = " << likelihood << std::endl 01351 << " p-value = " << fPValue << std::endl 01352 << std::endl; 01353 } 01354 01355 // normalization 01356 if (fNormalization > 0) 01357 std::cout << " Normalization : " << fNormalization << std::endl; 01358 }
double BCModel::Probability | ( | std::vector< double > | parameter | ) | [inline] |
Returns the a posteriori probability given a set of parameter values
parameters | A set of parameter values |
Definition at line 373 of file BCModel.h.
00374 { 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 | A set of parameter values |
Definition at line 359 of file BCModel.h.
00360 { return exp( this->LogProbabilityNN(parameter) ); };
int BCModel::ReadErrorBandFromFile | ( | const char * | file | ) |
Read
Definition at line 853 of file BCModel.cxx.
00854 { 00855 TFile * froot = new TFile(file); 00856 if(!froot->IsOpen()) 00857 { 00858 BCLog::OutWarning(Form("BCModel::ReadErrorBandFromFile. Couldn't open file %s.",file)); 00859 return 0; 00860 } 00861 00862 int r=0; 00863 00864 TH2D * h2 = (TH2D*)froot->Get("errorbandxy"); 00865 if(h2) 00866 { 00867 h2->SetDirectory(0); 00868 h2->SetName(TString::Format("errorbandxy_%d",BCLog::GetHIndex())); 00869 this->SetErrorBandHisto(h2); 00870 r=1; 00871 } 00872 else 00873 BCLog::OutWarning(Form( 00874 "BCModel::ReadErrorBandFromFile : Couldn't read histogram \"errorbandxy\" from file %s.", 00875 file)); 00876 00877 froot->Close(); 00878 00879 return r; 00880 }
int BCModel::ReadMarginalizedFromFile | ( | const char * | file | ) |
Read
Definition at line 789 of file BCModel.cxx.
00790 { 00791 TFile * froot = new TFile(file); 00792 if(!froot->IsOpen()) 00793 { 00794 BCLog::OutError(Form("BCModel::ReadMarginalizedFromFile : Couldn't open file %s.",file)); 00795 return 0; 00796 } 00797 00798 // We reset the MCMCEngine here for the moment. 00799 // In the future maybe we only want to do this if the engine 00800 // wans't initialized at all or when there were some changes 00801 // in the model. 00802 // But maybe we want reset everything since we're overwriting 00803 // the marginalized distributions anyway. 00804 this -> MCMCInitialize(); 00805 00806 int k=0; 00807 int n=this->GetNParameters(); 00808 for(int i=0;i<n;i++) 00809 { 00810 BCParameter * a = this -> GetParameter(i); 00811 TKey * key = froot -> GetKey(TString::Format("hist_%s_%s", this -> GetName().data(), a -> GetName().data())); 00812 if(key) 00813 { 00814 TH1D * h1 = (TH1D*) key -> ReadObjectAny(TH1D::Class()); 00815 h1->SetDirectory(0); 00816 if(this->SetMarginalized(i,h1)) 00817 k++; 00818 } 00819 else 00820 BCLog::OutWarning(Form( 00821 "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s\" from file %s.", 00822 this -> GetName().data(), a -> GetName().data(), file)); 00823 } 00824 00825 for(int i=0;i<n-1;i++) 00826 { 00827 for(int j=i+1;j<n;j++) 00828 { 00829 BCParameter * a = this -> GetParameter(i); 00830 BCParameter * b = this -> GetParameter(j); 00831 TKey * key = froot -> GetKey(TString::Format("hist_%s_%s_%s", this -> GetName().data(), a -> GetName().data(), b -> GetName().data())); 00832 if(key) 00833 { 00834 TH2D * h2 = (TH2D*) key -> ReadObjectAny(TH2D::Class()); 00835 h2->SetDirectory(0); 00836 if(this->SetMarginalized(i,j,h2)) 00837 k++; 00838 } 00839 else 00840 BCLog::OutWarning(Form( 00841 "BCModel::ReadMarginalizedFromFile : Couldn't read histogram \"hist_%s_%s_%s\" from file %s.", 00842 this -> GetName().data(), a -> GetName().data(), b -> GetName().data(), file)); 00843 } 00844 } 00845 00846 froot->Close(); 00847 00848 return k; 00849 }
int BCModel::ReadMode | ( | const char * | file | ) |
Read mode from file created by WriteMode() call
Definition at line 664 of file BCModel.cxx.
00665 { 00666 ifstream ifi(file); 00667 if(!ifi.is_open()) 00668 { 00669 BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.",file)); 00670 return 0; 00671 } 00672 00673 int npar = fParameterSet -> size(); 00674 std::vector <double> mode; 00675 00676 int i=0; 00677 while (i<npar && !ifi.eof()) 00678 { 00679 double a; 00680 ifi>>a; 00681 mode.push_back(a); 00682 i++; 00683 } 00684 00685 if(i<npar) 00686 { 00687 BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.",file)); 00688 BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.",npar,i)); 00689 return 0; 00690 } 00691 00692 BCLog::OutSummary(Form("# Read in best fit parameters (mode) for model \'%s\' from file %s:",fName.data(),file)); 00693 this->SetMode(mode); 00694 for(int j=0 ; j<npar; j++) 00695 BCLog::OutSummary(Form("# -> Parameter %d : %s = %e", j, fParameterSet->at(j)->GetName().data(), fBestFitParameters[j])); 00696 00697 BCLog::OutWarning("# ! Best fit values obtained before this call will be overwritten !"); 00698 00699 return npar; 00700 }
double BCModel::SamplingFunction | ( | std::vector< double > | parameters | ) | [virtual] |
Sampling function used for importance sampling. Method needs to be overloaded by the user.
parameters | A set of parameter values |
Definition at line 531 of file BCModel.cxx.
00532 { 00533 double probability = 1.0; 00534 for (std::vector<BCParameter*>::const_iterator it = fParameterSet -> begin(); it != fParameterSet -> end(); ++it) 00535 probability *= 1.0 / ((*it) -> GetUpperLimit() - (*it) -> GetLowerLimit()); 00536 return probability; 00537 }
void BCModel::SetDataBoundaries | ( | unsigned int | index, | |
double | lowerboundary, | |||
double | upperboundary, | |||
bool | fixed = false | |||
) |
Definition at line 371 of file BCModel.cxx.
00372 { 00373 // check if data set exists 00374 if (!fDataSet) 00375 { 00376 BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first."); 00377 return; 00378 } 00379 00380 // check if index is within range 00381 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 00382 { 00383 BCLog::OutError("BCModel::SetDataBoundaries : Index out of range."); 00384 return; 00385 } 00386 00387 // check if boundary data points exist 00388 if (!fDataPointLowerBoundaries) 00389 fDataPointLowerBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues()); 00390 00391 if (!fDataPointUpperBoundaries) 00392 fDataPointUpperBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues()); 00393 00394 if (fDataFixedValues.size() == 0) 00395 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false); 00396 00397 // set boundaries 00398 fDataPointLowerBoundaries -> SetValue(index, lowerboundary); 00399 fDataPointUpperBoundaries -> SetValue(index, upperboundary); 00400 fDataFixedValues[index] = fixed; 00401 }
void BCModel::SetDataSet | ( | BCDataSet * | dataset | ) | [inline] |
Sets the data set.
dataset | A data set |
Definition at line 273 of file BCModel.h.
00274 { fDataSet = dataset; fNormalization = -1.0; };
void BCModel::SetErrorBandContinuous | ( | bool | flag | ) |
Sets the error band flag to continuous function
Definition at line 405 of file BCModel.cxx.
00406 { 00407 fErrorBandContinuous = flag; 00408 00409 if (flag) 00410 return; 00411 00412 // clear x-values 00413 fErrorBandX.clear(); 00414 00415 // copy data x-values 00416 for (unsigned int i = 0; i < fDataSet -> GetNDataPoints(); ++i) 00417 fErrorBandX.push_back(fDataSet -> GetDataPoint(i) -> GetValue(fFitFunctionIndexX)); 00418 }
void BCModel::SetGoFNChains | ( | int | n | ) | [inline] |
void BCModel::SetGoFNIterationsMax | ( | int | n | ) | [inline] |
void BCModel::SetGoFNIterationsRun | ( | int | n | ) | [inline] |
void BCModel::SetIndex | ( | int | index | ) | [inline] |
Sets the index of the model within the BCModelManager.
index | The index of the model |
Definition at line 240 of file BCModel.h.
00241 { fIndex = index; };
void BCModel::SetModelAPosterioriProbability | ( | double | probability | ) | [inline] |
Sets the a posteriori probability for a model.
model | The model | |
probability | The a posteriori probability |
Definition at line 260 of file BCModel.h.
00261 { fModelAPosteriori = probability; };
void BCModel::SetModelAPrioriProbability | ( | double | probability | ) | [inline] |
Sets the a priori probability for a model.
model | The model | |
probability | The a priori probability |
Definition at line 253 of file BCModel.h.
00254 { fModelAPriori = probability; };
void BCModel::SetName | ( | const char * | name | ) | [inline] |
void BCModel::SetNbins | ( | const char * | parname, | |
int | nbins | |||
) |
Set the number of bins for the marginalized distribution of a parameter.
parname | The name of the parameter in the parameter set | |
nbins | Number of bins (default = 100) |
Definition at line 166 of file BCModel.cxx.
00167 { 00168 BCParameter * p = this -> GetParameter(parname); 00169 if(!p) 00170 { 00171 BCLog::OutWarning(Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname)); 00172 return; 00173 } 00174 00175 this -> BCIntegrate::SetNbins(nbins, p -> GetIndex()); 00176 }
void BCModel::SetNDataPointsMaximum | ( | unsigned int | maximum | ) | [inline] |
Sets the maximum number of data points.
Definition at line 290 of file BCModel.h.
00291 { fNDataPointsMaximum = maximum; };
void BCModel::SetNDataPointsMinimum | ( | unsigned int | minimum | ) | [inline] |
Sets the minimum number of data points.
Definition at line 285 of file BCModel.h.
00286 { fNDataPointsMinimum = minimum; };
void BCModel::SetNormalization | ( | double | norm | ) | [inline] |
Sets the normalization of the likelihood. The normalization is the integral of likelihood over all parameters.
norm | The normalization of the likelihood |
Definition at line 267 of file BCModel.h.
00268 { fNormalization = norm; };
void BCModel::SetParameterSet | ( | BCParameterSet * | parset | ) | [inline] |
Set all parameters of the model using a BCParameterSet container.
Definition at line 246 of file BCModel.h.
00247 { fParameterSet = parset; };
void BCModel::SetSingleDataPoint | ( | BCDataSet * | dataset, | |
unsigned int | index | |||
) |
Definition at line 361 of file BCModel.cxx.
00362 { 00363 if (index < 0 || index > dataset -> GetNDataPoints()) 00364 return; 00365 00366 this -> SetSingleDataPoint(dataset -> GetDataPoint(index)); 00367 }
void BCModel::SetSingleDataPoint | ( | BCDataPoint * | datapoint | ) |
Sets a single data point as data set.
datapoint | A data point |
Definition at line 346 of file BCModel.cxx.
00347 { 00348 // create new data set consisting of a single data point 00349 BCDataSet * dataset = new BCDataSet(); 00350 00351 // add the data point 00352 dataset -> AddDataPoint(datapoint); 00353 00354 // set this new data set 00355 this -> SetDataSet(dataset); 00356 00357 }
BCDataPoint * BCModel::VectorToDataPoint | ( | std::vector< double > | data | ) | [private] |
Converts a vector of doubles into a BCDataPoint
Definition at line 1567 of file BCModel.cxx.
01568 { 01569 int sizeofvector = int(data.size()); 01570 BCDataPoint * datapoint = new BCDataPoint(sizeofvector); 01571 datapoint -> SetValues(data); 01572 return datapoint; 01573 }
void BCModel::WriteMode | ( | const char * | file | ) |
Write mode into file
Definition at line 623 of file BCModel.cxx.
00624 { 00625 ofstream ofi(file); 00626 if(!ofi.is_open()) 00627 { 00628 std::cerr<<"Couldn't open file "<<file<<std::endl; 00629 return; 00630 } 00631 00632 int npar = fParameterSet -> size(); 00633 for (int i=0; i<npar; i++) 00634 ofi<<fBestFitParameters.at(i)<<std::endl; 00635 00636 ofi<<std::endl; 00637 ofi<<"#######################################################################"<<std::endl; 00638 ofi<<"#"<<std::endl; 00639 ofi<<"# This file was created automatically by BCModel::WriteMode() call."<<std::endl; 00640 ofi<<"# It can be read in by call to BCModel::ReadMode()."<<std::endl; 00641 ofi<<"# Do not modify it unless you know what you're doing."<<std::endl; 00642 ofi<<"#"<<std::endl; 00643 ofi<<"#######################################################################"<<std::endl; 00644 ofi<<"#"<<std::endl; 00645 ofi<<"# Best fit parameters (mode) for model:"<<std::endl; 00646 ofi<<"# \'"<<fName.data()<<"\'"<<std::endl; 00647 ofi<<"#"<<std::endl; 00648 ofi<<"# Number of parameters: "<<npar<<std::endl; 00649 ofi<<"# Parameters ordered as above:"<<std::endl; 00650 00651 for (int i=0; i<npar; i++) 00652 { 00653 ofi<<"# "<<i<<": "; 00654 ofi<<fParameterSet->at(i)->GetName().data()<<" = "; 00655 ofi<<fBestFitParameters.at(i)<<std::endl; 00656 } 00657 00658 ofi<<"#"<<std::endl; 00659 ofi<<"########################################################################"<<std::endl; 00660 }
BCDataSet* BCModel::fDataSet [protected] |
int BCModel::fGoFNChains [protected] |
int BCModel::fGoFNIterationsMax [protected] |
int BCModel::fGoFNIterationsRun [protected] |
int BCModel::fIndex [protected] |
bool BCModel::flag_ConditionalProbabilityEntry [protected] |
double BCModel::fModelAPosteriori [protected] |
double BCModel::fModelAPriori [protected] |
std::string BCModel::fName [protected] |
unsigned int BCModel::fNDataPointsMaximum [protected] |
unsigned int BCModel::fNDataPointsMinimum [protected] |
double BCModel::fNormalization [private] |
BCParameterSet* BCModel::fParameterSet [protected] |
double BCModel::fPValue [protected] |