#include <BCModel.h>
Inherits BCIntegrate.
Inherited by BCEfficiencyFitter, BCGoFTest, BCGraphFitter, and BCHistogramFitter.
Inheritance diagram for BCModel:
Public Member Functions | |
Member functions (get) | |
| std::string | GetName () |
| int | GetIndex () |
| double | GetModelAPrioriProbability () |
| double | GetModelAPosterioriProbability () |
| double | GetNormalization () |
| BCDataSet * | GetDataSet () |
| BCDataPoint * | GetDataPointLowerBoundaries () |
| BCDataPoint * | GetDataPointUpperBoundaries () |
| double | GetDataPointLowerBoundary (unsigned int index) |
| double | GetDataPointUpperBoundary (unsigned int index) |
| bool | GetFlagBoundaries () |
| int | GetNDataPoints () |
| BCDataPoint * | GetDataPoint (unsigned int index) |
| unsigned int | GetNDataPointsMinimum () |
| unsigned int | GetNDataPointsMaximum () |
| unsigned int | GetNParameters () |
| BCParameter * | GetParameter (int index) |
| BCParameter * | GetParameter (const char *name) |
| BCParameterSet * | GetParameterSet () |
| double | GetBestFitParameter (unsigned int index) |
| double | GetBestFitParameterError (unsigned int index) |
| std::vector< double > | GetBestFitParameters () |
| std::vector< double > | GetBestFitParameterErrors () |
| double | GetBestFitParameterMarginalized (unsigned int index) |
| std::vector< double > | GetBestFitParametersMarginalized () |
| TH2D * | GetErrorBandXY () |
| TH2D * | GetErrorBandXY_yellow (double level=.68, int nsmooth=0) |
| std::vector< double > | GetErrorBand (double level) |
| TGraph * | GetErrorBandGraph (double level1, double level2) |
| TGraph * | GetFitFunctionGraph (std::vector< double > parameters) |
| TGraph * | GetFitFunctionGraph () |
| TGraph * | GetFitFunctionGraph (std::vector< double > parameters, double xmin, double xmax, int n=1000) |
| bool | GetFixedDataAxis (unsigned int index) |
Member functions (set) | |
| void | SetName (const char *name) |
| void | SetIndex (int index) |
| void | SetParameterSet (BCParameterSet *parset) |
| void | SetModelAPrioriProbability (double probability) |
| void | SetModelAPosterioriProbability (double probability) |
| void | SetNormalization (double norm) |
| void | SetDataSet (BCDataSet *dataset) |
| void | SetSingleDataPoint (BCDataPoint *datapoint) |
| void | SetSingleDataPoint (BCDataSet *dataset, unsigned int index) |
| void | SetNDataPointsMinimum (unsigned int minimum) |
| void | SetNDataPointsMaximum (unsigned int maximum) |
| void | SetDataBoundaries (unsigned int index, double lowerboundary, double upperboundary, bool fixed=false) |
| void | SetErrorBandContinuous (bool flag) |
| void | SetNbins (const char *parname, int nbins) |
Member functions (miscellaneous methods) | |
| int | AddParameter (const char *name, double lowerlimit, double upperlimit) |
| int | AddParameter (BCParameter *parameter) |
| double | APrioriProbability (std::vector< double > parameters) |
| virtual double | LogAPrioriProbability (std::vector< double > parameters) |
| double | Likelihood (std::vector< double > parameter) |
| virtual double | LogLikelihood (std::vector< double > parameter) |
| double | ProbabilityNN (std::vector< double > parameter) |
| double | LogProbabilityNN (std::vector< double > parameter) |
| double | Probability (std::vector< double > parameter) |
| double | LogProbability (std::vector< double > parameter) |
| double | ConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters) |
| virtual double | LogConditionalProbabilityEntry (BCDataPoint *datapoint, std::vector< double > parameters) |
| virtual double | SamplingFunction (std::vector< double > parameters) |
| double | Eval (std::vector< double > parameters) |
| double | LogEval (std::vector< double > parameters) |
| double | EvalSampling (std::vector< double > parameters) |
| double | Normalize () |
| int | CheckParameters (std::vector< double > parameters) |
| void | FindMode (std::vector< double > start=std::vector< double >(0)) |
| void | FindModeMinuit (std::vector< double > start=std::vector< double >(0), int printlevel=1) |
| void | WriteMode (const char *file) |
| int | ReadMode (const char *file) |
| int | ReadMarginalizedFromFile (const char *file) |
| int | ReadErrorBandFromFile (const char *file) |
| int | MarginalizeAll () |
| BCH1D * | GetMarginalized (BCParameter *parameter) |
| BCH1D * | GetMarginalized (const char *name) |
| BCH2D * | GetMarginalized (BCParameter *parameter1, BCParameter *parameter2) |
| BCH2D * | GetMarginalized (const char *name1, const char *name2) |
| int | PrintAllMarginalized1D (const char *filebase) |
| int | PrintAllMarginalized2D (const char *filebase) |
| int | PrintAllMarginalized (const char *file, unsigned int hdiv=1, unsigned int ndiv=1) |
| virtual void | CorrelateDataPointValues (std::vector< double > &x) |
| double | GetPvalueFromChi2 (std::vector< double > par, int sigma_index) |
| double | GetPvalueFromChi2Johnson (std::vector< double > par) |
| double | GetChi2Johnson (std::vector< double > par, const int nBins=-1) |
| double | GetAvalueFromChi2Johnson (TTree *tree, TH1D *distribution=0) |
| double | GetPvalueFromChi2NDoF (std::vector< double > par, int sigma_index) |
| BCH1D * | CalculatePValue (std::vector< double > par, bool flag_histogram=false) |
| double | GetPValue () |
| double | GetPValueNDoF () |
| double | GetChi2NDoF () |
| void | SetGoFNIterationsMax (int n) |
| void | SetGoFNIterationsRun (int n) |
| void | SetGoFNChains (int n) |
| double | HessianMatrixElement (BCParameter *parameter1, BCParameter *parameter2, std::vector< double > point) |
| void | PrintSummary () |
| void | PrintResults (const char *file) |
| void | PrintShortFitSummary (int chi2flag=0) |
| void | PrintHessianMatrix (std::vector< double > parameters) |
| void | FixDataAxis (unsigned int index, bool fixed) |
| virtual double | CDF (const std::vector< double > ¶meters, int index, bool lower=false) |
Protected Attributes | |
| int | fIndex |
| std::string | fName |
| double | fModelAPriori |
| double | fModelAPosteriori |
| BCParameterSet * | fParameterSet |
| BCDataSet * | fDataSet |
| unsigned int | fNDataPointsMinimum |
| unsigned int | fNDataPointsMaximum |
| bool | flag_ConditionalProbabilityEntry |
| double | fPValue |
| double | fChi2NDoF |
| double | fPValueNDoF |
| bool | flag_discrete |
| int | fGoFNIterationsMax |
| int | fGoFNIterationsRun |
| int | fGoFNChains |
Private Member Functions | |
| BCDataPoint * | VectorToDataPoint (std::vector< double > data) |
| int | CompareStrings (const char *string1, const char *string2) |
| int | NumberBins () |
Private Attributes | |
| double | fNormalization |
Definition at line 50 of file BCModel.h.
| BCModel::BCModel | ( | ) |
The default constructor.
Definition at line 70 of file BCModel.cxx.
00070 : BCIntegrate() 00071 { 00072 fNormalization = -1.0; 00073 fDataSet = 0; 00074 fParameterSet = new BCParameterSet(); 00075 00076 fIndex = -1; 00077 fPValue = -1; 00078 fPValueNDoF = -1; 00079 fChi2NDoF = -1; 00080 00081 fName = "model"; 00082 fDataPointUpperBoundaries = 0; 00083 fDataPointLowerBoundaries = 0; 00084 00085 flag_ConditionalProbabilityEntry = true; 00086 00087 fGoFNChains = 5; 00088 fGoFNIterationsMax = 100000; 00089 fGoFNIterationsRun = 2000; 00090 00091 flag_discrete=false; 00092 00093 }
| BCModel::BCModel | ( | const char * | name | ) |
A constructor.
| name | The name of the model |
Definition at line 42 of file BCModel.cxx.
00042 : BCIntegrate() 00043 { 00044 fNormalization = -1.0; 00045 fDataSet = 0; 00046 fParameterSet = new BCParameterSet; 00047 00048 fIndex = -1; 00049 fPValue = -1; 00050 fPValueNDoF = -1; 00051 fChi2NDoF = -1; 00052 00053 fName = (char *) name; 00054 flag_ConditionalProbabilityEntry = true; 00055 00056 fDataPointUpperBoundaries = 0; 00057 fDataPointLowerBoundaries = 0; 00058 00059 fErrorBandXY = 0; 00060 00061 fGoFNChains = 5; 00062 fGoFNIterationsMax = 100000; 00063 fGoFNIterationsRun = 2000; 00064 00065 flag_discrete=false; 00066 }
| BCModel::~BCModel | ( | ) | [virtual] |
The default destructor.
Definition at line 97 of file BCModel.cxx.
00098 { 00099 delete fParameterSet; 00100 00101 if (fDataPointLowerBoundaries) 00102 delete fDataPointLowerBoundaries; 00103 00104 if (fDataPointUpperBoundaries) 00105 delete fDataPointUpperBoundaries; 00106 }
| BCModel::BCModel | ( | ) |
The default constructor.
Definition at line 70 of file BCModel.cxx.
00070 : BCIntegrate() 00071 { 00072 fNormalization = -1.0; 00073 fDataSet = 0; 00074 fParameterSet = new BCParameterSet(); 00075 00076 fIndex = -1; 00077 fPValue = -1; 00078 fPValueNDoF = -1; 00079 fChi2NDoF = -1; 00080 00081 fName = "model"; 00082 fDataPointUpperBoundaries = 0; 00083 fDataPointLowerBoundaries = 0; 00084 00085 flag_ConditionalProbabilityEntry = true; 00086 00087 fGoFNChains = 5; 00088 fGoFNIterationsMax = 100000; 00089 fGoFNIterationsRun = 2000; 00090 00091 flag_discrete=false; 00092 00093 }
| BCModel::BCModel | ( | const char * | name | ) |
A constructor.
| name | The name of the model |
Definition at line 42 of file BCModel.cxx.
00042 : BCIntegrate() 00043 { 00044 fNormalization = -1.0; 00045 fDataSet = 0; 00046 fParameterSet = new BCParameterSet; 00047 00048 fIndex = -1; 00049 fPValue = -1; 00050 fPValueNDoF = -1; 00051 fChi2NDoF = -1; 00052 00053 fName = (char *) name; 00054 flag_ConditionalProbabilityEntry = true; 00055 00056 fDataPointUpperBoundaries = 0; 00057 fDataPointLowerBoundaries = 0; 00058 00059 fErrorBandXY = 0; 00060 00061 fGoFNChains = 5; 00062 fGoFNIterationsMax = 100000; 00063 fGoFNIterationsRun = 2000; 00064 00065 flag_discrete=false; 00066 }
| BCModel::~BCModel | ( | ) | [virtual] |
The default destructor.
Definition at line 97 of file BCModel.cxx.
00098 { 00099 delete fParameterSet; 00100 00101 if (fDataPointLowerBoundaries) 00102 delete fDataPointLowerBoundaries; 00103 00104 if (fDataPointUpperBoundaries) 00105 delete fDataPointUpperBoundaries; 00106 }
| std::string BCModel::GetName | ( | ) | [inline] |
| int BCModel::GetIndex | ( | ) | [inline] |
| double BCModel::GetModelAPrioriProbability | ( | ) | [inline] |
Definition at line 88 of file BCModel.h.
00089 { return fModelAPriori; };
| double BCModel::GetModelAPosterioriProbability | ( | ) | [inline] |
Definition at line 93 of file BCModel.h.
00094 { return fModelAPosteriori; };
| double BCModel::GetNormalization | ( | ) | [inline] |
Definition at line 98 of file BCModel.h.
00099 { return fNormalization; };
| BCDataSet* BCModel::GetDataSet | ( | ) | [inline] |
| BCDataPoint* BCModel::GetDataPointLowerBoundaries | ( | ) | [inline] |
Definition at line 108 of file BCModel.h.
00109 { return fDataPointLowerBoundaries; };
| BCDataPoint* BCModel::GetDataPointUpperBoundaries | ( | ) | [inline] |
Definition at line 113 of file BCModel.h.
00114 { return fDataPointUpperBoundaries; };
| 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); };
| 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); };
| bool BCModel::GetFlagBoundaries | ( | ) |
Definition at line 340 of file BCModel.cxx.
00341 { 00342 if (!fDataPointLowerBoundaries) 00343 return false; 00344 00345 if (!fDataPointUpperBoundaries) 00346 return false; 00347 00348 if (fDataPointLowerBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues()) 00349 return false; 00350 00351 if (fDataPointUpperBoundaries -> GetNValues() != fDataSet -> GetDataPoint(0) -> GetNValues()) 00352 return false; 00353 00354 return true; 00355 }
| int BCModel::GetNDataPoints | ( | ) |
Definition at line 110 of file BCModel.cxx.
00111 { 00112 int npoints = 0; 00113 if (fDataSet) 00114 npoints = fDataSet -> GetNDataPoints(); 00115 else 00116 { 00117 BCLog::OutWarning("BCModel::GetNDataPoints() : No data set defined."); 00118 return ERROR_NOEVENTS; 00119 } 00120 00121 return npoints; 00122 }
| BCDataPoint * BCModel::GetDataPoint | ( | unsigned int | index | ) |
| index | The index of the data point. |
Definition at line 126 of file BCModel.cxx.
00127 { 00128 if (fDataSet) 00129 return fDataSet -> GetDataPoint(index); 00130 00131 BCLog::OutWarning("BCModel::GetDataPoint. No data set defined."); 00132 return 0; 00133 }
| unsigned int BCModel::GetNDataPointsMinimum | ( | ) | [inline] |
Definition at line 144 of file BCModel.h.
00145 { return fNDataPointsMinimum; };
| unsigned int BCModel::GetNDataPointsMaximum | ( | ) | [inline] |
Definition at line 149 of file BCModel.h.
00150 { return fNDataPointsMaximum; };
| unsigned int BCModel::GetNParameters | ( | ) | [inline] |
Definition at line 154 of file BCModel.h.
00155 { return fParameterSet ? fParameterSet -> size() : 0; };
| BCParameter * BCModel::GetParameter | ( | int | index | ) |
| index | The index of the parameter in the parameter set. |
Definition at line 137 of file BCModel.cxx.
00138 { 00139 if (!fParameterSet) 00140 return 0; 00141 00142 if (index < 0 || index >= (int)this -> GetNParameters()) 00143 { 00144 BCLog::OutWarning( 00145 Form("BCModel::GetParameter. Parameter index %d not within range.", index)); 00146 return 0; 00147 } 00148 00149 return fParameterSet -> at(index); 00150 }
| BCParameter * BCModel::GetParameter | ( | const char * | name | ) |
| name | The name of the parameter in the parameter set. |
Definition at line 154 of file BCModel.cxx.
00155 { 00156 if (!fParameterSet) 00157 return 0; 00158 00159 int index = -1; 00160 for (unsigned int i = 0; i < this->GetNParameters(); i++) 00161 if (name == this -> GetParameter(i) -> GetName()) 00162 index = i; 00163 00164 if (index<0) 00165 { 00166 BCLog::OutWarning(Form( 00167 "BCModel::GetParameter : Model %s has no parameter named '%s'", 00168 (this -> GetName()).data(), name 00169 ) 00170 ); 00171 return 0; 00172 } 00173 00174 return this->GetParameter(index); 00175 }
| BCParameterSet* BCModel::GetParameterSet | ( | ) | [inline] |
| 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::GetBestFitParameterError | ( | unsigned int | index | ) | [inline] |
| 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 187 of file BCModel.h.
00188 { return fBestFitParameters; };
| std::vector<double> BCModel::GetBestFitParameterErrors | ( | ) | [inline] |
| 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 198 of file BCModel.h.
00199 { return fBestFitParametersMarginalized[index]; };
| 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 205 of file BCModel.h.
00206 { return fBestFitParametersMarginalized; };
| TH2D* BCModel::GetErrorBandXY | ( | ) | [inline] |
Definition at line 210 of file BCModel.h.
00211 { return fErrorBandXY; };
| TH2D * BCModel::GetErrorBandXY_yellow | ( | double | level = .68, |
|
| int | nsmooth = 0 | |||
| ) |
Definition at line 250 of file BCModel.cxx.
00251 { 00252 if (!fErrorBandXY) 00253 return 0; 00254 00255 int nx = fErrorBandXY -> GetNbinsX(); 00256 int ny = fErrorBandXY -> GetNbinsY(); 00257 00258 // copy existing histogram 00259 TH2D * hist_tempxy = (TH2D*) fErrorBandXY -> Clone(TString::Format("%s_sub_%f.2",fErrorBandXY->GetName(),level)); 00260 hist_tempxy -> Reset(); 00261 hist_tempxy -> SetFillColor(kYellow); 00262 00263 // loop over x bins 00264 for (int ix = 1; ix < nx; ix++) 00265 { 00266 BCH1D * hist_temp = new BCH1D(); 00267 00268 TH1D * hproj = fErrorBandXY -> ProjectionY("temphist", ix, ix); 00269 if(nsmooth>0) 00270 hproj->Smooth(nsmooth); 00271 00272 hist_temp -> SetHistogram(hproj); 00273 00274 TH1D * hist_temp_yellow = hist_temp -> GetSmallestIntervalHistogram(level); 00275 00276 for (int iy = 1; iy <= ny; ++iy) 00277 hist_tempxy -> SetBinContent(ix, iy, hist_temp_yellow -> GetBinContent(iy)); 00278 00279 delete hist_temp_yellow; 00280 delete hist_temp; 00281 } 00282 00283 return hist_tempxy; 00284 }
| 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 193 of file BCModel.cxx.
00194 { 00195 std::vector <double> errorband; 00196 00197 if (!fErrorBandXY) 00198 return errorband; 00199 00200 int nx = fErrorBandXY -> GetNbinsX(); 00201 errorband.assign(nx, 0.0); 00202 00203 // loop over x and y bins 00204 for (int ix = 1; ix <= nx; ix++) 00205 { 00206 TH1D * temphist = fErrorBandXY -> ProjectionY("temphist", ix, ix); 00207 00208 int nprobSum = 1; 00209 double q[1]; 00210 double probSum[1]; 00211 probSum[0] = level; 00212 00213 temphist -> GetQuantiles(nprobSum, q, probSum); 00214 00215 errorband[ix-1] = q[0]; 00216 } 00217 00218 return errorband; 00219 }
| TGraph * BCModel::GetErrorBandGraph | ( | double | level1, | |
| double | level2 | |||
| ) |
Definition at line 223 of file BCModel.cxx.
00224 { 00225 if (!fErrorBandXY) 00226 return 0; 00227 00228 // define new graph 00229 int nx = fErrorBandXY -> GetNbinsX(); 00230 00231 TGraph * graph = new TGraph(2 * nx); 00232 graph -> SetFillStyle(1001); 00233 graph -> SetFillColor(kYellow); 00234 00235 // get error bands 00236 std::vector <double> ymin = this -> GetErrorBand(level1); 00237 std::vector <double> ymax = this -> GetErrorBand(level2); 00238 00239 for (int i = 0; i < nx; i++) 00240 { 00241 graph -> SetPoint(i, fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1), ymin.at(i)); 00242 graph -> SetPoint(nx + i, fErrorBandXY -> GetXaxis() -> GetBinCenter(nx - i), ymax.at(nx - i - 1)); 00243 } 00244 00245 return graph; 00246 }
| TGraph * BCModel::GetFitFunctionGraph | ( | std::vector< double > | parameters | ) |
Definition at line 288 of file BCModel.cxx.
00289 { 00290 if (!fErrorBandXY) 00291 return 0; 00292 00293 // define new graph 00294 int nx = fErrorBandXY -> GetNbinsX(); 00295 TGraph * graph = new TGraph(nx); 00296 00297 // loop over x values 00298 for (int i = 0; i < nx; i++) 00299 { 00300 double x = fErrorBandXY -> GetXaxis() -> GetBinCenter(i + 1); 00301 00302 std::vector <double> xvec; 00303 xvec.push_back(x); 00304 double y = this -> FitFunction(xvec, parameters); 00305 xvec.clear(); 00306 00307 graph -> SetPoint(i, x, y); 00308 } 00309 00310 return graph; 00311 }
| 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, | |
| double | xmin, | |||
| double | xmax, | |||
| int | n = 1000 | |||
| ) |
Definition at line 315 of file BCModel.cxx.
00316 { 00317 // define new graph 00318 TGraph * graph = new TGraph(n+1); 00319 00320 double dx = (xmax-xmin)/(double)n; 00321 00322 // loop over x values 00323 for (int i = 0; i <= n; i++) 00324 { 00325 double x = (double)i*dx+xmin; 00326 std::vector <double> xvec; 00327 xvec.push_back(x); 00328 double y = this -> FitFunction(xvec, parameters); 00329 00330 xvec.clear(); 00331 00332 graph -> SetPoint(i, x, y); 00333 } 00334 00335 return graph; 00336 }
| bool BCModel::GetFixedDataAxis | ( | unsigned int | index | ) |
Definition at line 1560 of file BCModel.cxx.
01561 { 01562 // check if index is within range 01563 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 01564 { 01565 BCLog::OutWarning("BCModel::GetFixedDataAxis. Index out of range."); 01566 return false; 01567 } 01568 01569 return fDataFixedValues.at(index); 01570 }
| void BCModel::SetName | ( | const char * | name | ) | [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 246 of file BCModel.h.
00247 { fIndex = index; };
| void BCModel::SetParameterSet | ( | BCParameterSet * | parset | ) | [inline] |
Set all parameters of the model using a BCParameterSet container.
Definition at line 252 of file BCModel.h.
00253 { fParameterSet = parset; };
| 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 259 of file BCModel.h.
00260 { fModelAPriori = probability; };
| 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 266 of file BCModel.h.
00267 { fModelAPosteriori = probability; };
| 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 273 of file BCModel.h.
00274 { fNormalization = norm; };
| void BCModel::SetDataSet | ( | BCDataSet * | dataset | ) | [inline] |
Sets the data set.
| dataset | A data set |
Definition at line 279 of file BCModel.h.
00280 { fDataSet = dataset; fNormalization = -1.0; };
| void BCModel::SetSingleDataPoint | ( | BCDataPoint * | datapoint | ) |
Sets a single data point as data set.
| datapoint | A data point |
Definition at line 359 of file BCModel.cxx.
00360 { 00361 // create new data set consisting of a single data point 00362 BCDataSet * dataset = new BCDataSet(); 00363 00364 // add the data point 00365 dataset -> AddDataPoint(datapoint); 00366 00367 // set this new data set 00368 this -> SetDataSet(dataset); 00369 00370 }
| void BCModel::SetSingleDataPoint | ( | BCDataSet * | dataset, | |
| unsigned int | index | |||
| ) |
Definition at line 374 of file BCModel.cxx.
00375 { 00376 if (index < 0 || index > dataset -> GetNDataPoints()) 00377 return; 00378 00379 this -> SetSingleDataPoint(dataset -> GetDataPoint(index)); 00380 }
| 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::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::SetDataBoundaries | ( | unsigned int | index, | |
| double | lowerboundary, | |||
| double | upperboundary, | |||
| bool | fixed = false | |||
| ) |
Definition at line 384 of file BCModel.cxx.
00385 { 00386 // check if data set exists 00387 if (!fDataSet) 00388 { 00389 BCLog::OutError("BCModel::SetDataBoundaries : Need to define data set first."); 00390 return; 00391 } 00392 00393 // check if index is within range 00394 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 00395 { 00396 BCLog::OutError("BCModel::SetDataBoundaries : Index out of range."); 00397 return; 00398 } 00399 00400 // check if boundary data points exist 00401 if (!fDataPointLowerBoundaries) 00402 fDataPointLowerBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues()); 00403 00404 if (!fDataPointUpperBoundaries) 00405 fDataPointUpperBoundaries = new BCDataPoint(fDataSet -> GetDataPoint(0) -> GetNValues()); 00406 00407 if (fDataFixedValues.size() == 0) 00408 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false); 00409 00410 // set boundaries 00411 fDataPointLowerBoundaries -> SetValue(index, lowerboundary); 00412 fDataPointUpperBoundaries -> SetValue(index, upperboundary); 00413 fDataFixedValues[index] = fixed; 00414 }
| void BCModel::SetErrorBandContinuous | ( | bool | flag | ) |
Sets the error band flag to continuous function
Definition at line 418 of file BCModel.cxx.
00419 { 00420 fErrorBandContinuous = flag; 00421 00422 if (flag) 00423 return; 00424 00425 // clear x-values 00426 fErrorBandX.clear(); 00427 00428 // copy data x-values 00429 for (unsigned int i = 0; i < fDataSet -> GetNDataPoints(); ++i) 00430 fErrorBandX.push_back(fDataSet -> GetDataPoint(i) -> GetValue(fFitFunctionIndexX)); 00431 }
| 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 179 of file BCModel.cxx.
00180 { 00181 BCParameter * p = this -> GetParameter(parname); 00182 if(!p) 00183 { 00184 BCLog::OutWarning(Form("BCModel::SetNbins : Parameter '%s' not found so Nbins not set",parname)); 00185 return; 00186 } 00187 00188 this -> BCIntegrate::SetNbins(nbins, p -> GetIndex()); 00189 }
| 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 435 of file BCModel.cxx.
00436 { 00437 // create new parameter 00438 BCParameter * parameter = new BCParameter(name, lowerlimit, upperlimit); 00439 00440 int flag_ok = this -> AddParameter(parameter); 00441 if (flag_ok) 00442 delete parameter; 00443 00444 return flag_ok; 00445 }
| int BCModel::AddParameter | ( | BCParameter * | parameter | ) |
Adds a parameter to the model.
| parameter | A model parameter |
Definition at line 449 of file BCModel.cxx.
00450 { 00451 // check if parameter set exists 00452 if (!fParameterSet) 00453 { 00454 BCLog::OutError("BCModel::AddParameter : Parameter set does not exist"); 00455 return ERROR_PARAMETERSETDOESNOTEXIST; 00456 } 00457 00458 // check if parameter with same name exists 00459 int flag_exists = 0; 00460 for (unsigned int i = 0; i < this -> GetNParameters(); i++) 00461 if (this -> CompareStrings(parameter -> GetName().data(), this -> GetParameter(i) -> GetName().data()) == 0) 00462 flag_exists = -1; 00463 00464 if (flag_exists < 0) 00465 { 00466 BCLog::OutError( 00467 Form("BCModel::AddParameter : Parameter with name %s exists already. ", parameter -> GetName().data())); 00468 return ERROR_PARAMETEREXISTSALREADY; 00469 } 00470 00471 // define index of new parameter 00472 unsigned int index = fParameterSet -> size(); 00473 parameter -> SetIndex(index); 00474 00475 // add parameter to parameter container 00476 fParameterSet -> push_back(parameter); 00477 00478 // add parameters to integation methods 00479 this -> SetParameters(fParameterSet); 00480 00481 return 0; 00482 }
| double BCModel::APrioriProbability | ( | std::vector< double > | parameters | ) | [inline] |
| 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 344 of file BCModel.h.
| double BCModel::Likelihood | ( | std::vector< double > | parameter | ) | [inline] |
| 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 514 of file BCModel.cxx.
00515 { 00516 double logprob = 0.; 00517 00518 // add log of conditional probabilities event-by-event 00519 for (unsigned int i=0;i<fDataSet -> GetNDataPoints();i++) 00520 { 00521 BCDataPoint * datapoint = this -> GetDataPoint(i); 00522 logprob += this -> LogConditionalProbabilityEntry(datapoint, parameters); 00523 } 00524 00525 return logprob; 00526 }
| double BCModel::ProbabilityNN | ( | std::vector< double > | parameter | ) | [inline] |
| 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 486 of file BCModel.cxx.
00487 { 00488 // add log of conditional probability 00489 double logprob = this -> LogLikelihood(parameters); 00490 00491 // add log of prior probability 00492 logprob += this -> LogAPrioriProbability(parameters); 00493 00494 return logprob; 00495 }
| double BCModel::Probability | ( | std::vector< double > | parameter | ) | [inline] |
| 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 499 of file BCModel.cxx.
00500 { 00501 // check if normalized 00502 if (fNormalization<0. || fNormalization==0.) 00503 { 00504 BCLog::Out(BCLog::warning, BCLog::warning, 00505 "BCModel::LogProbability. Normalization not available or zero."); 00506 return 0.; 00507 } 00508 00509 return this -> LogProbabilityNN(parameters) - log(fNormalization); 00510 }
| 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 395 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 405 of file BCModel.h.
00406 { flag_ConditionalProbabilityEntry = false; return 0.; };
| 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 544 of file BCModel.cxx.
00545 { 00546 double probability = 1.0; 00547 for (std::vector<BCParameter*>::const_iterator it = fParameterSet -> begin(); it != fParameterSet -> end(); ++it) 00548 probability *= 1.0 / ((*it) -> GetUpperLimit() - (*it) -> GetLowerLimit()); 00549 return probability; 00550 }
| 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.
| double BCModel::LogEval | ( | std::vector< double > | parameters | ) | [virtual] |
Overloaded function to evaluate integral.
Reimplemented from BCIntegrate.
Definition at line 530 of file BCModel.cxx.
00531 { 00532 return this -> LogProbabilityNN(parameters); 00533 }
| double BCModel::EvalSampling | ( | std::vector< double > | parameters | ) | [virtual] |
Overloaded function to evaluate integral.
Reimplemented from BCIntegrate.
Definition at line 537 of file BCModel.cxx.
00538 { 00539 return this -> SamplingFunction(parameters); 00540 }
| double BCModel::Normalize | ( | ) |
Integrates over the un-normalized probability and updates fNormalization.
Definition at line 554 of file BCModel.cxx.
00555 { 00556 BCLog::OutSummary(Form("Model \'%s\': Normalizing probability",this->GetName().data())); 00557 00558 unsigned int n = this -> GetNvar(); 00559 00560 // initialize BCIntegrate if not done already 00561 if (n == 0) 00562 { 00563 this->SetParameters(fParameterSet); 00564 n = this->GetNvar(); 00565 } 00566 00567 // integrate and get best fit parameters 00568 // maybe we have to remove the mode finding from here in the future 00569 fNormalization = this -> Integrate(); 00570 00571 BCLog::OutDetail(Form(" --> Normalization factor : %.6g", fNormalization)); 00572 00573 return fNormalization; 00574 }
| 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 578 of file BCModel.cxx.
00579 { 00580 // check if vectors are of equal size 00581 if (!parameters.size() == fParameterSet -> size()) 00582 return ERROR_INVALIDNUMBEROFPARAMETERS; 00583 00584 // check if parameters are within limits 00585 for (unsigned int i = 0; i < fParameterSet -> size(); i++) 00586 { 00587 BCParameter * modelparameter = fParameterSet -> at(i); 00588 00589 if (modelparameter -> GetLowerLimit() > parameters.at(i) || 00590 modelparameter -> GetUpperLimit() < parameters.at(i)) 00591 { 00592 BCLog::OutError( 00593 Form("BCModel::CheckParameters : Parameter %s not within limits.", fParameterSet -> at(i) -> GetName().data())); 00594 return ERROR_PARAMETERNOTWITHINRANGE; 00595 } 00596 } 00597 00598 return 0; 00599 }
| 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 603 of file BCModel.cxx.
00604 { 00605 // this implementation is CLEARLY not good we have to work on this. 00606 00607 BCLog::OutSummary(Form("Model \'%s\': Finding mode", this -> GetName().data())); 00608 00609 // synchronize parameters in BCIntegrate 00610 this -> SetParameters(fParameterSet); 00611 00612 switch(this -> GetOptimizationMethod()) 00613 { 00614 case BCIntegrate::kOptSA: 00615 this -> FindModeSA(start); 00616 return; 00617 00618 case BCIntegrate::kOptMinuit: 00619 { 00620 int printlevel = -1; 00621 if (BCLog::GetLogLevelScreen() <= BCLog::detail) 00622 printlevel = 0; 00623 if (BCLog::GetLogLevelScreen() <= BCLog::debug) 00624 printlevel = 1; 00625 this -> BCIntegrate::FindModeMinuit(start, printlevel); 00626 return; 00627 } 00628 00629 case BCIntegrate::kOptMetropolis: 00630 this -> MarginalizeAll(); 00631 return; 00632 } 00633 00634 BCLog::OutError( 00635 Form("BCModel::FindMode : Invalid mode finding method: %d", 00636 this->GetOptimizationMethod())); 00637 00638 return; 00639 }
| 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.
| start | point in parameter space from which the mode finding is started. | |
| printlevel | The print level. |
Reimplemented from BCIntegrate.
Definition at line 643 of file BCModel.cxx.
00644 { 00645 // synchronize parameters in BCIntegrate 00646 this -> SetParameters(fParameterSet); 00647 00648 this -> BCIntegrate::FindModeMinuit(start,printlevel); 00649 }
| void BCModel::WriteMode | ( | const char * | file | ) |
Write mode into file
Definition at line 653 of file BCModel.cxx.
00654 { 00655 ofstream ofi(file); 00656 if(!ofi.is_open()) 00657 { 00658 std::cerr<<"Couldn't open file "<<file<<std::endl; 00659 return; 00660 } 00661 00662 int npar = fParameterSet -> size(); 00663 for (int i=0; i<npar; i++) 00664 ofi<<fBestFitParameters.at(i)<<std::endl; 00665 00666 ofi<<std::endl; 00667 ofi<<"#######################################################################"<<std::endl; 00668 ofi<<"#"<<std::endl; 00669 ofi<<"# This file was created automatically by BCModel::WriteMode() call."<<std::endl; 00670 ofi<<"# It can be read in by call to BCModel::ReadMode()."<<std::endl; 00671 ofi<<"# Do not modify it unless you know what you're doing."<<std::endl; 00672 ofi<<"#"<<std::endl; 00673 ofi<<"#######################################################################"<<std::endl; 00674 ofi<<"#"<<std::endl; 00675 ofi<<"# Best fit parameters (mode) for model:"<<std::endl; 00676 ofi<<"# \'"<<fName.data()<<"\'"<<std::endl; 00677 ofi<<"#"<<std::endl; 00678 ofi<<"# Number of parameters: "<<npar<<std::endl; 00679 ofi<<"# Parameters ordered as above:"<<std::endl; 00680 00681 for (int i=0; i<npar; i++) 00682 { 00683 ofi<<"# "<<i<<": "; 00684 ofi<<fParameterSet->at(i)->GetName().data()<<" = "; 00685 ofi<<fBestFitParameters.at(i)<<std::endl; 00686 } 00687 00688 ofi<<"#"<<std::endl; 00689 ofi<<"########################################################################"<<std::endl; 00690 }
| int BCModel::ReadMode | ( | const char * | file | ) |
Read mode from file created by WriteMode() call
Definition at line 694 of file BCModel.cxx.
00695 { 00696 ifstream ifi(file); 00697 if(!ifi.is_open()) 00698 { 00699 BCLog::OutError(Form("BCModel::ReadMode : Couldn't open file %s.",file)); 00700 return 0; 00701 } 00702 00703 int npar = fParameterSet -> size(); 00704 std::vector <double> mode; 00705 00706 int i=0; 00707 while (i<npar && !ifi.eof()) 00708 { 00709 double a; 00710 ifi>>a; 00711 mode.push_back(a); 00712 i++; 00713 } 00714 00715 if(i<npar) 00716 { 00717 BCLog::OutError(Form("BCModel::ReadMode : Couldn't read mode from file %s.",file)); 00718 BCLog::OutError(Form("BCModel::ReadMode : Expected %d parameters, found %d.",npar,i)); 00719 return 0; 00720 } 00721 00722 BCLog::OutSummary(Form("# Read in best fit parameters (mode) for model \'%s\' from file %s:",fName.data(),file)); 00723 this->SetMode(mode); 00724 for(int j=0 ; j<npar; j++) 00725 BCLog::OutSummary(Form("# -> Parameter %d : %s = %e", j, fParameterSet->at(j)->GetName().data(), fBestFitParameters[j])); 00726 00727 BCLog::OutWarning("# ! Best fit values obtained before this call will be overwritten !"); 00728 00729 return npar; 00730 }
| 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::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::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 734 of file BCModel.cxx.
00735 { 00736 BCLog::OutSummary(Form("Running MCMC for model \'%s\'",this->GetName().data())); 00737 00738 // prepare function fitting 00739 double dx = 0.0; 00740 double dy = 0.0; 00741 00742 if (fFitFunctionIndexX >= 0) 00743 { 00744 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX)) 00745 / (double)fErrorBandNbinsX; 00746 00747 dy = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) - fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY)) 00748 / (double)fErrorBandNbinsY; 00749 00750 fErrorBandXY = new TH2D( 00751 TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "", 00752 fErrorBandNbinsX, 00753 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - .5 * dx, 00754 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + .5 * dx, 00755 fErrorBandNbinsY, 00756 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - .5 * dy, 00757 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + .5 * dy); 00758 fErrorBandXY -> SetStats(kFALSE); 00759 00760 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix) 00761 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy) 00762 fErrorBandXY -> SetBinContent(ix, iy, 0.0); 00763 } 00764 00765 this -> MCMCMetropolis(); 00766 this -> FindModeMCMC(); 00767 00768 // PrintResults(Form("%s.txt", this -> GetName().data())); 00769 00770 return 1; 00771 }
| 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 776 of file BCModel.cxx.
00777 { 00778 if (!parameter) 00779 { 00780 BCLog::OutError("BCModel::GetMarginalized : Parameter does not exist."); 00781 return 0; 00782 } 00783 00784 int index = parameter -> GetIndex(); 00785 if (fMCMCFlagsFillHistograms[index] == false) 00786 { 00787 BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' not filled.",parameter->GetName().data())); 00788 return 0; 00789 } 00790 00791 // get histogram 00792 TH1D * hist = this -> MCMCGetH1Marginalized(index); 00793 if(!hist) 00794 return 0; 00795 00796 // set axis labels 00797 hist->SetName(Form("hist_%s_%s", GetName().data(), parameter->GetName().data())); 00798 hist->SetXTitle(parameter->GetName().data()); 00799 hist->SetYTitle(Form("p(%s|data)", parameter->GetName().data())); 00800 hist->SetStats(kFALSE); 00801 00802 // set histogram 00803 BCH1D * hprob = new BCH1D(); 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 < GetNParameters(); i++) 00811 fBestFitParametersMarginalized.push_back(0.); 00812 00813 fBestFitParametersMarginalized[index] = bestfit; 00814 00815 hprob->SetGlobalMode(fBestFitParameters.at(index)); 00816 00817 return hprob; 00818 }
| BCH1D* BCModel::GetMarginalized | ( | const char * | name | ) | [inline] |
Definition at line 495 of file BCModel.h.
00496 { return this -> GetMarginalized(this -> GetParameter(name)); };
| BCH2D * BCModel::GetMarginalized | ( | BCParameter * | parameter1, | |
| BCParameter * | parameter2 | |||
| ) |
If MarginalizeAll method was used, the individual marginalized distributions with respect to otwo parameters can be retrieved using this method.
| parameter1 | First parameter | |
| parameter2 | Second parameter |
Definition at line 1174 of file BCModel.cxx.
01175 { 01176 int index1 = par1->GetIndex(); 01177 int index2 = par2->GetIndex(); 01178 01179 if (fMCMCFlagsFillHistograms[index1]==false || fMCMCFlagsFillHistograms[index2]==false ) 01180 { 01181 BCLog::OutError(Form("BCModel::GetMarginalized : Distribuion for '%s' and/or '%s' not filled.", 01182 par1->GetName().data(),par2->GetName().data())); 01183 return 0; 01184 } 01185 01186 if (index1 == index2) 01187 { 01188 BCLog::OutError("BCModel::GetMarginalized : Provided parameters are identical. Distribution not available."); 01189 return 0; 01190 } 01191 01192 BCParameter * npar1 = par1; 01193 BCParameter * npar2 = par2; 01194 01195 if (index1>index2) 01196 { 01197 npar1 = par2; 01198 npar2 = par1; 01199 01200 int itmp = index1; 01201 index1 = index2; 01202 index2 = itmp; 01203 } 01204 01205 // get histogram 01206 TH2D * hist = this -> MCMCGetH2Marginalized(index1, index2); 01207 01208 if(hist==0) 01209 return 0; 01210 01211 BCH2D * hprob = new BCH2D(); 01212 01213 // set axis labels 01214 hist->SetName(Form("hist_%s_%s_%s", GetName().data(), npar1->GetName().data(), npar2->GetName().data())); 01215 hist->SetXTitle(Form("%s", npar1->GetName().data())); 01216 hist->SetYTitle(Form("%s", npar2->GetName().data())); 01217 hist->SetStats(kFALSE); 01218 01219 double gmode[] = { fBestFitParameters.at(index1), fBestFitParameters.at(index2) }; 01220 hprob->SetGlobalMode(gmode); 01221 01222 // set histogram 01223 hprob->SetHistogram(hist); 01224 01225 return hprob; 01226 }
| 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)); };
| 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 }
| 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(!fMCMCFlagFillHistograms) 00979 { 00980 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not filled."); 00981 return 0; 00982 } 00983 00984 int npar = GetNParameters(); 00985 00986 if(fMCMCH1Marginalized.size()==0 || (fMCMCH2Marginalized.size()==0 && npar>1)) 00987 { 00988 BCLog::OutError("BCModel::PrintAllMarginalized : Marginalized distributions not available."); 00989 return 0; 00990 } 00991 00992 // if there's only one parameter, we just want to call Print() 00993 if (fMCMCH1Marginalized.size()==1 && fMCMCH2Marginalized.size()==0) 00994 { 00995 BCParameter * a = GetParameter(0); 00996 GetMarginalized(a)->Print(file); 00997 return 1; 00998 } 00999 01000 int c_width=600; // default canvas width 01001 int c_height=600; // default canvas height 01002 01003 int type = 112; // landscape 01004 01005 if (hdiv > vdiv) 01006 { 01007 if(hdiv>3) 01008 { 01009 c_width=1000; 01010 c_height=700; 01011 } 01012 else 01013 { 01014 c_width=800; 01015 c_height=600; 01016 } 01017 } 01018 else if(hdiv < vdiv) 01019 { 01020 if(hdiv>3) 01021 { 01022 c_height=1000; 01023 c_width=700; 01024 } 01025 else 01026 { 01027 c_height=800; 01028 c_width=600; 01029 } 01030 type=111; 01031 } 01032 01033 // calculate number of plots 01034 int nplots2d = npar * (npar-1)/2; 01035 int nplots = npar + nplots2d; 01036 01037 // give out warning if too many plots 01038 BCLog::OutSummary(Form( 01039 "Printing all marginalized distributions (%d x 1D + %d x 2D = %d) into file %s", 01040 npar,nplots2d,nplots,file)); 01041 if(nplots>100) 01042 BCLog::OutDetail("This can take a while..."); 01043 01044 // setup the canvas and postscript file 01045 TCanvas * c = new TCanvas( "c","canvas",c_width,c_height); 01046 01047 TPostScript * ps = new TPostScript(file,type); 01048 01049 if(type==112) 01050 ps->Range(24,16); 01051 else 01052 ps->Range(16,24); 01053 01054 // draw all 1D distributions 01055 ps->NewPage(); 01056 c->cd(); 01057 c->Clear(); 01058 c->Divide(hdiv,vdiv); 01059 01060 int n=0; 01061 for(int i=0;i<npar;i++) 01062 { 01063 // get corresponding parameter 01064 BCParameter * a = GetParameter(i); 01065 01066 // check if histogram exists 01067 if ( !GetMarginalized(a) ) 01068 continue; 01069 01070 // check if histogram is filled 01071 if ( GetMarginalized(a)->GetHistogram()->Integral() <= 0) 01072 continue; 01073 01074 // if current page is full, switch to new page 01075 if(i!=0 && i%(hdiv*vdiv)==0) 01076 { 01077 c->Update(); 01078 ps->NewPage(); 01079 c->cd(); 01080 c->Clear(); 01081 c->Divide(hdiv,vdiv); 01082 } 01083 01084 // go to next pad 01085 c->cd(i%(hdiv*vdiv)+1); 01086 01087 this -> GetMarginalized(a) -> Draw(); 01088 n++; 01089 01090 if(n%100==0) 01091 BCLog::OutDetail(Form(" --> %d plots done",n)); 01092 } 01093 01094 if (n > 0) 01095 { 01096 c->Update(); 01097 01098 // draw all the 2D distributions 01099 ps->NewPage(); 01100 c->cd(); 01101 c->Clear(); 01102 } 01103 01104 c->Divide(hdiv,vdiv); 01105 01106 int k=0; 01107 for(int i=0;i<npar-1;i++) 01108 { 01109 if(!fMCMCFlagsFillHistograms[i]) 01110 continue; 01111 for(int j=i+1;j<npar;j++) 01112 { 01113 if(!fMCMCFlagsFillHistograms[j]) 01114 continue; 01115 01116 // get corresponding parameters 01117 BCParameter * a = GetParameter(i); 01118 BCParameter * b = GetParameter(j); 01119 01120 // check if histogram exists 01121 if ( !GetMarginalized(a,b) ) 01122 continue; 01123 01124 // check if histogram is filled 01125 if ( GetMarginalized(a,b)->GetHistogram()->Integral() <= 0 ) 01126 continue; 01127 01128 // if current page is full, switch to new page 01129 if(k!=0 && k%(hdiv*vdiv)==0) 01130 { 01131 c->Update(); 01132 ps->NewPage(); 01133 c->cd(); 01134 c->Clear(); 01135 c->Divide(hdiv,vdiv); 01136 } 01137 01138 // go to next pad 01139 c->cd(k%(hdiv*vdiv)+1); 01140 01141 double meana = (a->GetLowerLimit() + a->GetUpperLimit()) / 2.; 01142 double deltaa = (a->GetUpperLimit() - a->GetLowerLimit()); 01143 if (deltaa <= 1e-7 * meana) 01144 continue; 01145 01146 double meanb = (b->GetLowerLimit() + b->GetUpperLimit()) / 2.; 01147 double deltab = (b->GetUpperLimit() - b->GetLowerLimit()); 01148 if (deltab <= 1e-7 * meanb) 01149 continue; 01150 01151 GetMarginalized(a,b)->Draw(52); 01152 k++; 01153 01154 if((n+k)%100==0) 01155 BCLog::OutDetail(Form(" --> %d plots done",n+k)); 01156 } 01157 } 01158 01159 if( (n+k)>100 && (n+k)%100 != 0 ) 01160 BCLog::OutDetail(Form(" --> %d plots done",n+k)); 01161 01162 c->Update(); 01163 ps->Close(); 01164 01165 delete c; 01166 delete ps; 01167 01168 // return total number of drawn histograms 01169 return n+k; 01170 }
| void BCModel::CorrelateDataPointValues | ( | std::vector< double > & | x | ) | [virtual] |
Constrains a data point
| x | A vector of double |
Definition at line 1489 of file BCModel.cxx.
| 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 1230 of file BCModel.cxx.
01231 { 01232 double ll = this -> LogLikelihood(par); 01233 int n = this -> GetNDataPoints(); 01234 01235 double sum_sigma=0; 01236 for (int i=0;i<n;i++) 01237 sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index)); 01238 01239 double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma); 01240 01241 fPValue = TMath::Prob(chi2,n); 01242 01243 return fPValue; 01244 }
| double BCModel::GetPvalueFromChi2Johnson | ( | std::vector< double > | par | ) |
Calculate p-value from asymptotic Chi2 distribution for arbitrary problems using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).
| par | Parameter set for the calculation of the likelihood |
Definition at line 1247 of file BCModel.cxx.
01247 { 01248 double chi2 = GetChi2Johnson(par); 01249 // look up corresponding p value 01250 fPValue = TMath::Prob(chi2, NumberBins() - 1); 01251 return fPValue; 01252 }
| double BCModel::GetChi2Johnson | ( | std::vector< double > | par, | |
| const int | nBins = -1 | |||
| ) |
Calculate Chi2 (also called R^{B}) for arbitrary problems with binned data using the definition (3) from Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).
| par | Parameter set for the calculation of the likelihood | |
| nBins | how many bins to use for the data, for negative an adapted rule \ of thumb by Wald(1942) is used, with at least three bins |
Definition at line 1255 of file BCModel.cxx.
01256 { 01257 01258 typedef unsigned int uint; 01259 01260 // number of observations 01261 int n = this -> GetNDataPoints(); 01262 01263 if (nBins < 0) 01264 nBins = NumberBins(); 01265 01266 // fixed width quantiles, including final point! 01267 std::vector<double> a; 01268 for (int i = 0; i <= nBins; i++) { 01269 a.push_back(i * 1.0 / double(nBins)); 01270 } 01271 01272 // determine the bin counts and fill the histogram with data using the CDF 01273 TH1I* hist = new TH1I("h1", "h1 title", nBins, 0.0, 1.0); 01274 01275 //discrete case requires randomization to allocate counts of bins that cover more 01276 // than one quantile 01277 if (flag_discrete) { 01278 //loop over observations, each may have different likelihood and CDF 01279 for (int j = 0; j < n; j++) { 01280 //actual value 01281 double CDF = this->CDF(par, j); 01282 //for the bin just before 01283 double CDFlower = this->CDF(par, j, true); 01284 01285 //what quantiles q are covered, count from q_1 to q_{nBins} 01286 int qMax = 1; 01287 for (int i = 0; i < nBins; i++) { 01288 if (CDF > a.at(i)) 01289 qMax = i + 1; 01290 else 01291 break; 01292 } 01293 int qMin = 1; 01294 for (int i = 0; i < nBins; i++) { 01295 if (CDFlower > a.at(i)) 01296 qMin = i + 1; 01297 else 01298 break; 01299 } 01300 01301 // simplest case: observation bin entirely contained in one quantile 01302 if (qMin == qMax) { 01303 //watch out for overflow because CDF exactly 1 01304 if (CDF < 1) 01305 hist -> Fill(CDF); 01306 else 01307 hist -> AddBinContent(qMax); 01308 01309 continue;//this observation finished 01310 } 01311 01312 // if more than quantile is covered need more work: 01313 //determine probabilities of this observation to go for each quantile covered 01314 // as follows: If each quantile has size 0.25 and the CDF(integral of likelihood) 01315 //for current observation gives gives 0.27, but for observation-1 we would have 01316 // 0.20, then 5/7 of the 7% go for first quantile and 2/7 for the second. 01317 //This extend to bins covering more than two quantiles 01318 std::vector<double> prob; 01319 //normalization 01320 double norm = 1 / double(CDF - CDFlower); 01321 01322 for (int i = 0; i < (qMax - qMin + 1); i++) { 01323 if (i == 0) { 01324 prob.push_back(norm * (a.at(qMin) - CDFlower)); 01325 continue; 01326 } 01327 if (i == (qMax - qMin)) { 01328 prob.push_back(norm * (CDF - a.at(qMax - 1))); 01329 continue; 01330 } 01331 //default case 01332 prob.push_back(norm * (a.at(i) - a.at(i - 1))); 01333 } 01334 //have distribution, use inverse-transform method 01335 double U = fRandom -> Rndm(); 01336 //build up the integral (CDF) 01337 for (uint i = 1; i < prob.size(); i++) 01338 prob.at(i) += prob.at(i - 1); 01339 // and search with linear comput. complexity 01340 for (uint i = 0; i < prob.size(); i++) { 01341 //we finally allocate the count, as center of quantile 01342 if (U <= prob.at(i)) { 01343 hist -> Fill((a.at(qMin + i - 1) + a.at(qMin + i)) / 2.0); 01344 break; 01345 } 01346 } 01347 } 01348 } else { //continuous case is simple 01349 for (int j = 0; j < n; j++) { 01350 hist -> Fill(this->CDF(par, j)); 01351 } 01352 } 01353 01354 // calculate chi^2 01355 double chi2 = 0.0; 01356 double mk, pk; 01357 double N = double(n); 01358 for (int i = 1; i <= nBins; i++) { 01359 mk = hist->GetBinContent(i); 01360 pk = a.at(i) - a.at(i - 1); 01361 chi2 += (mk - N * pk) * (mk - N * pk) / (N * pk); 01362 } 01363 01364 delete hist; 01365 01366 return chi2; 01367 01368 }
| double BCModel::GetAvalueFromChi2Johnson | ( | TTree * | tree, | |
| TH1D * | distribution = 0 | |||
| ) |
Calculate the A-value, a summary statistic. It computes the frequency that a Chi2 value determined from the data by Johnson's binning prescription is larger than a value sampled from the reference chi2 distribution. They out from one chain is used. A=1/2 provides no evidence against the null hypothesis. Compare Johnson, V.E. A Bayesian chi2 Test for Goodness-of-Fit. The Annals of Statistics 32, 2361-2384(2004).
| par | tree contains the samples of posterior of the parameters | |
| par | histogram filled by function with distribution of p values |
Definition at line 1370 of file BCModel.cxx.
01370 { 01371 01372 //model parameters 01373 int nPar = (int)this->GetNParameters(); 01374 std::vector<double> param(nPar); 01375 01376 //parameters saved in branches should be the same 01377 int nParBranches=-1; 01378 tree->SetBranchAddress("fNParameters", &nParBranches); 01379 01380 //assume all events have same number of parameters, so check only first 01381 tree->GetEntry(0); 01382 if(nParBranches != nPar){ 01383 BCLog::OutError(Form("Cannot compute A value: number of parameters in tree (%d)" 01384 "doesn't match number of parameters in model (%d)",nParBranches,nPar)); 01385 return -1.0; 01386 } 01387 01388 01389 //buffer to construct correct branchnames for parameters, e.g. "fParameter3" 01390 char* branchName = new char[10+nPar]; 01391 //set up variables filled for each sample of parameters 01392 //assume same order as in model 01393 for(int i=0;i<(int)nPar;i++){+ 01394 sprintf(branchName, "fParameter%d",i); 01395 tree->SetBranchAddress(branchName, ¶m[i]); 01396 } 01397 01398 // get the p value from Johson's definition for each param from posterior 01399 long nEntries = tree->GetEntries(); 01400 01401 //RN ~ chi2 with K-1 DoF needed for comparison 01402 std::vector<double> randoms(nEntries); 01403 int K = NumberBins(); 01404 BCMath::RandomChi2(randoms,K-1); 01405 01406 //number of Johnson chi2 values bigger than reference 01407 int nBigger=0; 01408 for(int i=0;i<nEntries;i++){ 01409 tree->GetEntry(i); 01410 double chi2 = this->GetChi2Johnson(param); 01411 if(distribution!=0) 01412 distribution->Fill(chi2); 01413 // compare to set of chi2 variables 01414 if(chi2>randoms.at(i) ) 01415 nBigger++; 01416 } 01417 01418 return nBigger/double(nEntries); 01419 }
| double BCModel::GetPvalueFromChi2NDoF | ( | std::vector< double > | par, | |
| int | sigma_index | |||
| ) |
Definition at line 1423 of file BCModel.cxx.
01424 { 01425 double ll = this -> LogLikelihood(par); 01426 int n = this -> GetNDataPoints(); 01427 int npar = this -> GetNParameters(); 01428 01429 double sum_sigma=0; 01430 for (int i=0;i<n;i++) 01431 sum_sigma += log(this -> GetDataPoint(i) -> GetValue(sigma_index)); 01432 01433 double chi2 = -2.*(ll + (double)n/2. * log(2.*M_PI) + sum_sigma); 01434 01435 fChi2NDoF = chi2/double(n-npar); 01436 fPValueNDoF = TMath::Prob(chi2,n-npar); 01437 01438 return fPValueNDoF; 01439 }
| BCH1D * BCModel::CalculatePValue | ( | std::vector< double > | par, | |
| bool | flag_histogram = false | |||
| ) |
Definition at line 1443 of file BCModel.cxx.
01444 { 01445 BCH1D * hist = 0; 01446 01447 // print log 01448 BCLog::OutSummary("Do goodness-of-fit-test"); 01449 01450 // create model test 01451 BCGoFTest * goftest = new BCGoFTest("modeltest"); 01452 01453 // set this model as the model to be tested 01454 goftest -> SetTestModel(this); 01455 01456 // set the point in parameter space which is tested an initialize 01457 // the model testing 01458 if (!goftest -> SetTestPoint(par)) 01459 return 0; 01460 01461 // disable the creation of histograms to save _a lot_ of memory 01462 // (histograms are not needed for p-value calculation) 01463 goftest->MCMCSetFlagFillHistograms(false); 01464 01465 // set parameters of the MCMC for the GoFTest 01466 goftest -> MCMCSetNChains(fGoFNChains); 01467 goftest -> MCMCSetNIterationsMax(fGoFNIterationsMax); 01468 goftest -> MCMCSetNIterationsRun(fGoFNIterationsRun); 01469 01470 // get p-value 01471 fPValue = goftest -> GetCalculatedPValue(flag_histogram); 01472 01473 // get histogram 01474 if (flag_histogram) 01475 { 01476 hist = new BCH1D(); 01477 hist -> SetHistogram(goftest -> GetHistogramLogProb()); 01478 } 01479 01480 // delete model test 01481 delete goftest; 01482 01483 // return histogram 01484 return hist; 01485 }
| double BCModel::GetPValue | ( | ) | [inline] |
| double BCModel::GetPValueNDoF | ( | ) | [inline] |
| double BCModel::GetChi2NDoF | ( | ) | [inline] |
| void BCModel::SetGoFNIterationsMax | ( | int | n | ) | [inline] |
| void BCModel::SetGoFNIterationsRun | ( | int | n | ) | [inline] |
| void BCModel::SetGoFNChains | ( | int | n | ) | [inline] |
| 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 1496 of file BCModel.cxx.
01497 { 01498 // check number of entries in vector 01499 if (point.size() != this -> GetNParameters()) 01500 { 01501 BCLog::OutWarning("BCModel::HessianMatrixElement. Invalid number of entries in the vector."); 01502 return -1; 01503 } 01504 01505 // define steps 01506 double nsteps = 1e5; 01507 double dx1 = par1 -> GetRangeWidth() / nsteps; 01508 double dx2 = par2 -> GetRangeWidth() / nsteps; 01509 01510 // define points at which to evaluate 01511 std::vector<double> xpp = point; 01512 std::vector<double> xpm = point; 01513 std::vector<double> xmp = point; 01514 std::vector<double> xmm = point; 01515 01516 int idx1 = par1 -> GetIndex(); 01517 int idx2 = par2 -> GetIndex(); 01518 01519 xpp[idx1] += dx1; 01520 xpp[idx2] += dx2; 01521 01522 xpm[idx1] += dx1; 01523 xpm[idx2] -= dx2; 01524 01525 xmp[idx1] -= dx1; 01526 xmp[idx2] += dx2; 01527 01528 xmm[idx1] -= dx1; 01529 xmm[idx2] -= dx2; 01530 01531 // calculate probability at these points 01532 double ppp = this -> Likelihood(xpp); 01533 double ppm = this -> Likelihood(xpm); 01534 double pmp = this -> Likelihood(xmp); 01535 double pmm = this -> Likelihood(xmm); 01536 01537 // return derivative 01538 return (ppp + pmm - ppm - pmp) / (4. * dx1 * dx2); 01539 }
| void BCModel::PrintSummary | ( | ) |
Prints a summary on the screen.
Definition at line 1574 of file BCModel.cxx.
01575 { 01576 int nparameters = this -> GetNParameters(); 01577 01578 // model summary 01579 std::cout 01580 << std::endl 01581 << " ---------------------------------" << std::endl 01582 << " Model : " << fName.data() << std::endl 01583 << " ---------------------------------"<< std::endl 01584 << " Index : " << fIndex << std::endl 01585 << " Number of parameters : " << nparameters << std::endl 01586 << std::endl 01587 << " - Parameters : " << std::endl 01588 << std::endl; 01589 01590 // parameter summary 01591 for (int i=0; i<nparameters; i++) 01592 fParameterSet -> at(i) -> PrintSummary(); 01593 01594 // best fit parameters 01595 if (this -> GetBestFitParameters().size() > 0) 01596 { 01597 std::cout 01598 << std::endl 01599 << " - Best fit parameters :" << std::endl 01600 << std::endl; 01601 01602 for (int i=0; i<nparameters; i++) 01603 { 01604 std::cout 01605 << " " << fParameterSet -> at(i) -> GetName().data() 01606 << " = " << this -> GetBestFitParameter(i) 01607 << " (overall)" << std::endl; 01608 if ((int)fBestFitParametersMarginalized.size() == nparameters) 01609 std::cout 01610 << " " << fParameterSet -> at(i) -> GetName().data() 01611 << " = " << this -> GetBestFitParameterMarginalized(i) 01612 << " (marginalized)" << std::endl; 01613 } 01614 } 01615 01616 std::cout << std::endl; 01617 01618 // model testing 01619 if (fPValue >= 0) 01620 { 01621 double likelihood = this -> Likelihood(this -> GetBestFitParameters()); 01622 std::cout 01623 << " - Model testing:" << std::endl 01624 << std::endl 01625 << " p(data|lambda*) = " << likelihood << std::endl 01626 << " p-value = " << fPValue << std::endl 01627 << std::endl; 01628 } 01629 01630 // normalization 01631 if (fNormalization > 0) 01632 std::cout << " Normalization : " << fNormalization << std::endl; 01633 }
| void BCModel::PrintResults | ( | const char * | file | ) |
Prints a summary of the Markov Chain Monte Carlo to a file.
Definition at line 1637 of file BCModel.cxx.
01638 { 01639 // print summary of Markov Chain Monte Carlo 01640 01641 // open file 01642 ofstream ofi(file); 01643 01644 // check if file is open 01645 if(!ofi.is_open()) 01646 { 01647 std::cerr << "Couldn't open file " << file <<std::endl; 01648 return; 01649 } 01650 01651 // number of parameters and chains 01652 int npar = MCMCGetNParameters(); 01653 int nchains = MCMCGetNChains(); 01654 01655 // check convergence 01656 bool flag_conv = MCMCGetNIterationsConvergenceGlobal() > 0; 01657 01658 ofi 01659 << std::endl 01660 << " -----------------------------------------------------" << std::endl 01661 << " Summary" << std::endl 01662 << " -----------------------------------------------------" << std::endl 01663 << std::endl; 01664 01665 ofi 01666 << " Model summary" << std::endl 01667 << " =============" << std::endl 01668 << " Model: " << fName.data() << std::endl 01669 << " Number of parameters: " << npar << std::endl 01670 << " List of Parameters and ranges:" << std::endl; 01671 for (int i = 0; i < npar; ++i) 01672 ofi 01673 << " (" << i << ") Parameter \"" 01674 << fParameterSet -> at(i) -> GetName().data() << "\"" 01675 << ": " << fParameterSet -> at(i) -> GetLowerLimit() 01676 << " - " 01677 << fParameterSet -> at(i) -> GetUpperLimit() << std::endl; 01678 ofi << std::endl; 01679 01680 // give warning if MCMC did not converge 01681 if (!flag_conv && fMCMCFlagRun) 01682 ofi 01683 << " WARNING: the Markov Chain did not converge!" << std::endl 01684 << " Be cautious using the following results!" << std::endl 01685 << std::endl; 01686 01687 // print results of marginalization (if MCMC was run) 01688 if ( fMCMCFlagRun && fMCMCFlagFillHistograms) 01689 { 01690 ofi 01691 << " Results of the marginalization" << std::endl 01692 << " ==============================" << std::endl 01693 << " List of parameters and properties of the marginalized" << std::endl 01694 << " distributions:" << std::endl; 01695 for (int i = 0; i < npar; ++i) 01696 { 01697 if (!fMCMCFlagsFillHistograms.at(i)) 01698 continue; 01699 01700 BCH1D * bch1d = GetMarginalized(fParameterSet -> at(i)); 01701 01702 ofi 01703 << " (" << i << ") Parameter \"" 01704 << fParameterSet->at(i)->GetName().data() << "\"" << std::endl 01705 << " Mean +- RMS: " 01706 << std::setprecision(4) << bch1d->GetMean() 01707 << " +- " 01708 << std::setprecision(4) << bch1d->GetRMS() << std::endl 01709 << " Median +- sigma: " 01710 << std::setprecision(4) << bch1d->GetMedian() 01711 << " + " << std::setprecision(4) << bch1d->GetQuantile(0.84) - bch1d->GetMedian() 01712 << " - " << std::setprecision(4) << bch1d->GetMedian() - bch1d->GetQuantile(0.16) << std::endl 01713 << " (Marginalized) mode: " << bch1d->GetMode() << std::endl 01714 << " Smallest interval(s) containing 68% and local modes:" << std::endl; 01715 01716 std::vector <double> v; 01717 v = bch1d->GetSmallestIntervals(0.68); 01718 int ninter = int(v.size()); 01719 01720 for (int j = 0; j < ninter; j+=5) 01721 ofi 01722 << " " << v.at(j) << " - " << v.at(j+1) 01723 << " (local mode at " << v.at(j+3) 01724 << " with rel. height " << v.at(j+2) 01725 << "; rel. area " << v.at(j+4) << ")" 01726 << std::endl; 01727 01728 ofi 01729 << " 5% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.05) << std::endl 01730 << " 10% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.10) << std::endl 01731 << " 16% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.16) << std::endl 01732 << " 84% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.85) << std::endl 01733 << " 90% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.90) << std::endl 01734 << " 95% quantile: " << std::setprecision(4) << bch1d->GetQuantile(0.95) << std::endl 01735 << std::endl; 01736 } 01737 ofi << std::endl; 01738 } 01739 01740 ofi 01741 << " Results of the optimization" << std::endl 01742 << " ===========================" << std::endl 01743 << " Optimization algorithm used: "; 01744 switch(GetOptimizationMethodMode()) 01745 { 01746 case BCIntegrate::kOptSA: 01747 ofi << " Simulated Annealing" << std::endl; 01748 break; 01749 case BCIntegrate::kOptMinuit: 01750 ofi << " Minuit" << std::endl; 01751 break; 01752 case BCIntegrate::kOptMetropolis: 01753 ofi << " MCMC" << std::endl; 01754 break; 01755 } 01756 01757 if (int(fBestFitParameters.size()) > 0) 01758 { 01759 ofi << " List of parameters and global mode:" << std::endl; 01760 for (int i = 0; i < npar; ++i) 01761 ofi 01762 << " (" << i << ") Parameter \"" 01763 << fParameterSet->at(i)->GetName().data() << "\": " 01764 << fBestFitParameters.at(i) << " +- " << fBestFitParameterErrors.at(i) << std::endl; 01765 ofi << std::endl; 01766 } 01767 else 01768 { 01769 ofi << " No best fit information available." << std::endl; 01770 ofi << std::endl; 01771 } 01772 01773 if (fPValue >= 0.) 01774 { 01775 ofi 01776 << " Results of the model test" << std::endl 01777 << " =========================" << std::endl 01778 << " p-value at global mode: " << fPValue << std::endl 01779 << std::endl; 01780 } 01781 01782 if ( fMCMCFlagRun ) 01783 { 01784 ofi 01785 << " Status of the MCMC" << std::endl 01786 << " ==================" << std::endl 01787 << " Convergence reached: " << (flag_conv ? "yes" : "no") << std::endl; 01788 01789 if (flag_conv) 01790 ofi << " Number of iterations until convergence: " << MCMCGetNIterationsConvergenceGlobal() << std::endl; 01791 else 01792 ofi 01793 << " WARNING: the Markov Chain did not converge! Be\n" 01794 << " cautious using the following results!" << std::endl 01795 << std::endl; 01796 ofi 01797 << " Number of chains: " << MCMCGetNChains() << std::endl 01798 << " Number of iterations per chain: " << MCMCGetNIterationsRun() << std::endl 01799 << " Average efficiencies:" << std::endl; 01800 01801 std::vector <double> efficiencies; 01802 efficiencies.assign(npar, 0.); 01803 01804 for (int ipar = 0; ipar < npar; ++ipar) 01805 for (int ichain = 0; ichain < nchains; ++ichain) 01806 { 01807 int index = ichain * npar + ipar; 01808 efficiencies[ipar] += double(MCMCGetNTrialsTrue().at(index)) / double(MCMCGetNTrialsTrue().at(index) + MCMCGetNTrialsFalse().at(index)) / double(nchains) * 100.; 01809 } 01810 01811 for (int ipar = 0; ipar < npar; ++ipar) 01812 ofi 01813 << " (" << ipar << ") Parameter \"" 01814 << fParameterSet->at(ipar)->GetName().data() << "\": " 01815 << efficiencies.at(ipar) << "%" << std::endl; 01816 ofi << std::endl; 01817 } 01818 01819 ofi 01820 << " -----------------------------------------------------" << std::endl 01821 << std::endl 01822 << " Notes" << std::endl 01823 << " =====" << std::endl 01824 << " (i) Median +- sigma denotes the median, m, of the" << std::endl 01825 << " marginalized distribution and the intervals from" << std::endl 01826 << " m to the 16% and 84% quantiles." << std::endl 01827 << " -----------------------------------------------------" << std::endl; 01828 01829 // close file 01830 // ofi.close; 01831 }
| void BCModel::PrintShortFitSummary | ( | int | chi2flag = 0 |
) |
Prints a short summary of the fit results on the screen.
Definition at line 1835 of file BCModel.cxx.
01836 { 01837 BCLog::OutSummary("---------------------------------------------------"); 01838 BCLog::OutSummary(Form("Fit summary for model \'%s\':", GetName().data())); 01839 BCLog::OutSummary(Form(" Number of parameters: Npar = %i", GetNParameters())); 01840 BCLog::OutSummary(Form(" Number of data points: Ndata = %i", GetNDataPoints())); 01841 BCLog::OutSummary(" Number of degrees of freedom:"); 01842 BCLog::OutSummary(Form(" NDoF = Ndata - Npar = %i", GetNDataPoints()-GetNParameters())); 01843 01844 BCLog::OutSummary(" Best fit parameters (global):"); 01845 for (unsigned int i = 0; i < GetNParameters(); ++i) 01846 BCLog::OutSummary(Form(" %s = %.3g", GetParameter(i)->GetName().data(), GetBestFitParameter(i))); 01847 01848 BCLog::OutSummary(" Goodness-of-fit test:"); 01849 BCLog::OutSummary(Form(" p-value = %.3g", GetPValue())); 01850 if(chi2flag) 01851 { 01852 BCLog::OutSummary(Form(" p-value corrected for NDoF = %.3g", GetPValueNDoF())); 01853 BCLog::OutSummary(Form(" chi2 / NDoF = %.3g", GetChi2NDoF())); 01854 } 01855 BCLog::OutSummary("---------------------------------------------------"); 01856 }
| 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 1860 of file BCModel.cxx.
01861 { 01862 // check number of entries in vector 01863 if (parameters.size() != GetNParameters()) 01864 { 01865 BCLog::OutError("BCModel::PrintHessianMatrix : Invalid number of entries in the vector"); 01866 return; 01867 } 01868 01869 // print to screen 01870 std::cout 01871 << std::endl 01872 << " Hessian matrix elements: " << std::endl 01873 << " Point: "; 01874 01875 for (int i = 0; i < int(parameters.size()); i++) 01876 std::cout << parameters.at(i) << " "; 01877 std::cout << std::endl; 01878 01879 // loop over all parameter pairs 01880 for (unsigned int i = 0; i < GetNParameters(); i++) 01881 for (unsigned int j = 0; j < i; j++) 01882 { 01883 // calculate Hessian matrix element 01884 double hessianmatrixelement = HessianMatrixElement(fParameterSet->at(i), 01885 fParameterSet->at(j), parameters); 01886 01887 // print to screen 01888 std::cout << " " << i << " " << j << " : " << hessianmatrixelement << std::endl; 01889 } 01890 }
| void BCModel::FixDataAxis | ( | unsigned int | index, | |
| bool | fixed | |||
| ) |
Definition at line 1543 of file BCModel.cxx.
01544 { 01545 // check if index is within range 01546 if (index < 0 || index > fDataSet -> GetDataPoint(0) -> GetNValues()) 01547 { 01548 BCLog::OutWarning("BCModel::FixDataAxis. Index out of range."); 01549 return; 01550 } 01551 01552 if (fDataFixedValues.size() == 0) 01553 fDataFixedValues.assign(fDataSet -> GetDataPoint(0) -> GetNValues(), false); 01554 01555 fDataFixedValues[index] = fixed; 01556 }
| virtual double BCModel::CDF | ( | const std::vector< double > & | parameters, | |
| int | index, | |||
| bool | lower = false | |||
| ) | [inline, virtual] |
1dim cumulative distribution function of the probability to get the data f(x|param) for a single measurement, assumed to be identical for all measurements
| parameters | The parameter values at which point to compute the cdf | |
| index | The data point index starting at 0,1...N-1 | |
| lower | only needed for discrete distributions! Return the CDF for the count one less than actually observed, e.g. in Poisson process, if 3 actually observed, then CDF(2) is returned |
Reimplemented in BCHistogramFitter.
Definition at line 626 of file BCModel.h.
| BCDataPoint * BCModel::VectorToDataPoint | ( | std::vector< double > | data | ) | [private] |
Converts a vector of doubles into a BCDataPoint
Definition at line 1894 of file BCModel.cxx.
01895 { 01896 int sizeofvector = int(data.size()); 01897 BCDataPoint * datapoint = new BCDataPoint(sizeofvector); 01898 datapoint -> SetValues(data); 01899 return datapoint; 01900 }
| int BCModel::CompareStrings | ( | const char * | string1, | |
| const char * | string2 | |||
| ) | [private] |
Compares to strings
Definition at line 1904 of file BCModel.cxx.
01905 { 01906 int flag_same = 0; 01907 01908 if (strlen(string1) != strlen(string2)) 01909 return -1; 01910 01911 for (int i = 0; i < int(strlen(string1)); i++) 01912 if (string1[i] != string2[i]) 01913 flag_same = -1; 01914 01915 return flag_same; 01916 }
| int BCModel::NumberBins | ( | ) | [inline, private] |
rule of thumb for good number of bins (Wald1942, Johnson2004) to group observations updated so minimum is three bins (for 1-5 observations)!
Definition at line 715 of file BCModel.h.
00715 { 00716 return (int)(exp(0.4 * log(this -> GetNDataPoints())) + 2); 00717 }
int BCModel::fIndex [protected] |
std::string BCModel::fName [protected] |
double BCModel::fModelAPriori [protected] |
double BCModel::fModelAPosteriori [protected] |
BCParameterSet* BCModel::fParameterSet [protected] |
BCDataSet* BCModel::fDataSet [protected] |
unsigned int BCModel::fNDataPointsMinimum [protected] |
unsigned int BCModel::fNDataPointsMaximum [protected] |
bool BCModel::flag_ConditionalProbabilityEntry [protected] |
double BCModel::fPValue [protected] |
double BCModel::fChi2NDoF [protected] |
double BCModel::fPValueNDoF [protected] |
bool BCModel::flag_discrete [protected] |
int BCModel::fGoFNIterationsMax [protected] |
int BCModel::fGoFNIterationsRun [protected] |
int BCModel::fGoFNChains [protected] |
double BCModel::fNormalization [private] |
1.5.1